1 前言
地理坐标系，是使用三维球面来定义地球表面位置，以实现通过经纬度对地球表面点位引用的坐标系。 地理坐标系经过地图投影操作后就变成了投影坐标系。而地图投影是按照一定的数学法则将地球椭球面上点的经维度坐标转换到平面上的直角坐标。
2 流程
2.1 矢量数据地理坐标转投影坐标
要素feature的形状geometry是由一系列点坐标构成的，将每个要素的形状点一一进行坐标转换即可。
下面以WGS84坐标转UTM投影为例：
from osgeo import ogr,osr,gdal
import glob
import osdef 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.RasterYSizelr_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 = NoneUTM_raster_ds= None来源：应用推广部
供稿：技术研发部
编辑：方梅
