矢量数据投影转换

矢量数据投影转换

矢量数据投影转换

作者：阿振

案例说明

接着上一篇博文中，我们得到了WGS84坐标系下的中国省区图，而我们一般中国地图中使用的是割圆锥投影。

由于我国位于中纬度地区，中国地图和分省地图经常采用割圆锥投影，中国地图的中央经线常位于东经105度，两条标准纬线分别为北纬25度和北纬47度，而各省的参数可根据地理位置和轮廓形状初步加以判定。

在​​SpatialReference​​​中查到我们一般使用的中国地图投影为：​​+lat_1=25 +lat_2=47 +lat_0=30 +lon_0=105 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs​​

使用该投影，我们祖国雄鸡才会变得雄赳赳气昂昂，更好地展现我们神州大地的风采。

方法介绍

跟栅格数据投影转换一样，使用GDAL库，我们有两种方法进行矢量数据的重投影：

使用命令工具及其对应的命令行API接口进行转换（简单，准确，实践中一定要用这种方法）GDAL提供了​​ogr2ogr​​命令行工具进行矢量数据投影转换，命令如下：​​ogr2ogr -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=30 +lon_0=105 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs " China_Projected.shp China.shp​​​​-t_srs​​​选项制定输出数据投影，当然可以是ESPG，也可以是PROJ4或者OGC WKT格式的投影定义都OKGDAL对该命令的封装的C/C++函数是​​GDALVectorTranslate()​​,Python中是​​gdal.VectorTranslate()​​使用GDAL提供的基本API进行实现如果要自己利用基本API函数实现的话，基本思路如下：

利用​​osgeo.ogr.Driver.CreateDataSource()​​创建输出数据根据源文件创建目标文件的属性字段定义利用​​osgeo.osr.CoordinateTransformation​​对象将源文件中的Geometry对象转为目标文件中的Geometry对象（其实质是进行不同投影系统下空间几何体的坐标转换）遍历源文件，依次将所有几何体的Geometry及其属性写入目标文件

代码实现

调用​​gdal.VectorTranslate()​​命令行工具的包装函数实现：

from osgeo import gdalimport osos.environ['SHAPE_ENCODING'] = "utf-8"src_file = 'China.shp'dst_file = 'China_Reprojected.shp'# 使用命令行API转换# 输出数据投影定义，参考资料：= """+proj=aea +lat_1=25 +lat_2=47 +lat_0=30 +lon_0=105 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs """gdal.VectorTranslate(dst_file, src_file, dstSRS=srs_def, reproject=True)src_ds = ogr.Open(src_file)src_layer = src_ds.GetLayer(0)src_srs = src_layer.GetSpatialRef() # 输入数据投影

调用基本API函数实现

from osgeo import ogrfrom osgeo import osrimport osos.environ['SHAPE_ENCODING'] = "utf-8"src_file = 'China.shp'dst_file = 'China_Reprojected.shp'src_ds = ogr.Open(src_file)src_layer = src_ds.GetLayer(0)src_srs = src_layer.GetSpatialRef() # 输入数据投影# 输出数据投影定义，参考资料：= """+proj=aea +lat_1=25 +lat_2=47 +lat_0=30 +lon_0=105 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs """dst_srs = osr.SpatialReference()dst_srs.ImportFromProj4(srs_def)# 创建转换对象ctx = osr.CoordinateTransformation(src_srs, dst_srs)# 创建输出文件driver = ogr.GetDriverByName('ESRI Shapefile')dst_ds = driver.CreateDataSource(dst_file)dst_layer = dst_ds.CreateLayer('province', dst_srs, ogr.wkbPolygon)# 给输出文件图层添加属性定义layer_def = src_layer.GetLayerDefn()for i in range(layer_def.GetFieldCount()): field_def = layer_def.GetFieldDefn(i) dst_layer.CreateField(field_def)# 循环遍历源Shapefile中的几何体添加到目标文件中src_feature = src_layer.GetNextFeature()while src_feature: geometry = src_feature.GetGeometryRef() geometry.Transform(ctx) dst_feature = ogr.Feature(layer_def) dst_feature.SetGeometry(geometry) # 设置Geometry # 依次设置属性值 for i in range(layer_def.GetFieldCount()): field_def = layer_def.GetFieldDefn(i) field_name = field_def.GetName() dst_feature.SetField(field_name, src_feature.GetField(field_name)) dst_layer.CreateFeature(dst_feature) dst_feature = None src_feature = None src_feature = src_layer.GetNextFeature()dst_ds.FlushCache()del src_dsdel dst_ds# 创建Shapefile的prj投影文件dst_srs.MorphToESRI()(dst_path, dst_name) = os.path.split(dst_file)with open(dst_path + os.pathsep + dst_name + '.prj', 'w') as

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。