用于在 CRS 格式之间进行读取、写入和转换的 GIS 包。
项目描述
PyCRS
PyCRS 是一个纯 Python GIS 包,用于在各种通用坐标参考系统 (CRS) 字符串和数据源格式之间进行读取、写入和转换。
目录
介绍
Python 应该有一个独立的 GIS 库,专注于坐标参考系统元数据。也就是说,一个专注于用于存储和表示 crs 定义的各种格式的库,包括 OGC WKT、ESRI WKT、Proj4 以及由 EPSG、ESRI 和 SR-ORG 等组织定义的各种短代码。在许多类型的 GIS 工作中,正确解析和转换这些格式是必不可少的。例如,当尝试使用 PyProj 从非 proj4 crs 格式转换坐标时。或者当想要将 crs 从 GeoJSON 文件转换为 .prj 文件时。或者只是将 crs 定义添加到以前缺少的文件中时。
当我创建 PyCRS 时,在 crs 格式之间读取和转换的唯一方法是使用广泛的 Python GDAL 套件及其 srs 子模块,但某些应用程序的要求可能会排除使用 GDAL。还有一些在线网站/服务,但这些仅允许部分查找或从一种格式单向转换为另一种格式。因此,我希望 PyCRS 将使轻量级应用程序更容易读取更广泛的数据文件,并正确解释和可能转换它们的 crs 定义。完全用 Python 编写,我也希望它能帮助澄清各种格式之间的差异,并让更多人更容易帮助它保持最新和无错误。
地位
目前,支持的格式为 OGC WKT (v1)、ESRI WKT、Proj4 以及可从 spatialreference.org 获得的任何 EPSG、ESRI 或 SR-ORG 代码。将来我希望添加对 OGC URN 标识符字符串和 GeoTIFF 文件标签的支持。
在某些情况下,PyCRS 不能完美地解析或在所有 crs 格式之间进行转换,基于(希望主要是外观)与 GDAL 等其他解析器结果的差异。在源代码库中有一个 tester.py 脚本,它使用http://www.remotesensing.org/geotiff/proj_list/上列出的一系列常用 crs 。目前,三种主要格式之间的加载和转换的总体成功率为 80-90%,并且使用每个 crs 渲染世界的视觉检查通常看起来是正确的。但是,从数学的角度来看,转换后的 crs 字符串是否在逻辑上彼此等价,这需要更详细的质量检查。
变化
1.0.2 (2020-03-09)
- 添加搜索功能以查找 epsg.io
- 根据 prj2epsg.org (#47, @fitoprincipe) 添加 to_epsg_code() 方法
- 查找 epsg 代码时切换到使用 epsg.io
- 通过切换到 SSL 修复损坏的空间参考链接
- 添加了 Clarke 1880 椭球
- 未指定椭球时使用默认基准椭球
- 找不到 proj4 椭球时更好的警告和错误
- 杂项错误修复
1.0.1 (2019-03-07)
- 其他文档字符串更改
- 更灵活的解析,以防投影不是第二个元素
- 更多信息异常
1.0.0 (2019-02-04)
- API 更改:
- 删除 CRS 包装类,而不是直接处理 GeogCS 或 ProjCS
- 将模块名称更改为 load.py 和 parse.py
- 将 Ellipsoid、Datum 和 Projection 移动到各自的模块
- 允许 proj4 输入和输出为 dict
- 修复:
- 添加了更多文档
- 修复 +f 解释
- 包括读取+rf参数
- 更好地支持 +a 和 +b 的椭球体
- 修复 Python 3 错误
- 修复 proj4 标准平行被忽略
- 解析之前被忽略的 WKT CS 名称
- 支持本初子午线城市名称
0.1.3 (2016-06-25)
- 修复了各种错误
- 适用于 Mac 和 Linux 的 Pip 安装修复程序
- Python 3 兼容性
0.1.2 (2015-08-05)
- 第一次正式发布
依赖项
纯 Python,无依赖。Python 2 和 3 兼容。
安装
PyCRS 是从命令行使用 pip 安装的:
pip install pycrs
它也可以将“pycrs”包文件夹放在可导入的位置,如“PythonXX/Lib/site-packages”。
例子
首先导入 pycrs 模块:
>>> import pycrs
创建 CS 实例
PyCRS 使用不同类型的 CS 类来表示和处理所有坐标参考系统。要创建一个,您可以从源加载它,从字符串解析它,从 CRS 代码中查找,或者从头开始构建它。让我们回顾一下创建 CS 类型实例的这些不同方法。
从外部源加载
如果您知道 crs 信息位于某个外部源中,PyCRS 提供了一些方便的函数来加载这些信息,所有这些都位于“pycrs.load”模块中。
从 Shapefile 加载
在大多数情况下,这意味着读取 shapefile 附带的 ESRI .prj 文件。PyCRS 有一个方便的功能来做到这一点:
>>> crs = pycrs.load.from_file("testfiles/natearth.prj")
从 GeoJSON 加载
同样的功能还支持从 GeoJSON 文件中读取 crs:
>>> crs = pycrs.load.from_file("testfiles/cshapes.geo.json")
从 URL 加载
如果您的 crs 未在文件中定义,而是在网页上定义为纯文本,则还有一个功能:
>>> crs = pycrs.load.from_url("http://spatialreference.org/ref/esri/54030/ogcwkt/")
从文本字符串解析
但是,在许多情况下,您的代码中可能已经有了字符串表示。这可能是如果您正在与其他库进行互操作,或者您已经从某个外部源读取它。在这些情况下,您可以使用“pycrs.parse”模块中可用的函数来创建 CS 类型实例。
从 proj4 字符串或 dict 解析
要从 proj4 字符串创建 CS 类型实例,您可以这样做:
>>> proj4 = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
>>> crs = pycrs.parse.from_proj4(proj4)
或者,如果您的 proj4 字符串表示为 dict:
>>> proj4_as_dict = dict(proj='robin', lon_0=0, x_0=0, y_0=0, ellps='WGS84', datum='WGS84', units='m')
>>> crs = pycrs.parse.from_proj4(proj4_as_dict)
从 ESRI WKT 字符串解析
ESRI WKT 格式是通常在 shapefile 的 .prj 文件中找到的格式。如果你已经从一个文件中加载了它,你可以像这样解析它:
>>> esri_wkt = 'PROJCS["World_Robinson",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Robinson"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],UNIT["Meter",1]]'
>>> crs = pycrs.parse.from_esri_wkt(esri_wkt)
从 OGC WKT 字符串解析
开放地理空间联盟 (OGC) WKT 格式是 ESRI WKT 的更新变体。只有微小的差异,但将来可能会得到更多支持。如果你已经把它作为一个字符串,你可以像这样解析它:
>>> ogc_wkt = 'PROJCS["World_Robinson",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Robinson"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],UNIT["Meter",1],AUTHORITY["EPSG","54030"]]'
>>> crs = pycrs.parse.from_ogc_wkt(ogc_wkt)
从未知字符串解析
最后,如果你不知道crs字符串的格式,你也可以让PyCRS为你自动检测和解析crs类型:
>>> for unknown in [proj4, esri_wkt, ogc_wkt]:
... crs = pycrs.parse.from_unknown_text(unknown)
查找坐标系代码
存储坐标系的另一种常用方法是通过可用于许多更常用的查找代码的查找代码。多个不同的机构已经定义了自己的代码集。
查找 EPSG 代码
要查找 EPSG 定义的代码:
>>> crs = pycrs.parse.from_epsg_code(4326)
查找 ESRI 代码
要查找 ESRI 定义的代码:
>>> crs = pycrs.parse.from_esri_code(54030)
查找 SR 代码
要查找由 spatialreference.org 定义的代码:
>>> crs = pycrs.parse.from_sr_code(42)
按名称或区域搜索坐标系
最后,如果您不知道特定 crs 的具体代码或定义,还有搜索功能可以在 epsg.io 网站上进行查找。例如,如果您想使用 Robinson 投影但不确定如何加载它,您可以交互式搜索和检查潜在匹配项:
>>> for match in pycrs.utils.search_name('robinson'):
... # do something
... pass
同样,如果您需要知道哪个投影适合特定国家/地区,您也可以搜索:
>>> for match in pycrs.utils.search_area('brazil'):
... # do something
... pass
您还可以使用各种说明符进行自定义搜索,例如特定区域的投影(说明符的完整列表可在https://github.com/maptiler/epsg.io获得):
>>> for match in pycrs.utils.search('name:utm area:belgium'):
... # do something
... pass
从搜索结果加载
搜索匹配结果包括 epsg 代码、proj4 和 wkt 表示,其中任何一个都可用于加载 CS 实例:
>>> topmatch = next(pycrs.utils.search_name('wgs 84'))
>>> crs = pycrs.parse.from_proj4(topmatch['proj4'])
检查 CS 实例
加载、解析、查找或创建坐标参考系统后,您最终会得到 pycrs.cs 模块中的一个 CS 类型实例。所有 CS 类型都是 pycrs.CS 的子类:
>>> isinstance(crs, pycrs.CS)
True
PyCRS 目前仅支持两种最常见的 CS 类型,地理或投影。
地理CS
地理参考系统将坐标保存在经纬度空间中,我们指定它的原因是因为定义地球形状的方式不同。举个例子,我们加载常用的WGS84地理坐标系:
>>> crs = pycrs.parse.from_epsg_code(4326)
如果它加载正确,您现在应该有一个 pycrs.GeogCS 实例,您还可以通过cs_type
属性检查它:
>>> isinstance(crs, pycrs.GeogCS)
True
>>> crs.cs_type
'Geographic'
通过 GeogCS 实例,我们可以进一步访问它的子组件和参数。例如,如果我们想检查数据,我们可以这样做:
>>> datum = crs.datum
>>> isinstance(datum, pycrs.elements.datums.WGS84)
True
或者椭球的逆展平因子:
>>> ellips = crs.datum.ellips
>>> ellips.inv_flat.value
298.257223563
有关如何检查 GeogCS 实例的更多想法,以下概述给出了地理 CS 的组成和属性的概念:
- pycrs.GeogCS
name
-> 字符串datum
-> 来自 pycrs.elements.datums 的数据name
-> pycrs.elements.datums.DatumNameproj4
-> 字符串esri_wkt
-> 字符串ogc_wkt
-> 字符串
ellips
-> 来自 pycrs.elements.ellipsoids 的椭球体name
-> pycrs.elements.ellipsoids.EllipsoidNameproj4
-> 字符串esri_wkt
-> 字符串ogc_wkt
-> 字符串
semimaj_ax
-> pycrs.elements.parameters.SemiMajorRadiusvalue
-> 浮动
semimin_ax
->(可选)pycrs.elements.parameters.SemiMinorRadiusvalue
-> 浮动
flat
->(可选)pycrs.elements.parameters.Flatningvalue
-> 浮动
inv_flat
->(可选)pycrs.elements.parameters.InverseFlatningvalue
-> 浮动
datumshift
->(可选)pycrs.elements.parameters.DatumShift 或无
prime_mer
-> pycrs.elements.parameters.PrimeMeridianvalue
-> 浮动
angunit
-> 来自 pycrs.elements.units 的角度单位,例如度数unitname
-> pycrs.elements.units.UnitNameproj4
-> 字符串esri_wkt
-> 字符串ogc_wkt
-> 字符串
unitmultiplier
-> pycrs.elements.units.UnitMultipliervalue
-> 浮动
twin_ax
-> 元组- 1:来自 pycrs.elements.directions 的命名指南针方向(东西方向)
- 2:来自 pycrs.elements.directions 的命名指南针方向(南北)
预计 CS
投影参考系统将坐标保存在投影 xy 空间中。除了通过子 GeogCS 定义地球的形状外,投影参考系统还定义了一些附加参数,以便将坐标转换为各种地图类型。我们以常用的 World Robinson 投影坐标系为例:
>>> crs = pycrs.parse.from_esri_code(54030)
如果它加载正确,您现在应该有一个 pycrs.ProjCS 实例,您还可以通过cs_type
属性检查它:
>>> isinstance(crs, pycrs.ProjCS)
True
>>> crs.cs_type
'Projected'
通过 ProjCS 实例,我们可以进一步访问它的子组件和参数。例如,如果我们想检查命名投影,我们可以这样做:
>>> proj = crs.proj
>>> isinstance(proj, pycrs.elements.projections.Robinson)
True
或者查看坐标单位的类型:
>>> unit = crs.unit
>>> isinstance(unit, pycrs.elements.units.Meter)
True
有关如何检查 ProjCS 实例的更多想法,以下概述给出了投影 CS 的组成和属性的想法:
- pycrs.ProjCS
name
-> 字符串geogcs
-> pycrs.GeogCS(参见地理 CS 部分...)proj
-> 来自 pycrs.elements.projections 的投影name
-> pycrs.elements.projections.ProjNameproj4
-> 字符串esri_wkt
-> 字符串ogc_wkt
-> 字符串
params
-> 列表- 1:来自 pycrs.elements.parameters 的命名参数
- 2:来自 pycrs.elements.parameters 的命名参数
- 3:...
- n:来自 pycrs.elements.parameters 的命名参数
unit
-> 来自 pycrs.elements.units 的单位unitname
-> pycrs.elements.units.UnitNameproj4
-> 字符串esri_wkt
-> 字符串ogc_wkt
-> 字符串
unitmultiplier
-> pycrs.elements.units.UnitMultipliervalue
-> 浮动
twin_ax
-> 元组- 1:来自 pycrs.elements.directions 的命名指南针方向(东西方向)
- 2:来自 pycrs.elements.directions 的命名指南针方向(南北)
转换为其他 CRS 格式
读取原始数据源的 crs 后,您可能希望将其转换为其他一些 crs 格式。PyCRS 允许转换为以下 CRS 格式:
转换为 Proj4
# as a string
>>> crs.to_proj4()
'+proj=robin +datum=WGS84 +ellps=WGS84 +a=6378137.0 +rf=298.257223563 +pm=0 +lon_0=0 +x_0=0 +y_0=0 +units=m +axis=enu +no_defs'
# or as a dict
>>> isinstance(crs.to_proj4(as_dict=True), dict)
True
转换为 ESRI WKT
>>> crs.to_esri_wkt()
'PROJCS["Unknown", GEOGCS["Unknown", DATUM["D_WGS_1984", SPHEROID["WGS_1984", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0], UNIT["Degree", 0.017453292519943295], AXIS["Lon", EAST], AXIS["Lat", NORTH]], PROJECTION["Robinson"], PARAMETER["Central_Meridian", 0], PARAMETER["False_Easting", 0], PARAMETER["False_Northing", 0], UNIT["Meter", 1.0], AXIS["X", EAST], AXIS["Y", NORTH]]'
转换为 OGC WKT
>>> crs.to_ogc_wkt()
'PROJCS["Unknown", GEOGCS["Unknown", DATUM["WGS_1984", SPHEROID["WGS_1984", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0], UNIT["degree", 0.017453292519943295], AXIS["Lon", EAST], AXIS["Lat", NORTH]], PROJECTION["Robinson"], PARAMETER["Central_Meridian", 0], PARAMETER["false_easting", 0], PARAMETER["false_northing", 0], UNIT["Meters", 1.0], AXIS["X", EAST], AXIS["Y", NORTH]]'
表示为坐标系代码
正如可以从各种机构定义的预定义坐标系统代码加载 crs 一样,有时还需要查找 crs 并将其表示为坐标系统代码(如果存在)。
表示为 EPSG 代码
目前,这仅用于查找 EPSG 代码。这将搜索 prj2epsg.org 上的 crs wkt 表示并返回第一个结果的 EPSG 代码,如果 wkt 没有关联的 EPSG 代码,则返回 None。
# to epsg code
>>> crs = pycrs.parse.from_epsg_code(4326)
>>> crs.to_epsg_code()
4326
因为 crs 定义可以有很多变化,查找它的坐标系统代码可能会产生多个可能的匹配,在这种情况下会引发警告。为了更好地控制选择正确的匹配,还可以使用 wkt 搜索实用程序返回所有可能的匹配元数据:
# more flexible utility search for epsg code
>>> crs = pycrs.parse.from_sr_code(42)
>>> results = pycrs.utils.wkt_to_epsg(crs.to_esri_wkt())
# display the top match info
>>> topmatch = results['codes'][0]
>>> topmatch['code']
'3857'
>>> topmatch['name']
'WGS 84 / Pseudo-Mercator'
食谱
修改 CS 实例
在大多数情况下,您只需要加载 CRS 并将其转换为某种格式。但是,有时您可能想要调整 CS 类型实例的参数。了解您的 CS 类型实例的组成,这就像设置/替换所需的属性一样简单。
让我们使用 World Robinson 投影来演示一些示例:
>>> crs = pycrs.parse.from_esri_code(54030) # Robinson projection from esri code
>>> crs.to_ogc_wkt()
'PROJCS["Unknown", GEOGCS["Unknown", DATUM["WGS_1984", SPHEROID["WGS_1984", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0], UNIT["degree", 0.017453292519943295], AXIS["Lon", EAST], AXIS["Lat", NORTH]], PROJECTION["Robinson"], PARAMETER["Central_Meridian", 0], PARAMETER["false_easting", 0], PARAMETER["false_northing", 0], UNIT["Meters", 1.0], AXIS["X", EAST], AXIS["Y", NORTH]]'
这是默认 Robinson 投影的地图:
假设我们想将其基准从 WGS84 切换到 NAD83,我们可以这样做:
>>> crs.geogcs.datum = pycrs.elements.datums.NAD83()
>>> crs.to_ogc_wkt()
'PROJCS["Unknown", GEOGCS["Unknown", DATUM["North_American_Datum_1983", SPHEROID["GRS_1980", 6378137.0, 298.257222101]], PRIMEM["Greenwich", 0], UNIT["degree", 0.017453292519943295], AXIS["Lon", EAST], AXIS["Lat", NORTH]], PROJECTION["Robinson"], PARAMETER["Central_Meridian", 0], PARAMETER["false_easting", 0], PARAMETER["false_northing", 0], UNIT["Meters", 1.0], AXIS["X", EAST], AXIS["Y", NORTH]]'
或者假设我们想改变它的本初子午线,使经度轴的中心更靠近太平洋而不是 Greenwhich:
>>> crs.geogcs.prime_mer.value = 160.0
>>> crs.to_ogc_wkt()
'PROJCS["Unknown", GEOGCS["Unknown", DATUM["North_American_Datum_1983", SPHEROID["GRS_1980", 6378137.0, 298.257222101]], PRIMEM["Greenwich", 160], UNIT["degree", 0.017453292519943295], AXIS["Lon", EAST], AXIS["Lat", NORTH]], PROJECTION["Robinson"], PARAMETER["Central_Meridian", 0], PARAMETER["false_easting", 0], PARAMETER["false_northing", 0], UNIT["Meters", 1.0], AXIS["X", EAST], AXIS["Y", NORTH]]'
这是该地图的外观(由于多边形穿过子午线,奇怪的线条只是渲染问题):
或者,如果我们只是一起切换投影类型:
>>> crs.proj = pycrs.elements.projections.Sinusoidal()
>>> crs.to_ogc_wkt()
'PROJCS["Unknown", GEOGCS["Unknown", DATUM["North_American_Datum_1983", SPHEROID["GRS_1980", 6378137.0, 298.257222101]], PRIMEM["Greenwich", 160], UNIT["degree", 0.017453292519943295], AXIS["Lon", EAST], AXIS["Lat", NORTH]], PROJECTION["Sinusoidal"], PARAMETER["Central_Meridian", 0], PARAMETER["false_easting", 0], PARAMETER["false_northing", 0], UNIT["Meters", 1.0], AXIS["X", EAST], AXIS["Y", NORTH]]'
坐标变换
想要在 CRS 格式之间进行转换的一个常见原因是,如果您想要将坐标从一个坐标系转换到另一个坐标系。在 Python 中,这通常使用 PyProj 模块完成,该模块仅采用 proj4 格式。使用 PyCRS,我们可以轻松定义要转换的原始坐标系并获得其 proj4 表示:
>>> fromcrs = pycrs.parse.from_epsg_code(4326) # WGS84 projection from epsg code
>>> fromcrs_proj4 = fromcrs.to_proj4()
然后,我们可以使用 PyCRS 从您选择的格式定义我们的目标投影,然后将其转换为 PyProj 期望的 proj4 格式:
>>> tocrs = pycrs.parse.from_esri_code(54030) # Robinson projection from esri code
>>> tocrs_proj4 = tocrs.to_proj4()
使用 proj4 crs 格式定义的源和目标投影,我们准备使用 PyProj 转换我们的数据坐标:
>>> import pyproj
>>> fromproj = pyproj.Proj(fromcrs_proj4)
>>> toproj = pyproj.Proj(tocrs_proj4)
>>> lng,lat = -76.7075, 37.2707 # Williamsburg, Virginia :)
>>> pyproj.transform(fromproj, toproj, lng, lat)
(-6766170.001635834, 3985755.032695593)
编写 Shapefile .prj 文件
转换数据坐标后,您可能还希望将数据与新 crs 一起保存回文件。使用 PyCRS,您可以以各种 crs 格式执行此操作。例如,要编写 shapefile .prj 文件:
>>> with open("testfiles/shapefile.prj", "w") as writer:
... _ = writer.write(tocrs.to_esri_wkt())
测试
测试套件仍在进行中,并且分布在多个文件中。文件 testdocs.py(官方文档测试)和 testbatch.py(测试和渲染一批投影)可以从提示符运行:
python testdocs.py
python testbatch.py
测试文件有一些依赖的 python 包,需要安装它们才能完全工作:
执照
此代码可根据 MIT 许可证免费共享、使用、重用和修改,请参阅 license.txt
学分
- 卡里姆·巴加特
- 迈卡·科克兰
- 迈克·基特里奇
- 罗杰·卢
- 格雷戈里·哈尔沃森
- 克拉克
- 大卫·霍斯
- 罗德里戈·E·普林西比