python+gdal地理坐标转投影坐标

1 前言

地理坐标系，是使用三维球面来定义地球表面位置，以实现通过经纬度对地球表面点位引用的坐标系。?地理坐标系经过地图投影操作后就变成了投影坐标系。而地图投影是按照一定的数学法则将地球椭球面上点的经维度坐标转换到平面上的直角坐标。

2 流程

2.1 矢量数据地理坐标转投影坐标

要素feature的形状geometry是由一系列点坐标构成的，将每个要素的形状点一一进行坐标转换即可。

下面以WGS84坐标转UTM投影为例：

from osgeo import ogr,osr,gdal
import glob
import os

def vecter_WGS2UTM(shp_path, UTM_shp_path, longitude, is_north):
    
    ds = ogr.Open(shp_path)
    layer = ds.GetLayer(0)
    
    driver = ogr.GetDriverByName('ESRI Shapefile')
    # 创建输出文件
    if os.path.exists(UTM_shp_path):
        driver.DeleteDataSource(UTM_shp_path)
    
    out_ds = driver.CreateDataSource(UTM_shp_path)
    outlayer = out_ds.CreateLayer(UTM_shp_path[:-4],geom_type = ogr.wkbPolygon)
    
    
    # 当前地理参考
    spatialRef = layer.GetFeature(0).GetGeometryRef().GetSpatialReference()
    
    # 根据经度计算UTM区号,进而定义UTM投影
    zone = str(int(longitude/6) + 31)
    zone = int('326' + zone) if is_north else int('327' + zone)
    UTM_spatialRef = osr.SpatialReference()
    UTM_spatialRef.ImportFromEPSG(zone)
    
    # 投影转换
    coordinate_transfor = osr.CoordinateTransformation(spatialRef, UTM_spatialRef)
    
    
    # 定义输出属性表信息
    feature = layer.GetFeature(0)
    field_count = feature.GetFieldCount()
    field_names = []
    for attr in range(field_count):
        field_defn = feature.GetFieldDefnRef(attr)
        field_names.append(field_defn.GetName())
        outlayer.CreateField(field_defn)

    
    out_fielddefn = outlayer.GetLayerDefn()
    
    for feature in layer:
        
        geometry = feature.GetGeometryRef()
        geometry.Transform(coordinate_transfor)
        
        out_feature = ogr.Feature(out_fielddefn)
        out_feature.SetGeometry(geometry)
        for field_name in field_names:
            out_feature.SetField(field_name,feature.GetField(field_name))
        outlayer.CreateFeature(out_feature)
    
        feature.Destroy()
        out_feature.Destroy()

    # 清除缓存
    ds.Destroy()
    out_ds.Destroy()
    # 保存投影文件
    UTM_spatialRef.MorphFromESRI()
    prj_path = UTM_shp_path.replace(".shp",".prj")
    fn = open(prj_path,'w')
    fn.write(UTM_spatialRef.ExportToWkt())
    fn.close()

2.2 栅格数据地理坐标转投影坐标

栅格数据每个像素的地理/投影坐标是由仿射矩阵六参数和像素坐标计算得来的，所以先将仿射矩阵六参数进行转换，之后对栅格数据重采样即可。

下面以WGS84坐标转UTM投影为例：

def raster_WGS2UTM(raster_path, UTM_raster_path, longitude, is_north):
    
    raster_ds = gdal.Open(raster_path)
    raster_type = raster_ds.GetRasterBand(1).DataType
    
    # 栅格投影
    spatialRef = osr.SpatialReference()
    spatialRef.ImportFromWkt(raster_ds.GetProjection())
    
    # 根据经度计算UTM区号,进而定义UTM投影
    zone = str(int(longitude/6) + 31)
    zone = int('326' + zone) if is_north else int('327' + zone)
    UTM_spatialRef = osr.SpatialReference()
    UTM_spatialRef.ImportFromEPSG(zone)
    
    # 投影转换
    coordinate_transfor = osr.CoordinateTransformation(spatialRef, UTM_spatialRef)
    
    # 仿射矩阵六参数
    geotransform = raster_ds.GetGeoTransform()
    # 左上角upper left、右下角lower right坐标
    ul_x = geotransform[0]
    ul_y = geotransform[3]
    lr_x = geotransform[0]+geotransform[1]*raster_ds.RasterXSize+geotransform[2]*raster_ds.RasterYSize
    lr_y = geotransform[3]+geotransform[4]*raster_ds.RasterYSize+geotransform[5]*raster_ds.RasterYSize
    
    # 左上角、右下角在目标投影中的坐标
    (UTM_ul_x, UTM_ul_y, UTM_ul_z) = coordinate_transfor.TransformPoint(ul_y, ul_x)
    (UTM_lr_x, UTM_lr_y, UTM_lr_z) = coordinate_transfor.TransformPoint(lr_y, lr_x)
    
    
    # 创建目标图像文件
    driver = gdal.GetDriverByName("GTiff")
    UTM_raster_ds = driver.Create(UTM_raster_path, 
                                  raster_ds.RasterXSize,
                                  raster_ds.RasterYSize,
                                  raster_ds.RasterCount,
                                  raster_type)
    # 转换后图像的分辨率
    resolution = (UTM_lr_x-UTM_ul_x)/raster_ds.RasterXSize
    # 转换后图像的六个放射变换参数
    UTM_transform = [UTM_ul_x, resolution, 0, UTM_ul_y, 0, -resolution]
    
    UTM_raster_ds.SetGeoTransform(UTM_transform)
    UTM_raster_ds.SetProjection(UTM_spatialRef.ExportToWkt())
    # 投影转换后需要做重采样
    gdal.ReprojectImage(raster_ds, UTM_raster_ds, spatialRef.ExportToWkt(), 
                        UTM_spatialRef.ExportToWkt(), gdal.GRA_Bilinear)
    # 关闭
    raster_ds = None
    UTM_raster_ds= None

来源：应用推广部

供稿：技术研发部

编辑：方梅