python地理数据处理

python与地理空间分析（一）中简单介绍了地理空间分析对于数据分析和气象的重要作用，包含常用到的GIS数据类型和处理的python包的介绍，本篇文章书接上文，将对在GIS中常打交道的矢量数据的处理做简单介绍。考虑到在日常中对GIS常用的功能，包含的主题包括一些应用和数据的介绍：

•1 GIS中距离的计算•2 方位计算•3 坐标转换和投影转换

距离测量

距离测量是地理空间分析中的一个非常重要的功能，在气象数据处理中也会经常用到，例如查找最临近的气象站点、气象站点数据与其他数据匹配等操作。目前，针对不同的地球模型，计算地球上两点的距离，有三种不同的算法：

勾股定理

把地球当作一个没有曲率的平面模型，计算两点的距离即计算直线的距离，根据坐标利用勾股定理就可以计算，但是地球本身是具有曲率的，勾股定理的计算，比较简单和快速，在尺度上可以得到一个在可接受误差范围的距离，对精度有一定要求的并不能满足。#以横轴墨卡托投影下的地图坐标，单位为m

import math

x1=456456.23

y1=1279721.064

x2=576628.34

y2=1071740.33

x_dist=x1-x2

y_dist=y1-y2

dist_sq=x_dist**2+y_dist**2

distance=math.sqrt(dist_sq)

计算的结果应该是240202.668047278m#直接利用经纬度计算相同的点，需要转一下弧度

import math

lon1=-90.21

lat1=32.31

lon2=-88.95

lat2=30.43

lon_dist=math.radians(lon1-lon2)

lat_dist=math.radians(lat1-lat2)

dist_sq=lon_dist**2+lat_dist**2

distance=math.sqrt(dist_sq)*6371251.46 #地球半径，转换单位

计算的结果251664.4668769659m，两个相同的点，计算的结果会相差大约10km，这是两者的坐标和投影不一样，前者经过墨卡托投影将地球展成了平面。

半正矢公式

两点之间直线最短，但我们通过地图在看飞机航线时，往往并不是直线，而是成弧状，这就让我们非常疑惑。其实我们理解的直线就是利用勾股定理计算的地图上的两点间的距离，而实际中的航线是计算的大圆距离也是球面距离。大圆距离是指球体把桌面上两点之间的距离，球面上任意两点以及球心可以确定唯一的大圆，在这个大圆上连接这两点的较短的弧的长度就是大圆距离。计算大圆距离常用的算法就是半正矢公式。#点的坐标是经纬度

import math

lon1=-90.21

lat1=32.31

lon2=-88.95

lat2=30.43

lon_dist=math.radians(lon1-lon2)

lat_dist=math.radians(lat1-lat2)

lat1_rad=math.radians(lat1)

lat2_rad=math.radians(lat2)

a=math.sin(lat_dist/2)**2+math.sin(lon_dist/2)**2*math.cos(lat1_rad)*math.cos(lat2_rad)

c=2*math.asin(math.sqrt(a))

distance=c*6371251.46 #地球半径，转换单位

计算的结果240857.59366623117m，与经过墨卡托投影后利用计算的结果（240202.668047278m）相差0.5km，比单纯利用勾股定理计算，准确度大大提高。半正矢公式是最常用的距离计算公式，在一定精度保证条件下，代码简便。

Vincenty公式

大家学习地理时，都知道地球并不是标准的球形，因此单纯将地球简化为球形，来计算距离，也会存在误差。Vincenty公式就是基于椭球体地球模型的计算距离的公式。但是公式更复杂，且需要选择贴合本地的椭球模型参数。import math

distance=None

lon1=-90.21

lat1=32.31

lon2=-88.95

lat2=30.43

#椭球体参数，采用NAD83

a=6378137 #半长轴

f=1/298.257222101 #逆扁平率

b=abs((f*a)-a) #半短轴

L=math.radians(lon1-lon2)

U1=math.atan((1-f)*math.tan(math.radians(lat1)

U2=U1=math.atan((1-f)*math.tan(math.radians(lat2)

sinU1=math.sin(U1)

cosU1=math.cos(U1)

sinU2=math.sin(U2)

cosU2=math.cos(U2)

lam=L

for i in range(100):

sinLam=math.sin(lam)

cosLam=math.cos(lam)

sinSigma=math.sqrt((cosU2*sinLam)**2+(cosU1*sinU2-sinU1*cosU2*cosLam)**2)

if sinSigma==0:

distance=0 #重合点

break

cosSigma=sinU1*sinU2+cosU1*cosU2*cosLam

sigma=math.atan2(sinSigma,cosSigma)

sinAlpha=cosU1*cosU2*sinLam/sinSigma

cosSqAlpha=1-sinAlpha**2

cos2SigmaM=cosSigma-2*sinU1*sinU2/cosSqAlpha

if math.isnan(cos2SigmaM):

cosSigmaM=0 #赤道线

C=f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))

LP=lam

lam=L+(1-C)*f*sinAlpha*(sigma+C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))

if not abs(lam-LP) > 1e-12:

break

uSq=cosSqAlpha*(a**2-b**2)/b**2

A=1+uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))

