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

添加新云天气象小编QQ130188121，更快速度获取气象升学、就业、会议、征稿及学术动态信息！

最新热点文章：南京信息工程大学2020年招收攻读硕士学位研究生专业目录

学气象都知道NOAA，那COAA是啥？第八届COAA国际会议将在南信大气象楼召开啦！

“宇宙最强”科学联盟成立，国科大、北大、南信大等27所高校加入！
[h1][/h1][h1]python与地理空间分析（一）中简单介绍了地理空间分析对于数据分析和气象的重要作用，包含常用到的GIS数据类型和处理的python包的介绍，本篇文章书接上文，将对在GIS中常打交道的矢量数据的处理做简单介绍。考虑到在日常中对GIS常用的功能，包含的主题包括一些应用和数据的介绍：[/h1]1 GIS中距离的计算2 方位计算3 坐标转换和投影转换
[h1]距离测量[/h1]

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

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

#以横轴墨卡托投影下的地图坐标，单位为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，这是两者的坐标和投影不一样，前者经过墨卡托投影将地球展成了平面。
[h2]半正矢公式[/h2]

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

#点的坐标是经纬度

复制代码

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，比单纯利用勾股定理计算，准确度大大提高。半正矢公式是最常用的距离计算公式，在一定精度保证条件下，代码简便。
[h2]Vincenty公式[/h2]

大家学习地理时，都知道地球并不是标准的球形，因此单纯将地球简化为球形，来计算距离，也会存在误差。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算法实现。
[h1]方位计算[/h1]

气象中常用来计算风向，即风的方位，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

复制代码

[h1]坐标转换和重投影[/h1]

地理空间分析中，绕不开坐标和投影，在数据处理中，可能不同的数据源有着不同的坐标投影，这就需要把它们统一起来进行转换，然后再分析。
[h2]坐标转换[/h2]在气象数据中，常用到的投影是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')

复制代码

[h2]重投影[/h2]当对一个矢量文件或者栅格数据进行坐标统一时，就需要进行重投影，进行不同坐标系统的转换，例如将经常用的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

[h1]总结[/h1]本次文件介绍了，地理空间分析中对矢量数据一些应用算法的介绍，下次的主题是对矢量数据（主要是shapefile格式文件）的处理
往期文章
python与地理空间分析（一）
还在用matplotlib画图？你out啦
如何使用Python处理shp文件如何使用Python处理HDF格式数据

往期回顾：
南信大2020招收硕士研究生《天气学(含天气分析)》考试大纲
夏日为何汽车成蒸笼？桑拿天为何总不走？南信大召开COAA会议，气象专家如是说

南京信息工程大学2020年招收攻读硕士学位研究生专业目录

学气象都知道NOAA，那COAA是啥？第八届COAA国际会议将在南信大气象楼召开啦！

2020大气科学类专业校招开始早！附已开始招聘单位简章汇总

2020校招第4波！中国电科第十四研究所招聘大气科学类启动！

“宇宙最强”科学联盟成立，国科大、北大、南信大等27所高校加入！

南京信息工程大学2020年招收博士研究生专业一览表

南京信息工程大学2020年招收硕士研究生专业一览表

“别摆博士的谱，刚进中央气象台，你就是个生瓜蛋子” | 老戚外传

2020校招新动态！| 兰州大学大气院走访气象局推进2020届气象专场招聘会
中国内地高校在2019世界一流学科(大气科学)排名中的表现

大气科学经典教科书《动力气象学引论》中译本出版

视频 | 南信大校长李北群：争做大气科学“世界单打冠军”

福州气象 | 第七届海青节·两岸青年生态旅游气象服务创新交流会等你来！

视频 | 盘点全国高校的海洋与气象学院(系)

北清华，南信大！南信大获2019世界一流学科(大气科学)排名第18名

高考志愿，到底该不该报气象专业？| 气象职场
欢迎添加"新云天气象"小编的QQ、微信和QQ群！
深圳远征技术有限公司招聘防雷及相关专业人才！
请及时加入气象升学就业QQ群并添加新云天气象小编个人微信

南信大《天气学和天气分析》《数学（理）》考研自命题科目远程高端小班招生简章
新云天气象2020考研讲堂PPT
"新云天气象"小编的QQ空间开通了！

欢迎各用人单位、研究生培养单位通过“新云天气象”微信公众号、QQ群发布招聘、招生信息！
用人（招生）单位：长按（扫描）下方二维码，加小编个人微信或QQ号，免费提交贵单位的招聘、招生信息！
应聘学生（考生）：长按（扫描）下方二维码，加小编个人微信或QQ号，免费获得招聘、招生信息！

QQ号：130188121
QQ空间网址：
https://130188121.qzone.qq.com/
友情提示：添加时，务必说明您是哪个单位，哪个部门，以及怎么称呼。
学生务必说明学校、专业年级和姓名。谢谢您的配合！