如何使用Python进行地图投影转换操作？

本文共计836个文字，预计阅读时间需要4分钟。

一、使用OSR和PROJ4

1. OSR转换示例pythonfrom osgeo import osr

创建源坐标系统和目标坐标系统src_srs=osr.SpatialReference()src_srs.ImportFromEPSG(4326) # 假设源坐标系统为WGS84

dst_srs=osr.SpatialReference()dst_srs.ImportFromEPSG(3857) # 目标坐标系统为Web Mercator

创建转换器transformer=osr.CoordinateTransformation(src_srs, dst_srs)

转换坐标src_point=(114.0647, 22.52847)dst_point=transformer.TransformPoint(src_point[0], src_point[1])

print(转换后的坐标：, dst_point)

2.PROJ4使用示例

pythonfrom osgeo import ogr

创建数据源driver=ogr.GetDriverByName('ESRI Shapefile')data_source=driver.CreateDataSource('output.shp')

创建图层layer=data_source.CreateLayer('points', srs=src_srs)

创建字段field_name=ogr.FieldDefn('name', ogr.OFTString)field_name.SetWidth(50)layer.CreateField(field_name)

创建特征feature=ogr.Feature(layer.GetLayerDefn())feature.SetField('name', 'Point 1')geometry=ogr.Geometry(ogr.wkbPoint)geometry.AddPoint(*src_point)feature.SetGeometry(geometry)layer.CreateFeature(feature)feature=None

保存数据源data_source=None

一、投影包osr与proj4的使用

1.osr投影转换示例

from osgeo import osr,ogr
#定义投影
#wgs84
source=osr.SpatialReference()
source.ImportFromEPSG(4326)
#google
target=osr.SpatialReference()
target.ImportFromEPSG(3857)
#简单投影转换
coordTrans=osr.CoordinateTransformation(source,target)
#点转换
coordTrans.TransformPoint(117,40)
#点数组转换
coordTrans.TransformPoints([(117,40),(117.5,39.5)])
#SF几何对象转换
g=ogr.CreateGeometryFromWkt("POINT(117 40)")
#转换前wgs84
print(g.ExportToWkt())
print(g.GetX(),g.GetY())
#转换后google
g.Transform(coordTrans)
print(g.ExportToWkt())
print(g.GetX(),g.GetY())
2.投影转换示例

from pyproj import Proj,Geod,transform
# projection 1: UTM zone 15, grs80 ellipse, NAD83 datum
# (defined by epsg code 26915)
p1 = Proj(init=‘epsg:26915‘)
# projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum
p2 = Proj(init=‘epsg:26715‘)
#点的转换(首先将地理坐标转换成p1投影坐标系下的平面直角坐标，再将x1,y1转换到p2投影坐标系下，最后将p2投影坐标系下的平面直角坐标转换成地理坐标)
x1,y1=p1(-92.199881,38.56694)
x2, y2 = transform(p1,p2,x1,y1)
print(‘%9.3f %11.3f‘ % (x1,y1))
print(‘%9.3f %11.3f‘ % (x2,y2))
print(‘%8.3f %5.3f‘ % p2(x2,y2,inverse=True))
#点数组的转换
lats = (38.83,39.32,38.75)
lons = (-92.22,-94.72,-90.37)
x1,y1=p1(lons,lats)
x2,y2=transform(p1,p2,x1,y1)
xy=x1+y1
print(‘%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f‘ % xy)
xy=x2+y2
print(‘%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f‘ % xy)
lons, lats = p2(x2,y2,inverse=True)
xy=lons+lats
print(‘%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f‘ % xy)
p1 = Proj(proj=‘latlong‘,datum=‘WGS84‘)
x1 = -111.5; y1 = 45.25919444444
p2 = Proj(proj="utm",zone=10,datum=‘NAD27‘)
x2, y2 = transform(p1, p2, x1, y1)
print("%s %s" % (str(x2)[:9],str(y2)[:9]))

二、栅格数据投影转换

#栅格数据投影转换from osgeo import gdal,osrfrom osgeo.gdalconst import *#源图像投影source=osr.SpatialReference()source.ImportFromEPSG(32650)#目标图像投影target=osr.SpatialReference()target.ImportFromEPSG(3857)coordTrans=osr.CoordinateTransformation(source,target)#打开源图像文件ds=gdal.Open("fdem.tif")#仿射矩阵六参数mat=ds.GetGeoTransform()#源图像的左上角与右下角像素，在目标图像中的坐标(ulx, uly, ulz)=coordTrans.TransformPoint(mat[0],mat[3])(lrx, lry, lrz ) = coordTrans.TransformPoint(mat[0] + mat[1]*ds.RasterXSize, mat[3] + mat[5]* ds.RasterYSize )#创建目标图像文件（空白图像），行列数、波段数以及数值类型仍等同原图像driver=gdal.GetDriverByName("GTiff")ts=driver.Create("fdem_lonlat.tif",ds.RasterXSize,ds.RasterYSize,1,GDT_UInt16)#转换后图像的分辨率resolution=(int)((lrx-ulx)/ds.RasterXSize)#转换后图像的六个放射变换参数mat2=[ulx, resolution,0,uly,0, -resolution]ts.SetGeoTransform(mat2)ts.SetProjection(target.ExportToWkt())#投影转换后需要做重采样gdal.ReprojectImage(ds, ts, source.ExportToWkt(), target.ExportToWkt(), gdal.GRA_Bilinear)#关闭ds = Nonets= None