python地理数据处理-python与地理空间分析（二）矢量数据

python地理数据处理_python与地理空间分析（⼆）⽮量数据python与地理空间分析（⼀）中简单介绍了地理空间分析对于数据分析和⽓象的重要作⽤，包含常⽤到的GIS数据类型和处理的python包的介绍，本篇⽂章书接上⽂，将对在GIS中常打交道的⽮量数据的处理做简单介绍。考虑到在⽇常中对GIS常⽤的功能，包含的主题包括⼀些应⽤和数据的介绍：·1GIS中距离的计算·2⽅位计算·3坐标转换和投影转换距离测量距离测量是地理空间分析中的⼀个⾮常重要的功能，在⽓象数据处理中也会经常⽤到，例如查找最临近的⽓象站点、⽓象站点数据与其他数据匹配等操作。⽬前，针对不同的地球模型，计算地球上两点的距离，有三种不同的算法：勾股定理把地球当作⼀个没有曲率的平⾯模型，计算两点的距离即计算直线的距离，根据坐标利⽤勾股定理就可以计算，但是地球本⾝是具有曲率的，勾股定理的计算，⽐较简单和快速，在尺度上可以得到⼀个在可接受误差范围的距离，对精度有⼀定要求的并不能满⾜。#以横轴墨卡托投影下的地图坐标，单位为mimportmathx1=456456.23y1=1279721.064x2=576628.34y2=1071740.33x_dist=x1-x2y_dist=y1-y2dist_sq=x_dist**2+y_dist**2distance=math.sqrt(dist_sq)计算的结果应该是240202.668047278m#直接利⽤经纬度计算相同的点，需要转⼀下弧度importmathlon1=-90.21lat1=32.31lon2=-88.95lat2=30.43lon_dist=math.radians(lon1-lon2)lat_dist=math.radians(lat1-lat2)dist_sq=lon_dist**2+lat_dist**2distance=math.sqrt(dist_sq)*6371251.46#地球半径，转换单位计算的结果251664.4668769659m，两个相同的点，计算的结果会相差⼤约10km，这是两者的坐标和投影不⼀样，前者经过墨卡托投影将地球展成了平⾯。半正⽮公式

两点之间直线最短，但我们通过地图在看飞机航线时，往往并不是直线，⽽是成弧状，这就让我们⾮常疑惑。其实我们理解的直线就是利⽤勾股定理计算的地图上的两点间的距离，⽽实际中的航线是计算的⼤圆距离也是球⾯距离。⼤圆距离是指球体把桌⾯上两点之间的距离，球⾯上任意两点以及球⼼可以确定唯⼀的⼤圆，在这个⼤圆上连接这两点的较短的弧的长度就是⼤圆距离。计算⼤圆距离常⽤的算法就是半正⽮公式。#点的坐标是经纬度importmathlon1=-90.21lat1=32.31lon2=-88.95lat2=30.43lon_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公式就是基于椭球体地球模型的计算距离的公式。但是公式更复杂，且需要选择贴合本地的椭球模型参数。importmathdistance=Nonelon1=-90.21lat1=32.31lon2=-88.95lat2=30.43#椭球体参数，采⽤NAD83a=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=Lforiinrange(100):sinLam=math.sin(lam)cosLam=math.cos(lam)sinSigma=math.sqrt((cosU2*sinLam)**2+(cosU1*sinU2-sinU1*cosU2*cosLam)**2)ifsinSigma==0:distance=0#重合点breakcosSigma=sinU1*sinU2+cosU1*cosU2*cosLamsigma=math.atan2(sinSigma,cosSigma)sinAlpha=cosU1*cosU2*sinLam/sinSigmacosSqAlpha=1-sinAlpha**2cos2SigmaM=cosSigma-2*sinU1*sinU2/cosSqAlphaifmath.isnan(cos2SigmaM):cosSigmaM=0#⾚道线C=f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))LP=lamlam=L+(1-C)*f*sinAlpha*(sigma+C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))ifnotabs(lam-LP)>1e-12:breakuSq=cosSqAlpha*(a**2-b**2)/b**2A=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中有时也需要提供两点之间的朝向⽅位。frommathimportatan2,cos,sin,degreeslon1=-90.21

lat1=32.31lon2=-88.95lat2=30.43angle=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计算风向的代码importmathangleDegree=math.atan2(u,v)*180/math.piwin_Degree=angleDegree+180ifwin_Degree<0:win_Degree+=360elifwin_Degree==360:win_Degree=0#返回16⽅位val=int((win_deg/22.5)+.5)win_16=val%16坐标转换和重投影地理空间分析中，绕不开坐标和投影，在数据处理中，可能不同的数据源有着不同的坐标投影，这就需要把它们统⼀起来进⾏转换，然后再分析。坐标转换在⽓象数据中，常⽤到的投影是UTM投影，且⼀般是等距离投影，⽽⼀些数据中为了⽅便计算，常⽤等经纬度的投影，这就需要坐标之间的转换。常⽤的⼯具包是utm包：需要根据所在的投影带（上图）进⾏查找获取相应参数importutm#需要安装x=3578217.8414962334y=762684.723145958zone=15band='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的pythonAPI的帮助，也是GDAL的⼀部分。下⾯是⼀个简单⽰例，将⼀个shapfile⽂件进⾏重投影操作。#需要安装GDAL包fromosgeoimportogr

fromosgeoimportosrimportosimportshutil#简单的copy和重命名srcName="NYC_MUSEUMS_LAMBERT.shp"tgtName="NYC_MUSEUMS_GEO.shp"tgt_spatRef=osr.SpatialReference()tgt_spatRef.ImportFromEPSG(4326)driver=ogr.GetDriverByName("ESRIShapefile")src=driver.Open(srcName,0)srcLyr=src.GetLayer()src_spatRef=srcLyr.GetSpatialRef()ifos.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()whilesrcFeat: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.c