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來源：應用推廣部
供稿：技術研發部
編輯：方梅
——9)
)
