GDAL投影转换/坐标偏移笔记

1.OSR

#               OSR（矢量数据投影）
#作用：投影坐标系之间转换、地理坐标和投影坐标之间转换
#可用于几何对象和点（点属于几何对象）
from osgeo import gdal 
import osr
peters_sr = osr.SpatialReference()
peters_sr.ImportFromProj4('...')
ct = osr.CoordinateTransformation(web_mercator_sr, peters_sr)
#world对象为web_mercator_sr投影，但没有分配srs
world.Transform(ct)
#对几何对象的转换
ct.TransformPoint(x, y)
#对点的转换

#如果几何对象分配有srs，转换方法如下
world.TransformTo(web_mercator_sr)

#读取投影的地理基准，用于和地理坐标进行转换
geosrs = peters_sr.CloneGeogCS() 

2.PROJ

#                PROJ（矢量数据投影）
#作用：投影坐标之间转换、地理坐标和投影坐标之间转换
import pyproj
utm_proj = pyproj.Proj('+proj=utm +zone=31 +ellps=WGS84')
x, y = utm_proj(2.294, 48.858)  #地理坐标转化为投影坐标
x1, y1 = utm_proj(x, y, inverse=True)  #逆变换
wgs84 = pyproj.Proj('')
nad27 = pyproj.Proj('')
x, y = pyproj.transform(wgs84, nad27, 580744, 4504695) #投影坐标之间转换

3.栅格数据重投影

#               栅格数据重投影        
# 栅格数据也可以重投影，但比矢量数据投影更复杂。
# 栅格数据需要处理栅格数据中像元会弯曲和移动的事实，一对一的映射并不存在
# 通常用最近邻域插值法、双线性插值和三次卷积插值法进行插值
# 方法：1.gdalwarp 2.AutoCreateWarpedVRT 
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')  #UTM转无投影，即地理坐标系
old_ds = gdal.Open('nat_color.tif')
vrt_ds = gdal.AutoCreateWarpedVRT(old_ds, None, srs.ExportToWkt(), gdal.GRA_Bilinear) #第一个默认值None，使用源栅格数据的srs；第二个如果None，表示不发生重投影
gdal.GetDriverByName('gtiff').CreateCopy('nat_color_wgs84.tif', vrt_ds)      #该函数返回一个数据集对象，用CreateCopy函数保存为gtiff    

4.GetGeoTransform

#         GetGeoTransform，在真实坐标和栅格数据坐标具有相同srs情况下，计算坐标偏移。
#作用：图像坐标（行列号）和现实世界坐标（投影坐标或地理坐标）之间的转换。是仿射变换，不是投影转换，和上面不同。
#0、3 x\y坐标 起始点现实世界坐标  1、5 像素宽度和高度  2、4 x\y方向旋转角
gt = ds.GetGeoTransform() #正变换：图像坐标到现实世界坐标。正变换时输入行列号，输出的现实世界坐标是栅格图像srs下的坐标
inv_gt = gdal.InvGeoTransform(gt) #逆变换：现实世界坐标到图像坐标
offsets = gdal.ApplyGeoTransform(inv_gt, 465200, 5296000) #逆变换：输入的投影坐标具有和栅格图像的相同的srs
xoff, yoff = map(int, offsets)  #取整

5.gdal.Transformer

#        gdal.Transformer，可计算相同srs下的坐标偏移；不能用于不同srs投影转换
#作用：现实世界坐标（投影坐标或地理坐标）与图像坐标（行列号）之间的转换、两个栅格图像之间像素坐标偏移（行列号），如镶嵌
#原理就是在相同srs情况下，计算图1的像素坐标到现实世界坐标的偏移，再从现实世界坐标偏移到图2的像素坐标。其实就是两次仿射变换（正、逆），从而把图1的像素坐标偏移到图2的像素坐标。
#所以不能用于不同srs情况，因为该函数没有内置不同srs的投影转换公式。只能用于相同srs下，两个栅格数据集坐标的偏移。
#这里in_ds和out_ds具有相同srs。转换目的是为了把不同栅格数据的图像坐标（行列号）进行偏移，方便镶嵌
trans = gdal.Transformer(in_ds, out_ds, [])  #空白用于设置转换器选项
success, xyz = trans.TransformPoint(False, 0, 0) #False基于源数据计算目标栅格，反之为True。起始坐标为左上角 0，0
x,y,z = map(int, xyz)

#图像坐标和现实世界坐标之间转换
trans = gdal.Transformer(out_ds, None, [])
success, xyz = trans.TransformPoint(0, 1078, 648)

 

个人想法：

1.OSR \ PROJ \ AutoCreateWarpedVRT 分别用于矢量/点、栅格图像的投影转换，栅格与前两者有所不同。

2.GetGeoTransform是仿射变换，计算坐标偏移。（现实世界坐标（地理坐标或投影坐标）与图像坐标行列号），其中Geo[0/3]是现实世界坐标而不是行列号。

3.gdal.Transform本书中用于相同srs，不同栅格图像坐标（行列号）偏移的计算。书中用于镶嵌。

4.proj和osr都进行投影转换（可用于地理坐标和投影坐标之间的转换，是一种特例）。

5.计算栅格之间像素（行列号）偏移，这和仿射变换参数有关，而不是投影转换。

6.投影转换用OSR PROJ（矢量）/AutoCreateWarpedVRT(栅格)；Geotransform用于相同srs，图像坐标和现实世界坐标转换；gdal.Transformer用于相同srs，栅格图像

间图像坐标偏移，也可以计算图像坐标和现实世界坐标转换。

7.osr.CoordinateTransformation是不同坐标系的投影转换；

而gdal.Transformer是相同srs下，不同栅格图像之间的图像坐标偏移（也可计算现实世界坐标和图像坐标之间的转换）。

两者可以组合使用，完成不同投影的栅格图像之间的图像坐标偏移。