Python 使用gdal变换空间坐标系统

使用CoordinateTransformation函数转化不同的坐标系统

Gdal读取遥感图片

在使用GDAL库读取GeoTiff格式的遥感数据中，需要涉及到操作经纬度参数来提取图片中具体一块区域的图片。
读取到的一个dataset（对应GDALDataset类）是一个光栅数据以及和它有关系的信息的集合。特别地dataset包含了光栅数据的大小（像素、线等）。dataset同时也为对应的光栅数据指定了坐标系统。dataset本身还可以包含元数据，它们以一种键/值对 的方式来组织。GDAL的数据集是基于OpenGIS Grid Coverages的格式定义的。
dataset包含的地理空间坐标系统是不同的，有的是平面空间坐标，有的是大地空间坐标（如：经纬度），所以在不同的坐标空间之间转化在有的时候是非常必要的。

Dataset的坐标系统由OpenGIS WKT字符串定义，它包含了：

1. 一个全局的坐标系名称。 
2. 一个地理坐标系名称。 
3. 一个基准标识符。 
4. 椭球体的名字。长半轴(semi-major axis)和反扁率(inverse flattening)。 
5. 初子午线(prime meridian)名和其与格林威治子午线的偏移值。 
6. 投影方法类型（如横轴莫卡托）。 
7. 投影参数列表（如中央经线等）。 
8. 一个单位的名称和其到米和弧度单位的转换参数。 
9. 轴线的名称和顺序。
10. 在预定义的权威坐标系中的编码（如EPSG）。

在GDAL中，返回坐标系统的函数是GDALDataset::GetProjectionRef()。 它返回的坐标系统描述了地理参考坐标，暗含着仿射地理参考转换，这地理参考转换是由GDALDataset::GetGeoTransform()来返回。由GCPs地理参考坐标描述的坐标系统是由 GDALDataset::GetGCPProjection()返回的。
注意，返回的坐标系统字符串“”表示未知的地理参考坐标系统。
wkt对坐标转换很重要

