GDAL笔记--chapter8
#坐标系主要有大地坐标系、地理坐标系和投影坐标系 #其中地理坐标系使用三维球面定义地球表面位置，以实现经纬度对点位的应用。（椭球体参数和基准面）,可以叫做地理投影或者未投影 #投影坐标系是平面坐标，使用X、Y值的坐标系统描述点的位置。它由基准面和投影方法确定。 #等角、等积、等方位投影... from osgeo import ogr ds = ogr.Open('') sr = ds.GetLayer().GetSpatialRef()
使用空间参照系统(OSR/PROJ4)
1.从SR中获取信息
#从SR中获取信息 sr.GetAttrValue('PROJCS') #获取投影名称 sr.GetAttrValue('AUTHORITY') sr.GetAttrValue('AUTHORITY', 1) #获取SRID指定子级的参数 sr.GetAuthorityCode('DATUM') #获取指定SRID的代码，如基准面 sr.GetProjParm('osr.SRS_PP_FALSE_EASTING')
2.创建空间参考对象
#创建空间参考对象(EPSG/PROJ.4/WKT) from osgeo import osr sr = osr.SpatialReference() sr.ImportFromEPSG(26912) #NAD83 UTM 12N sr.GetAttrValue('PROJCS') #获取投影验证 sr.ImportFromProj4('''+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs''') sr = osr.SpatialReference(wkt) #或者ImportFromWkt,wkt字符串太长，这里省略 #自己构造空间参考对象 sr = osr.SpatialReference() sr.SetProjCS('USGS Albers') #命名 sr.SetWellKnownGeogCS('NAD83') #指定基准面 sr.SetACEA(29.5, 45.5, 23, -96, 0, 0) #设定阿尔伯斯等面积投影参数，标准平行线1/2、中央纬线/经线、东移、北移 sr.Fixup() #添加默认值 sr.Validate()
3.为数据分配srs
#为数据分配srs sr = osr.SpatialReferrence() sr.ImportFromEPSG(26912) driver = ogr.GetDriverByName('ESRI Shapefile') ds = driver.CreateDataSource(r'') lyr = ds.CreateLayer('countries', sr, ogr.wkbPolygon) #这里是先向图层分配srs，只是提供了相关信息。下面需要用UTM坐标系创建几何对象 geom.AssignSpatialReference(sr) #如果用的是单个几何对象而不是图层，那么需要为几何对象分配srs
4.几何对象重投影（OSR库）
（这里用OSR库对几何对象重投影，用于shapefiles文件多，栅格坐标也可以处理。与GDAL库中的gdal.Transformer、TransformPoint不同，后者输入为数据源而前者为srs；后者用于栅格坐标多而前者用于几何对象多）
Ps：还有个图像坐标和投影坐标的转换，和这里不同。那个需要用到GetGeoTransform，属于仿射变换。后面将会对这些进行总结。
#几何对象重投影 world = pb.get_shp_geom(r'...') #创建几何对象 tower = ogr.Geometry(wkt='POINT (2.294 48.858') tower.AssignSpatialReference(osr.SpatialReference(osr.SRS_WKT_WGS84)) from osgeo import gdal gdal.SetConfigOption('OGR_ENABLE_PARTIAL_PROJECTION', 'TRUE') #跳过不能成功的投影点 web_mercator_sr = osr.SpatialReference() #创建新srs web_mercator_sr.ImportFromEPSG(3857) #对已分配srs的几何对象，TransformTo即可(即已经有srs信息） world.TransformTo(web_mercator_sr) tower.TransformTo(web_mercator_sr) #对未分配srs的几何对象，但已知其srs的情况下 peters_sr = osr.SpatialReference() peters_sr.ImportFromProj4('...') #目标srs，而web_mercator_sr是已知的srs ct = osr.CoordinateTransformation(web_mercator_sr, peters_sr) world.Transform(ct)
5.更改基准
#更改基准:1.用网格位移文件 2.查询基准转换函数及值（如下） #注：必须确保原始和目标空间参考都有基准信息 sr = osr.SpatialReference() sr.SetWellKnownGeogCS('NAD27') sr.SetTOWGS84(-8, 160, 176)
6.重投影一个点图层
#重投影一个点图层 from osgeo import ogr,osr sr = osr.SpatialReference() #创建目标投影 sr.ImportFromProj4('...') ds = ogr.Open(r'...') in_lyr = ds.GetLayer('us_volcanos') out_lyr = ds.CreateLayer('us_volcanos_aea', sr, ogr.wkbPoint) #创建目标图层 out_lyr.CreateFields(in_lyr.schema) #复制图层属性列表 out_feat = ogr.Feature(out_lyr.GetLayerDefn()) #根据图层定义创建要素 for in_feat in in_lyr: geom = in_feat.geometry().Clone() #获取几何体并复制 geom.TransformTo(sr) #几何体投影转换 out_feat.SetGeometry(geom) #设置几何体 for i in range(in_feat.GetFieldCount()): #对每个要素，遍历字段 out_feat.SetField(i, in_feat.GetField(i)) out_lyr.CreateFeature(out_feat) #图层创建要素
7.pyproj模块在处理坐标值列表方面有优势，有两种方法转换坐标
#pyproj模块在处理坐标值列表方面有优势，有两种方法转换坐标 #1.pyproj类。示例为地理坐标和投影坐标的转换 import pyproj utm_proj = pyproj.Proj('+proj=utm +zone=31 +ellps=WGS84') #定义投影 x, y = utm_proj(2.29, 48.85) #经纬度转投影坐标 x1, y1 = utm_proj(x, y, inverse=True) #投影坐标转经纬度 #2.transform函数。投影坐标系之间转换，以基准转换为例 wgs84 = pyproj.Proj('+proj=utm +zone=18 +datum=WGS84') nad27 = pyproj.Proj('+proj=utm +zone=18 +datum=NAD27') x, y = pyproj.transform(wgs84, nad27, 580744, 4504695) #初始化Proj对象 p = pyproj.Proj('+proj=utm +zone=31 +ellps=WGS84') p = pyproj.Proj(proj='utm', zone=31, ellps='WGS84') p = pyproj.Proj(init='epsg:32631')
8.在地球上计算大圆距离，由始末点位计算向前、向后方向和距离
#计算大圆距离，由始末点位计算向前、向后方向和距离 la_lat, la_lon = 34.050, -118.250 berlin_lat, berlin_lon = 52.516, 13.383 geod = pyproj.Geod(ellps='WGS84') forward, back, dist = geod.inv(la_lon, la_lat, berlin_lon, berlin_lat) #根据来自的地方反算结束经纬度和方位 x,y,bearing = geod.fwd(berlin_lon,berlin_lat,back,dist) #传递始末点和点数，得到等距坐标 coords = geod.npts(la_lon,la_lat,berlin_lon,berlin_lat,100) for i in range(3): print(coords[i])
