GIS遥感数据工程【免费下载链接】gdalGDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.项目地址https://gitcode.com/gh_mirrors/gd/gdal点击查看免费下载空间参考系统Spatial Reference System, SRS是地理数据处理中绕不开的核心概念——无论是读取栅格投影信息、为矢量图层定义坐标系还是做坐标系转换都离不开它。本文以 GDAL 仓库中的 Python 空间参考 API 文档为主线系统讲解osgeo.osr模块下的SpatialReference、CoordinateTransformation、CoordinateTransformationOptions与CreateCoordinateTransformation四大核心对象并结合仓库源码SWIG 绑定与ogr/ogrct.cpp底层实现深入说明其工作原理。读完本文你将掌握如何在 GDAL Python 绑定中完成 CRS 的创建、查询、导出与坐标转换并理解轴序、坐标操作选择等底层机制。一、API 全景osr 模块中的空间参考对象官方 API 文档页面 spatial_ref_api.rst 通过 Sphinxautoclass/autofunction指令从 Python docstring 自动生成类、方法与函数的完整说明。该页面收录的内容包括对象说明对应底层 C 类osgeo.osr.SpatialReference空间参考系统CRS的 Python 代理OGRSpatialReferenceogr_spatialref.hosgeo.osr.CoordinateTransformation坐标转换对象OGRCoordinateTransformationogrct.cpposgeo.osr.CoordinateTransformationOptions坐标转换选项决定坐标操作的选取策略OGRCoordinateTransformationOptionsosgeo.osr.CreateCoordinateTransformation使用选项创建转换对象的模块级函数对应 C APIOCTNewCoordinateTransformationEx()这些对象的实际方法文档由 osr_spatialreference_docs.i、osr_coordinatetransformation_docs.i 和 osr_docs.i 中的%feature(docstring)提供并通过 osr.i 完成 SWIG 绑定。同一页面还通过exclude-members: thisown隐藏了 SWIG 内部的内存管理属性只暴露用户真正关心的地理语义方法。二、SpatialReference空间参考对象的创建与导入2.1 构造函数三种初始化方式SpatialReference的构造函数见 osr_spatialreference_docs.i 与 osr_python.i 中的__init__实现支持恰好一个参数四种写法from osgeo import osr # 方式一按 EPSG 代码创建内部转为 EPSG:xxxx 字符串交给 SetFromUserInput srs osr.SpatialReference(epsg5646) print(srs.GetName()) # NAD83 / Vermont (ftUS) # 方式二按用户输入字符串PROJ 字符串、WKT、ESRI 等均可 srs2 osr.SpatialReference(projutm zone18 datumWGS84) print(srs2.GetUTMZone()) # 18 # 方式三按 WKT 字符串 srs3 osr.SpatialReference(wktGEOGCS[WGS 84,DATUM[WGS_1984,...]]) # 方式四传入字典内部 json.dumps 后交给 SetFromUserInput srs4 osr.SpatialReference({proj: utm, zone: 18, datum: WGS84})从 osr_python.i 的源码可以看到构造时若同时传入多个位置参数会抛出ValueErrorepsg关键字会被拼成fEPSG:{epsg}字符串最终统一走SetFromUserInput()完成解析——这也是构造函数能接受如此多种输入格式的根本原因。2.2 导入方法从各种来源初始化 CRS方法输入返回值说明ImportFromEPSG(arg)intEPSG 代码OGRERR_NONE0或错误码从 PROJ 数据库中查询 EPSG 地理/投影/垂直 CRSGDAL 3.0 起与ImportFromEPSGA完全相同ImportFromEPSGA(arg)intEPSG 代码同上历史遗留接口GDAL 3.0 起等同ImportFromEPSGImportFromProj4(ppszInput)strPROJ 字符串OGRERR_NONE或OGRERR_CORRUPT_DATA解析projutm zone18 datumWGS84这类传统 PROJ.4 参数ImportFromWkt(ppszInput)strWKT 字符串OGRERR_NONE或OGRERR_CORRUPT_DATA从 WKT 1 字符串导入ImportFromUrl(url)strURLOGRERR_NONE或错误码下载 URL 指向的 SRS 定义并交给SetFromUserInputImportFromCF1(keyValues, units)dict错误码从 netCDF CF-1 conventions 的 grid mapping 导入Python 层做了 list 序列化预处理文档示例 srs osr.SpatialReference() srs.ImportFromEPSG(4326) 0 srs.ImportFromProj4(projutm zone18 datumWGS84) 0此外AutoIdentifyEPSG()可以为 CRS 中能够与 EPSG 标识符安全对应的部分自动补充权威代码成功返回OGRERR_NONE无法识别时返回OGRERR_UNSUPPORTED_SRS。2.3 查询方法读取 CRS 的关键属性名称与权威代码 srs osr.SpatialReference() srs.ImportFromEPSG(4326) 0 srs.GetAuthorityName(DATUM) # EPSG srs.GetAuthorityCode(DATUM) # 6326 srs.GetAuthorityCode(None) # 4326根节点GetAuthorityCode(target_key)/GetAuthorityName(target_key)接受节点路径如PROJCS、GEOGCS传None则取根元素。GetAttrValue(name, child0)按名称不区分大小写和 0 起始的子节点索引取属性值例如vt_sp.GetAttrValue(UNIT, 0)返回US survey foot。单位与椭球参数方法含义示例EPSG:5646GetAngularUnits()角度单位换算为弧度的系数EPSG:4326 下返回0.017453292519943295GetAngularUnitsName()角度单位名degreeGetLinearUnits()线性单位换算为米的系数0.30480060960121924US survey footGetLinearUnitsName()线性单位名US survey footGetTargetLinearUnits(target_key)指定节点如VERT_CS的线性单位—GetSemiMajor()/GetSemiMinor()椭球长短半轴GDAL 3.0 起以米为单位缺失时返回SRS_WGS84_SEMIMAJOR/SRS_WGS84_SEMIMINOR常量GetInvFlattening()椭球反扁率WGS84 为298.257223563NAD83 为298.257222101GetTOWGS84()/HasTOWGS84()获取/判断 TOWGS84 七参数—投影参数GetProjParm(name, default_val0.0)返回投影参数原始值GetNormProjParm(name, default_val0.0)返回归一化值线性参数归一到米、角度参数归一到度。参数名使用osr.SRS_PP_前缀常量 vt_sp osr.SpatialReference() vt_sp.ImportFromEPSG(5646) 0 vt_sp.GetProjParm(osr.SRS_PP_FALSE_EASTING) # 1640416.6667英尺原值 vt_sp.GetNormProjParm(osr.SRS_PP_FALSE_EASTING) # 500000.0000101601归一化为米轴信息与坐标纪元GDAL 3.x 动态 CRS 相关GetAxesCount()坐标轴数量。EPSG:4326返回 2而三维地理 CRSEPSG:4979返回 3。GetAxisName(target_key, iAxis)/GetAxisOrientation(target_key, iAxis)取轴名与朝向常量如osr.OAO_North、osr.OAO_East、osr.OAO_Up。GetAxisMappingStrategy()返回数据轴到 CRS 轴的映射策略——OAMS_TRADITIONAL_GIS_ORDER传统 GIS 顺序、OAMS_AUTHORITY_COMPLIANT与 CRS 轴一致、OAMS_CUSTOM由SetDataAxisToSRSAxisMapping自定义。GetDataAxisToSRSAxisMapping()返回实际映射元组。GetCoordinateEpoch()返回坐标纪元以十进制年份表示未设置时返回 0。GetAreaOfUse()返回 CRS 适用区域对象包含名称与四至经纬度 aou vt_sp.GetAreaOfUse() aou.name United States (USA) - Vermont - ... aou.west_lon_degree, aou.south_lat_degree, aou.east_lon_degree, aou.north_lat_degree (-73.44, 42.72, -71.5, 45.03)UTM 查询GetUTMZone()返回 UTM 带号南半球为负数、北半球为正数非 UTM 投影返回 0文档示例中osr.SpatialReference(projutm zone18 datumWGS84).GetUTMZone()得到 18。2.4 导出方法输出为不同格式方法输出格式说明ExportToWkt()WKT 1 字符串标准 WKT 1 导出ExportToPrettyWkt(simplifyFalse)美观排版的 WKT 1便于人工阅读可指定simplify简化ExportToProj4()PROJ.4 字符串官方明确警告不鼓励使用见OGRSpatialReference::exportToProj4ExportToPROJJSON(options)PROJJSON 字符串现代 JSON 格式options为 list/dict 控制输出细节ExportToCF1(options{})dict导出为 netCDF CF-1 定义Python 层会把字符串数值还原为 float 列表2.5 类型判断与比较类型判断IsProjected()、IsGeographic()、IsGeocentric()、IsVertical()、IsCompound()复合 CRS、IsDerivedGeographic()如旋转经纬度网格、IsDynamic()动态坐标 CRS、IsLocal()局部 CRS以及HasPointMotionOperation()是否关联点运动操作。等价比较IsSame(rhs, options)、IsSameGeogCS(rhs, options)、IsSameVertCS(rhs, options)返回 1 表示相同。注意比较时 TOWGS84 子句会被忽略见下文源码说明。复合 CRS 处理StripVertical()GDAL 3.6 新增将复合 CRS 转换为纯水平 CRS成功返回OGRERR_NONE。三、CoordinateTransformation坐标转换对象3.1 构造与 TransformPointCoordinateTransformation(src, dst)接受源与目标两个SpatialReference在 osr.i 的绑定中映射到底层OCTNewCoordinateTransformation()若传入第三个CoordinateTransformationOptions参数则走OCTNewCoordinateTransformationEx()。TransformPoint()是 Python 绑定中最常用的方法其逻辑在 osr_python.i 中以%feature(shadow)定义支持传入单个序列3 元素或 4 元素走_TransformPoint3Double/_TransformPoint4Double或按(x, y, z0.0)、(x, y, z, t)分别传参返回(x, y, z)或(x, y, z, t)元组 wgs84 osr.SpatialReference() wgs84.ImportFromEPSG(4326) 0 vt_sp osr.SpatialReference() vt_sp.ImportFromEPSG(5646) 0 ct osr.CoordinateTransformation(wgs84, vt_sp) # WGS84 经纬度 - Vermont State Plane 东/北坐标 ct.TransformPoint(44.26, -72.58) (1619458.11, 641509.19, 0.0) ct.TransformPoint(44.26, -72.58, 103) # 带高程 (1619458.11, 641509.19, 103.0)3.2 批量转换与边界转换TransformPoints(arg)批量转换arg可以是元组列表也可以是 2xN、3xN、4xN 的 NumPy 数组返回(x, y, z)或(x, y, z, t)元组列表 ct.TransformPoints([(44.26, -72.58), (44.26, -72.59)]) [(1619458.11, 641509.19, 0.0), (1616838.29, 641511.90, 0.0)] import numpy as np ct.TransformPoints(np.array([[44.26, -72.58], [44.26, -72.59]])) [(1619458.11, 641509.19, 0.0), (1616838.29, 641511.90, 0.0)]TransformBounds(minx, miny, maxx, maxy, densify_pts)转换包围盒沿边界加密采样以刻画非线性变换导致的边界弯曲densify_pts官方推荐取 21 ct.TransformBounds(44.2, -72.5, 44.3, -72.4, 21) (1640416.67, 619626.43, 1666641.49, 656096.76)TransformPointWithErrorCode(x, y, z, t)TransformPoint的变体额外返回错误码返回(x, y, z, t, error)元组对应 C APIOCTTransformEx()定义于 ogrct.cpp。3.3 CoordinateTransformationOptions精确控制坐标操作CoordinateTransformationOptions是 spatial_ref_api.rst 中专门收录的选项类对应 COGRCoordinateTransformationOptions提供以下方法SWIG 绑定见 osr.i方法作用SetAreaOfInterest(w, s, e, n)指定感兴趣区域四至经纬度用于在候选坐标操作中检索最合适者SetOperation(operation, inverseCTFalse)强制使用用户自定义的坐标操作PROJ 管道字符串一旦设置将无条件使用SetDesiredAccuracy(accuracy)设置期望精度米用于过滤候选操作SetBallparkAllowed(allowBallpark)是否允许使用大致精度的球面近似操作SetOnlyBest(onlyBest)是否只选择最佳操作使用方式options osr.CoordinateTransformationOptions() options.SetAreaOfInterest(-73.5, 42.7, -71.4, 45.1) # 限定在佛蒙特州范围 options.SetDesiredAccuracy(0.5) # 期望亚米级精度 ct osr.CoordinateTransformation(wgs84, vt_sp, options)模块级函数CreateCoordinateTransformation(src, dst, options)与直接构造等价同样对应OCTNewCoordinateTransformationEx见 osr_docs.i。四、源码视角坐标转换的底层机制4.1 创建流程OGRCreateCoordinateTransformationPython 层的构造最终汇聚到 ogr/ogrct.cpp 的OGRCreateCoordinateTransformation()。其核心流程是将源/目标 SRS 转为文本表示GetTextRepresentation先调用OGRProjCT::FindFromCache()尝试命中缓存未命中则new OGRProjCT()并通过Initialize()基于 PROJ 构建转换对象失败时返回nullptrPython 层表现为抛出异常或返回None取决于是否启用osr.UseExceptions()。从 ogrct.cpp 的注释可知几个关键行为源/目标 SRS 默认不应为NULL除非通过options提供了自定义坐标操作转换会遵循源与目标 SRS 声明的轴序及 data axis to SRS axis mapping如需 GDAL 3.0 之前的传统经纬度顺序行为可设置配置项OGR_CT_FORCE_TRADITIONAL_GIS_ORDERYES若options中定义了用户坐标操作管道则无条件使用若设置了感兴趣区域则用它检索最合适的转换且该转换会用于所有坐标变换即使点落在区域之外若未设置任何选项则检索候选操作列表并在每次Transform()调用时根据坐标集合的质心动态选择最优操作。4.2 操作选择策略OGR_CT_OP_SELECTIONGDAL 3.0.3 起OGR_CT_OP_SELECTION配置项可控制候选坐标操作的选取策略见 ogrct.cppPROJPROJ ≥ 6.3 时的默认值PROJproj_create_crs_to_crs()的默认行为对每个待转换点评估候选操作BEST_ACCURACY选择精度最佳的操作在单次Transform()调用坐标的平均点上决策当在 PROJ 8 下调用SetDesiredAccuracy()/SetBallparkAllowed()时会自动退化为该策略FIRST_MATCHING取候选列表中排在最前的操作精度未必最佳但适用区域通常更大GDAL 3.0.0–3.0.2 的默认行为。另外ogrct.cpp 还提到默认情况下若源/目标 SRS 通过代码引用了官方 CRS 且两者等价GDAL 会优先使用官方定义比较时忽略 TOWGS84 子句GDAL 3.4.1 起可用OGR_CT_PREFER_OFFICIAL_SRS_DEFNO关闭该行为强制使用用户给出的 SRS 定义。4.3 Transform 族 C API 映射批量转换最终落到 ogrct.cpp 的OCTTransform()单点与 ogrct.cpp 的OCTTransformBounds()边界转换边界加密策略即在此实现。TransformPointWithErrorCode则对应OCTTransformEx()允许在转换失败时获取具体错误码而非仅仅抛出异常。4.4 异常模式与测试验证与 GDAL 其它 Python 绑定一致osr模块支持UseExceptions()/DontUseExceptions()两种模式osr_python.i 中实现了_WarnIfUserHasNotSpecifiedIfUsingExceptions()在 GDAL 4.0 将默认启用异常前对未显式声明异常模式的用户发出FutureWarning。仓库的 autotest/osr 目录下存在大量与本文 API 一一对应的测试osr_basic.py基础 SRS 操作、osr_ct.py/osr_ct_proj.py坐标转换、osr_epsg.pyEPSG 导入、osr_proj4.py、osr_esri.py、osr_validate.py等可以作为学习各方法行为边界的参考样例。五、完整实战示例将上述 API 串联起来一个典型的定义 CRS → 查询属性 → 坐标转换工作流如下from osgeo import osr # 1. 定义源与目标 CRS wgs84 osr.SpatialReference() wgs84.ImportFromEPSG(4326) # WGS 84 经纬度 vt_sp osr.SpatialReference() vt_sp.ImportFromEPSG(5646) # NAD83 / Vermont (ftUS) # 2. 查询 CRS 关键信息 print(vt_sp.GetName()) # NAD83 / Vermont (ftUS) print(vt_sp.GetLinearUnitsName()) # US survey foot print(vt_sp.GetProjParm(osr.SRS_PP_FALSE_EASTING)) # 3. 建立转换并执行 ct osr.CoordinateTransformation(wgs84, vt_sp) easting, northing, z ct.TransformPoint(44.26, -72.58) print(fEasting: {easting:.2f}, Northing: {northing:.2f}) # 4. 可选用选项限制感兴趣区域并重新创建转换 opts osr.CoordinateTransformationOptions() opts.SetAreaOfInterest(-73.5, 42.7, -71.4, 45.1) ct2 osr.CreateCoordinateTransformation(wgs84, vt_sp, opts) # 5. 导出 WKT 用于持久化 wkt vt_sp.ExportToWkt()结语osgeo.osr的空间参考 API 是 GDAL Python 绑定中最高频使用的接口之一SpatialReference负责 CRS 的定义、解析、查询与导出CoordinateTransformation负责坐标系之间的转换而CoordinateTransformationOptions让高级用户可以精确控制坐标操作的选取策略。理解其背后的 SWIG 绑定结构osr.i、osr_python.i与 C 实现ogrct.cpp中的轴序处理、操作选择与缓存机制能帮助你在实际项目中规避经纬度顺序颠倒坐标操作选择不符合预期等常见陷阱写出正确可靠的地理坐标处理代码。赞分享GIS遥感数据工程【免费下载链接】gdalGDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.项目地址https://gitcode.com/gh_mirrors/gd/gdal点击查看免费下载相关推荐TVM Python Driver API 深度解析tvm.driver 命名空间与 tvm.compile 统一编译入口实战指南TVM Python Driver API 深度解析tvm.driver 命名空间与 tvm.compile 统一编译入口实战指南 本文围绕 TVMApac模型编译深度学习推理引擎autoskills 深度指南Cloudflare Sandbox SDK 完整 API 参考与实战解析autoskills 深度指南Cloudflare Sandbox SDK 完整 API 参考与实战解析 Cloudflare Sandbox SDK 在边缘Godot 动画系统深度解析AnimationNodeBlendSpace1D 一维混合空间实战指南Godot 动画系统深度解析AnimationNodeBlendSpace1D 一维混合空间实战指南 导读 AnimationNodeBlendSpace1D文档教程游戏开发上一篇10个使用openYuanrong最容易踩的坑安装、调用与调试FAQ大全下一篇一个端点服务整个团队MiniStack多账户与多区域隔离机制完全解析创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考