osr和Pyproj库的简单使用

osr和pyproj库的使用

坐标系统的构成

无论是矢量数据处理还是栅格数据的处理都需要用到坐标系统，坐标系统的构成主要是三部分，一、椭球体，可以用两个非常简单的参数来表示，a长半轴长度，b短半轴长度（也可以表述为长半轴a，短半轴b，扁率e三者知道任意两个也可表述该坐标系统），二、坐标系，为了描述地球上某一点的位置，人们想出了不同的坐标系统，如非常有名的空间直角坐标系，三、投影，为了制图以及量测的便利，人们想出了不同的投影办法去将椭球体展开到平面上。

使用osr对整个图层进行投影变换

思路：创建基于目标投影的图层，然后将原图层的几何要素一一投影过去，写入到目标图层中。

'''
Created on 2020年2月9日

@author: Sun Strong
'''
from osgeo import ogr,osr,gdal
gdal.SetConfigOption('OGR_ENABLE_PARTIAL_REPROJECTION','TRUE')#部分要素投影失败仍然正常投影
#定义原坐标系统和目标坐标系统
web_sr=osr.SpatialReference()
web_sr.ImportFromProj4('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs')
peters_sr=osr.SpatialReference()
peters_sr.ImportFromProj4('+proj=cea +ellps=WGS84 +datum=WGS84 +units=m +')
path=r"C:\Users\Sun Strong\Desktop\test"
old_lyrname='ne_110m_land_1p'#旧图层名
new_lyrname='ne_110m_land_1p'+'_peters'#新图层名
def TransformWholeLayer(sr1,sr2,path,old_lyrname,new_lyrname):
    ct=osr.CoordinateTransformation(sr1,sr2)#tranform的接收参数
    ds=ogr.Open(path,1)
    lyr=ds.GetLayer(old_lyrname)
    if ds.GetLayer(new_lyrname):
        ds.DeleteLayer(new_lyrname)
    peter_lyr=ds.CreateLayer(new_lyrname,sr2,ogr.wkbPolygon)
    for feat in lyr:#对图层的每个几何对象进行坐标转换，并写入到新图层
        p_feature=ogr.Feature(peter_lyr.GetLayerDefn())
        geom=feat.geometry().Clone()
        geom.Transform(ct)
        p_feature.SetGeometry(geom)
        peter_lyr.CreateFeature(p_feature)
    del ds
    print("This layer has been transformed successfully!")
TransformWholeLayer(web_sr,peters_sr,path,old_lyrname,new_lyrname)

使用pyproj进行单点坐标的转换

'''
Created on 2020年2月9日

@author: Sun Strong
'''

import pyproj
#地理坐标和投影坐标之间的转换
'''
utm_proj=pyproj.Proj('+proj=utm +zone=31 +ellps=WGS84')
utm_proj=pyproj.Proj(proj='utm',zone=31,ellps='WGS84')
'''
utm_proj=pyproj.Proj(init='epsg:32631')
x,y=utm_proj(2.294694,48.858093)
print(x,y)
a,b=utm_proj(x,y,inverse=True)
print(a,b)
#第二种转换方法，支持对投影坐标以及投影坐标之间的互相转换
wgs84=pyproj.Proj('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ')
proj1=pyproj.Proj('+proj=utm +zone=31 +ellps=WGS84')
x1,y1=pyproj.transform(proj1,wgs84,580744.32,4504695.26)
print(x1,y1)
x2,y2=utm_proj(580744.32,4504695.26,inverse=True)
print(x2,y2)

代码中共提供了两种坐标转换方法，其中第一种**仅仅限于投影坐标和其对应的地理投影坐标（相同基准）之间的转换，第二种转换方法则通用，投影坐标和投影坐标（不同基准）**也可以实现该转换。