Python GDAL | 对shapefile 数据进行重投影(坐标系转换)

31 篇文章 8 订阅
from osgeo import ogr, osr
from osgeo import gdal
import os
def reproject(inputfile,outputfile,layername,insrs,outsrs):
    gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")    
    gdal.SetConfigOption("SHAPE_ENCODING","") 
    driver = ogr.GetDriverByName('ESRI Shapefile')

# create the CoordinateTransformation
    coordTrans = osr.CoordinateTransformation(insrs, outsrs)

# get the input layer
    inDataSet = driver.Open(inputfile)
    inLayer = inDataSet.GetLayer()

# create the output layer
    outputShapefile = outputfile
    outDataSet = driver.CreateDataSource(outputShapefile)
    print inLayer.GetGeomType()
    outLayer = outDataSet.CreateLayer(layername,geom_type=inLayer.GetGeomType() )

    # add fields
    inLayerDefn = inLayer.GetLayerDefn()
    for i in range(0, inLayerDefn.GetFieldCount()):
        fieldDefn = inLayerDefn.GetFieldDefn(i)
        outLayer.CreateField(fieldDefn)

# get the output layer's feature definition
    outLayerDefn = outLayer.GetLayerDefn()

# loop through the input features
    inFeature = inLayer.GetNextFeature()
    while inFeature:
    # get the input geometry
        geom = inFeature.GetGeometryRef()
    # reproject the geometry
        geom.Transform(coordTrans)
    # create a new feature
        outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
        outFeature.SetGeometry(geom)
        for i in range(0, outLayerDefn.GetFieldCount()):
            outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
        outLayer.CreateFeature(outFeature)
    # destroy the features and get the next input feature
        outFeature.Destroy()
        inFeature.Destroy()
        inFeature = inLayer.GetNextFeature()

# close the shapefiles
    inDataSet.Destroy()
    outDataSet.Destroy()

 

  • 2
    点赞
  • 14
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
你可以使用GDAL库中的Python绑定来进行shp文件的投影转换。首先,确保你已经安装了GDAL库。然后,按照以下步骤进行操作: 1. 导入必要的库: ```python from osgeo import ogr, osr ``` 2. 打开原始的shp文件: ```python source = ogr.Open('path/to/source.shp') layer = source.GetLayer() ``` 3. 创建一个新的shp文件作为输出: ```python driver = ogr.GetDriverByName('ESRI Shapefile') output = driver.CreateDataSource('path/to/output.shp') outLayer = output.CreateLayer('output', geom_type=ogr.wkbPolygon) ``` 4. 定义原始坐标系和目标坐标系: ```python sourceSR = layer.GetSpatialRef() targetSR = osr.SpatialReference() targetSR.ImportFromEPSG(targetEPSG) # 替换targetEPSG为你想要的目标EPSG代码 ``` 5. 创建一个坐标转换器: ```python transform = osr.CoordinateTransformation(sourceSR, targetSR) ``` 6. 遍历原始图层中的要素,并进行投影转换: ```python for feature in layer: geom = feature.GetGeometryRef() geom.Transform(transform) # 创建新要素并将转换后的几何体添加到新图层中 newFeature = ogr.Feature(outLayer.GetLayerDefn()) newFeature.SetGeometry(geom) outLayer.CreateFeature(newFeature) newFeature = None source = None output = None ``` 这样,你就可以将原始shp文件中的几何体投影到目标坐标系,并保存为新的shp文件。 注意:在代码中,将`path/to/source.shp`和`path/to/output.shp`替换为你实际的文件路径,将`targetEPSG`替换为你想要的目标EPSG代码。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值