矢量数据投影转换
作者:阿振
邮箱:tanzhenyugis@163.com
博客:https://blog.csdn.net/theonegis/article/details/80089375
修改时间:2018-06-03
声明:本文为博主原创文章,转载请注明原文出处
案例说明
接着上一篇博文中,我们得到了WGS84坐标系下的中国省区图,而我们一般中国地图中使用的是割圆锥投影。
由于我国位于中纬度地区,中国地图和分省地图经常采用割圆锥投影,中国地图的中央经线常位于东经105度,两条标准纬线分别为北纬25度和北纬47度,而各省的参数可根据地理位置和轮廓形状初步加以判定。
在SpatialReference中查到我们一般使用的中国地图投影为:http://spatialreference.org/ref/sr-org/8657
PROJ4格式的定义为:+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库,我们有两种方法进行矢量数据的重投影:
使用命令工具及其对应的命令行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格式的投影定义都OK
GDAL对该命令的封装的C/C++函数是GDALVectorTranslate()
,Python中是gdal.VectorTranslate()
使用GDAL提供的基本API进行实现
如果要自己利用基本API函数实现的话,基本思路如下:
- 利用
osgeo.ogr.Driver.CreateDataSource()
创建输出数据
- 根据源文件创建目标文件的属性字段定义
- 利用
osgeo.osr.CoordinateTransformation
对象将源文件中的Geometry对象转为目标文件中的Geometry对象(其实质是进行不同投影系统下空间几何体的坐标转换)
- 遍历源文件,依次将所有几何体的Geometry及其属性写入目标文件
代码实现
- 调用
gdal.VectorTranslate()
命令行工具的包装函数实现:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
| from osgeo import gdal import os os.environ['SHAPE_ENCODING'] = "utf-8"
src_file = 'China.shp' dst_file = 'China_Reprojected.shp'
srs_def = """+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函数实现
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
| from osgeo import ogr from osgeo import osr import os os.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()
srs_def = """+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)
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) 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_ds del dst_ds
dst_srs.MorphToESRI() (dst_path, dst_name) = os.path.split(dst_file) with open(dst_path + os.pathsep + dst_name + '.prj', 'w') as f: f.write(dst_srs.ExportToWkt())
|