转化坐标系的步骤（个人理解）：

  • 定义源坐标系，定义目标坐标系。使用 ob = osr.SpatialReference()方式建立
  • 分别对2个坐标系定义参数。由于源坐标系来自遥感地图，所以要从要遥感数据中获得wkt字符串。 wkt = dataset.GetProjectionRef()可以从数据集中获取wkt字符串，然后使用ob_ori.ImportFromWkt(wkt)函数，把参数导入到源坐标系。
  • 目标坐标系在这我们选择美国定义的’WGS84’大地坐标系，由于这个坐标系是已知常用的坐标系，所以我们使用ob_for.SetWellKnownGeogCS(“WGS84”)来定义系统参数。
  • 到现在，2个坐标系统都被定义好了。下一步是建立它们的关系—— ct = osr.CoordinateTransformation(object,origin），其中，origin是源坐标系，object是目标坐标系。
  • 最后，就是传入坐标和获取坐标。

函数定义

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
'''把得到的平面直角坐标转化成WGS84经纬度坐标'''
'''原始坐标的投影系统用sr来表示，是由读取的TIF遥感图像决定的'''
def wkt_WGS84(xsize,ysize,wkt):
    '''首先建立2个投影坐标系，然后分别导入相应的坐标参数'''
    object_Long = osr.SpatialReference()
    object_xy = osr.SpatialReference()
    object_xy.ImportFromWkt(wkt)
    object_Long.SetWellKnownGeogCS("WGS84")
    '''建立转化方程，返回WGS84经纬度'''
    ct = osr.CoordinateTransformation(object_xy,object_Long)
    Latitude,longitude,High = ct.TransformPoint(xsize,ysize)
    print ('经过转化的WGS84经纬度为(%s,%s)' %(Latitude,longitude))
    return Latitude,longitude
'''把经纬度转化成平面直角坐标系统'''
def WGS84_wkt(Latitude,longitude,wkt):
    object_Long = osr.SpatialReference()
    object_xy = osr.SpatialReference()
    object_xy.ImportFromWkt(wkt)
    object_Long.SetWellKnownGeogCS("WGS84")
    ct = osr.CoordinateTransformation(object_Long,object_xy)
    xsize,ysize,zsize = ct.TransformPoint(Latitude,longitude,0.0)
    xsize = float(int(xsize+0.5))
    ysize = float(int(ysize+0.5))
    print ('经过转化的仿射坐标为(%s,%s)' %(xsize,ysize))
    return xsize,ysize

仿射地理变换

要想在2个坐标系之间变换，必须得得到原本遥感数据中的坐标值，这个坐标值我称之为仿射坐标。
GDAL数据集有两种方式描述栅格位置（用点/线坐标系）以及地理参考坐标系之间的关系。 第一种也是比较常用的是使用仿射转换，另一种则是GCPs。
仿射变换由6个参数构成，它们由GetGeoTransform()返回它们把点/线坐标， 用下面的关系转将点/线影射到地理坐标：
Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
假设影像上面为北方，GT2和GT4参数为0，而GT1是象元宽，GT5是象元高， （GT0，GT3）点位置是影像的左上角。

将地理坐标转换为像素点的操作为：
FileX = (int)floor((GeoX * GT[5] - GeoY * GT[2] - GT[0] * GT[5] + GT[2] * GT[3]) / (GT[1] * GT[5] - GT[2] * GT[4]));
FileY = (int)floor((GeoX * GT[4] - GeoY * GT[1] - GT[0] * GT[4] + GT[1] * GT[3]) / (GT[2] * GT[4] - GT[1] * GT[5]));

1
2
3
4
5
6
7
8
9
10
11
12
13
def Geo_Points(FileX,FileY,GT):
    GeoX = GT[0] + FileX * GT[1] + FileY * GT[2]
    GeoY = GT[3] + FileX * GT[4] + FileY * GT[5]
    print ('转化后的仿射坐标为(%s,%s)' %(GeoX,GeoY))
    return GeoX,GeoY

def File_points(GeoX,GeoY,GT):
    FileX = float((GeoX * GT[5] - GeoY * GT[2] -  GT[0] * GT[5] + GT[2] * GT[3]) /  (GT[1] * GT[5] - GT[2] * GT[4]))
    FileY = float((GeoX * GT[4] - GeoY * GT[1] -  GT[0] * GT[4] + GT[1] * GT[3]) /  (GT[2] * GT[4] - GT[1] * GT[5]))
    FileX = int(FileX+0.5)
    FileY = int(FileY+0.5)
    print ('转化后的图像坐标为(%s,%s)' %(FileX,FileY))
    return FileX,FileY

测试代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
dataset = gdal.Open(filePath_1)
'''获得tif的坐标，其中左上角的坐标为[0],[3]'''
adfGeoTransform = dataset.GetGeoTransform()
'''获得tif的投影坐标系统'''
wkt = dataset.GetProjectionRef()

''' 右下角的仿射坐标x,y'''
nXSize = dataset.RasterXSize #列数
nYSize = dataset.RasterYSize #行数
zero_x,zero_y = Geo_Points(0,0,adfGeoTransform)
rd_x,rd_y = Geo_Points(nXSize,nYSize,adfGeoTransform)

print ('所读取遥感图片的仿射系统参数为[{}]'.format(adfGeoTransform))

'''测试，把原点和右下角点转换成经纬度)'''
print ('原点的仿射坐标为(%s,%s)' %(zero_x,zero_y))
WGS84_zero_x,WGS84_zero_y= wkt_WGS84(zero_x,zero_y,wkt)
print ('右下角的仿射坐标为(%s,%s)' %(rd_x,rd_y))
WGS84_rd_x,WGS84_rd_y = wkt_WGS84(rd_x,rd_y,wkt)

'''测试，把得到的2个经纬度转化回仿射坐标'''
print ('原点的经纬度为(%s,%s)' %(WGS84_zero_x,WGS84_zero_y))
xsize0,ysize0 = WGS84_wkt(WGS84_zero_x,WGS84_zero_y,wkt)

'''测试，把仿射坐标对应到图像坐标'''
print ('右下角的图像坐标为(%s,%s)' %(nXSize,nYSize))
print ('右下角的仿射坐标为(%s,%s)' %(rd_x,rd_y))
file_rd_x,file_rd_y = File_points(rd_x,rd_y,adfGeoTransform)

测试结果

1
2
3
4
5
6
7
8
9
10
>>> 所读取遥感图片的仿射系统参数为[(126885.0, 30.0, 0.0, 3625515.0, 0.0, -30.0)]
>>> 原点的仿射坐标为(126885.0,3625515.0)
>>> 经过转化的WGS84经纬度为(119.020437123,32.7043061635)
>>> 右下角的仿射坐标为(365715.0,3409185.0)
>>> 经过转化的WGS84经纬度为(121.596244068,30.8081674204)
>>> 原点的经纬度为(119.020437123,32.7043061635)
>>> 经过转化的仿射坐标为(126885.0,3625515.0)
>>> 右下角的图像坐标为(7961,7211)
>>> 右下角的仿射坐标为(365715.0,3409185.0)
>>> 转化后的图像坐标为(7961,7211)