Python实现矢量与栅格数据投影转换及代码示例

宋时安

发布时间：2026-10-13| 3人浏览过

原创

本文详细介绍在Python中如何使用GDAL和OGR库进行矢量数据与栅格数据的投影转换。通过osr和proj4库实现坐标系统一，并演示从WGS84到Web墨卡托的转换过程，包含完整的代码示例和参数设置说明。

Python矢量与栅格数据投影转换概述

本文主要讲解利用Python语言将矢量数据转换为栅格数据的具体实现方法，文中提供了详细的代码示例，对读者的学习或实际工作具有参考价值。

在数据处理过程中，投影包osr与proj4的使用至关重要。

osr投影转换示例

以下是关于osr投影转换的具体示例。首先导入osr和ogr模块，定义源坐标系（如WGS84，EPSG:4326）和目标坐标系（如Google Maps使用的Web墨卡托，EPSG:3857）。通过osr.CoordinateTransformation创建转换对象，可以处理单个点的坐标转换，例如将经纬度(117, 40)转换为目标坐标系下的坐标。此外，还支持点数组的批量转换，如[(117,40),(117.5,39.5)]。对于OGR几何对象，可以通过CreateGeometryFromWkt创建几何体，利用Transform方法进行投影变换，并验证变换前后的WKT格式及坐标值。


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]))

栅格数据投影转换

接下来介绍栅格数据的投影转换方法。利用GDAL库，首先定义源图像的投影（如EPSG:32650）和目标投影（如EPSG:3857），并创建坐标转换对象。打开源TIFF文件后，获取其仿射变换矩阵（GeoTransform）。通过坐标转换对象计算源图像左上角和右下角在目标坐标系下的新坐标，从而确定新图像的边界。使用GTiff驱动创建新的空白栅格文件，保持行列数和波段数与原图一致，但数值类型可指定（如GDT_UInt16）。计算新图像的分辨率，并设置新的仿射变换矩阵和投影信息。最后，调用gdal.ReprojectImage函数，指定源和目标投影字符串及重采样方法（如双线性插值gdal.GRA_Bilinear），完成投影转换并关闭文件对象。


#栅格数据投影转换
from osgeo import gdal,osr
from 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 = None
ts= None

以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持本站。