B=uSq/1024*(256+uSq*(-128+uSq*(74-47*uSq)))

deltaSigma=B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))

s=b*A*(sigma-deltaSigma)

distance=s

计算结果240448.59665784411m，与投影变换的只相差200多米，虽然公式比较复杂，但更精准。总体来看，如果对于精度需求不是很高，而且数据量合适的情况下，可以选择半正矢公式算法，在保证一定精度条件下，高效率完成计算；如果计算量特别大，一种是将数据转换到投影坐标下，利用勾股定理进行计算，对精度没有太大要求，也可以直接利用勾股定理；在两极地区或者对精度要求非常高，可以采用Vincerty算法实现。

方位计算

气象中常用来计算风向，即风的方位，GIS中有时也需要提供两点之间的朝向方位。from math import atan2,cos,sin,degrees

lon1=-90.21

lat1=32.31

lon2=-88.95

lat2=30.43

angle=atan2(cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1),sin(lon2-lon1)*cos(lat2))

bearing=(degrees(angle)+360)%360 #避免负值

计算结果308.7992752836875 此外，补充利用风速U，V计算风向的代码import math

angleDegree = math.atan2(u, v) * 180 / math.pi

win_Degree = angleDegree + 180

if win_Degree < 0:

win_Degree += 360

elif win_Degree == 360:

win_Degree = 0

#返回16方位

val = int((win_deg / 22.5) + .5)

win_16=val % 16

坐标转换和重投影

地理空间分析中，绕不开坐标和投影，在数据处理中，可能不同的数据源有着不同的坐标投影，这就需要把它们统一起来进行转换，然后再分析。

坐标转换

在气象数据中，常用到的投影是UTM投影，且一般是等距离投影，而一些数据中为了方便计算，常用等经纬度的投影，这就需要坐标之间的转换。常用的工具包是utm包：

需要根据所在的投影带（上图）进行查找获取相应参数import utm #需要安装

x=3578217.8414962334

y=762684.723145958

zone=15

band='S'

utm.to_latlon(y,x,zone,band)

>> (32.30999990799516, -90.20999715032659)

utm.from_latlon(32.31,-90.21)

(762684.723145958, 3578217.8414962334, 15, 'S')

重投影

当对一个矢量文件或者栅格数据进行坐标统一时，就需要进行重投影，进行不同坐标系统的转换，例如将经常用的4326墨卡托投影转为互联网地图中的3857web墨卡托投影。重投影需要依靠OGR的python API的帮助，也是GDAL的一部分。下面是一个简单示例，将一个shapfile文件进行重投影操作。#需要安装GDAL包

from osgeo import ogr

from osgeo import osr

import os

import shutil #简单的copy和重命名

srcName="NYC_MUSEUMS_LAMBERT.shp"

tgtName="NYC_MUSEUMS_GEO.shp"

tgt_spatRef=osr.SpatialReference()

tgt_spatRef.ImportFromEPSG(4326)

driver=ogr.GetDriverByName("ESRI Shapefile")

src=driver.Open(srcName,0)

srcLyr=src.GetLayer()

src_spatRef=srcLyr.GetSpatialRef()

if os.path.exists(tgtName):

driver.DeleteDataSource(tgtName)

tgt=driver.CreateDataSource(tgtName)

lyrName=os.path.splitext(tgtName)[0]

#使用wkb格式声明几何图形

tgtLyr=tgt.CreateLayer(lyrName,geom_type=ogr.wkbPoint)

featDef=srcLyr.GetLayerDefn()

trans=osr.CoordinateTransformation(src_spatRef,tgt_spatRef)

srcFeat=srcLyr.GetNextFeature()

while srcFeat:

geom=srcFeat.GetGeometryRef()

geom.Transform(trans)

feature=ogr.Feature(featDef)

feature.SetGeometry(geom)

tgtLyr.CreateFeature(feature)

feature.Destroy()

srcFeat.Destroy()

srcFeat=srcLyr.GetNextFeature()

src.Destroy()

tgt.Destroy()

tgt_spatRef.MorphToESRI()

prj=open(lyrname+".prj",'w')

prj.write(tgt_spatRef.ExportToWkt())

prj.close()

srcDbf=os.path.splitext(srcName)[0]+'.dbf'

tgtDbf=lyrName+'.dbf'

shutil.copyfile(srcDbf,tgtDbf)

数据链接:https://pan.baidu.com/s/1kg5H-6P91xywcGG_innBhA 密码:m5ga。

google earth展示的结果：

通过查看shp文件的信息：

ogrinfo -al -so NYC_MUSEUMS_LAMBERT.shp

ogrinfo -al -so NYC_MUSEUMS_GEO.shp

总结

本次文件介绍了，地理空间分析中对矢量数据一些应用算法的介绍，下次的主题是对矢量数据（主要是shapefile格式文件）的处理