Python开发地理信息系统 python与地理信息系统

转载

mob64ca141a683a 2023-09-16 11:17:53

文章标签 Python开发地理信息系统 pycharm python qgis 坐标系统 文章分类 Python 后端开发

第5章 Python与地理信息系统

5.4重投影

借助osr实现数据的重投影一个Shapefile文件，不过现在这种使用方式很少了，可借助其他第三方库进行重投影。

在本示例中，将使用包含Lambert等角投影的纽约市博物馆和画廊位置的点Shapefile文件。可以把它转换成WGS84坐标系统。你可以通过如下地址获取该示例Shapefile文件的zip压缩格式。

from osgeo import ogr
from osgeo import osr
import os
import shutil

#1、定义文件路径，包括原始文件和重投影后的文件
#尽量命名绝对路径
srcName = r"D:\xxx\Jgy_PIE\3Python_Code\Python_proj\01Python_Learn\Learning_GeospatialAnalysiswithPython\NO_5_Python与地理信息系统\5-4重投影\NYC_MUSEUMS_LAMBERT\NYC_MUSEUMS_LAMBERT.shp"
tgtName = "GEO-WGS84.shp"

# 2、定义坐标系统
#osr.SpatialReference 和 osr.CoordinateTransformation 类
# 提供了用来描绘坐标系统（投影和基准面）以及坐标系统相互转换的服务。
tgt_spatRef = osr.SpatialReference()
# 定义要转换后的坐标 4326 是 WGS84
tgt_spatRef.ImportFromEPSG(4326)

#我国常用坐标系的对应编码
# spatialReference.importFromEPSG(4326)WGS84
# spatialReference.importFromEPSG(4214)BeiJing54
# spatialReference.importFromEPSG(4610)XIAN80
# spatialReference.importFromEPSG(4490)CGCS2000

#3、获取矢量文件
# 要读取某种类型的数据，必须要先载入数据驱动，也就是初始化一个对象，
# 让它“知道”某种数据结构。
driver = ogr.GetDriverByName("ESRI Shapefile")

# 其中update为0是只读，为1是可写,注意open(<filename>, <update>)中
# filename一定要写绝对路径！
src = driver.Open(srcName, 0)
srcLyr = src.GetLayer()

src_spatRef = srcLyr.GetSpatialRef()
# 如果文件存在，删除
if os.path.exists(tgtName):
    driver.DeleteDataSource(tgtName)
    
# 4、创建新文件，但是这个文件不能已经存在了，否则会出错
tgt = driver.CreateDataSource(tgtName)
lyrName = os.path.splitext(tgtName)[0]

#使用WKB格式声明几何图形
# wkt(OGC well-known text)和wkb(OGC well-known binary)是OGC制定的空间数据的组织规范，
# wkt是以文本形式描述，wkb是以二进制形式描述
# 使用wkt和wkb能够很好到和其他系统进行数据交换，目前大部分支持空间数据存储的数据库构造空间数据都采用这两种方式。

# 创建矢量，（图层名称，几何类型wkbPoint为点类型）
tgtLyr = tgt.CreateLayer(lyrName, geom_type=ogr.wkbPoint)
featDef = srcLyr.GetLayerDefn()
# 投影转换，将New York 3104 转换为WGS84
trans = osr.CoordinateTransformation(src_spatRef, tgt_spatRef)
srcFeat = srcLyr.GetNextFeature()

#  5、遍历每个点要素进行投影转换，遍历模板：
# feat = layer.GetNextFeature()  #读取下一个
# while feat:
#     feat = layer.GetNextFeature()

while srcFeat:
    geom = srcFeat.GetGeometryRef()
    geom.Transform(trans)
    feature = ogr.Feature(featDef)
    feature.SetGeometry(geom)
    tgtLyr.CreateFeature(feature)
    # 关闭数据源，相当于文件系统操作中的关闭文件,释放内存
    feature.Destroy()
    srcFeat.Destroy()
    srcFeat = srcLyr.GetNextFeature()
src.Destroy()
tgt.Destroy()

#6、导出投影文件将几何图形转换为Esri的WKT格式

tgt_spatRef.MorphToESRI()
prj = open(lyrName + ".prj", "w")
prj.write(tgt_spatRef.ExportToWkt())
prj.close()
srcDbf = os.path.splitext(srcName)[0] + ".dbf"
tgtDbf = lyrName + ".dbf"
shutil.copyfile(srcDbf, tgtDbf)

投影结果对比

原始影像：

Python开发地理信息系统 python与地理信息系统_qgis