管理项目#

坐标参考系#

坐标参考系(CRS)很重要,因为GeoSeries或GeoDataFrame对象中的几何形状只是任意空间中的坐标集合。CRS告诉Python这些坐标如何与地球上的位置相关。

您可以在以下位置找到最常用投影的代码 www.spatialreference.org

相同的CRS通常可以在许多方面被引用。例如,最常用的CRS之一是WGS84经纬度投影。这可以使用授权代码来参考 "EPSG:4326"

地貌熊猫 可以接受任何接受的东西 pyproj.CRS.from_user_input()

  • CRS WKT字符串

  • 授权字符串(即“epsg:4326”)

  • EPSG整数代码(如4326)

  • A pyproj.CRS

  • 具有to_wkt方法的对象。

  • 项目字符串

  • 项目参数词典

  • 参数的Proj关键字实参

  • 带有项目参数的JSON字符串

以下是几个非常常见的投影及其EPSG代码,以供参考:

  • WGS84纬度/经度: "EPSG:4326"

  • UTM区域(北部): "EPSG:32633"

  • UTM区域(南部): "EPSG:32733"

存储CRS信息的最佳格式是什么?#

通常,WKT或SRID比项目字符串更可取,因为它们可以包含有关给定CRS的更多信息。在大多数情况下,WKT和Proj字符串之间的转换将导致信息丢失,可能会导致错误的转换。如果可能,应使用WKT2。

欲了解更多详情,请访问https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems

设置投影#

投影有两个相关的操作:设置投影和重新投影。

出于某些原因,可能需要设置投影 地貌熊猫 具有坐标数据(x-y值),但没有关于这些坐标如何引用真实世界中的位置的信息。设置投影是一个人告诉你的方式 地貌熊猫 如何解释坐标。如果没有设置CRS, 地貌熊猫 几何图形操作仍可用,但不能进行坐标变换,并且其他软件可能无法正确解释导出的文件。

请注意, 大部分时间里 您不必设置投影。从信誉良好的来源加载的数据(使用 geopandas.read_file() 命令) 应该 始终包含投影信息。您可以通过查看对象当前的CRS GeoSeries.crs 属性。

但是,有时您可能会获得不包括投影的数据。在这种情况下,您必须将CRS设置为 地貌熊猫 知道如何解释坐标。

例如,如果手动将纬度和经度的电子表格转换为GeoSeries,则可以通过将WGS84经纬度CRS传递给 GeoSeries.set_crs() 方法(或通过设置 GeoSeries.crs 属性):

my_geoseries = my_geoseries.set_crs("EPSG:4326")
my_geoseries = my_geoseries.set_crs(epsg=4326)

重新投射#

重新投影是将位置的表示从一个坐标系更改到另一个坐标系的过程。将地球上的所有位置投影到一个二维平面中 are distortions ,则最适合您的应用程序的投影可能不同于与您导入的数据相关联的投影。在这些情况下,可以使用 GeoDataFrame.to_crs() 命令:

# load example data
In [1]: world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))

# Check original projection
# (it's Plate Carrée! x-y are long and lat)
In [2]: world.crs
Out[2]: 
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# Visualize
In [3]: ax = world.plot()

In [4]: ax.set_title("WGS84 (lat/lon)");

# Reproject to Mercator (after dropping Antartica)
In [5]: world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]

In [6]: world = world.to_crs("EPSG:3395") # world.to_crs(epsg=3395) would also work

In [7]: ax = world.plot()

In [8]: ax.set_title("Mercator");
../../_images/world_starting.png ../../_images/world_reproj.png

多个几何图形柱的投影#

GeoPandas 0.8实现了对分配给同一GeoDataFrame的不同几何列的不同投影的支持。现在,投影与每列的几何体一起存储(直接在GeometryArray层级上)。

请注意,如果GeometryArray具有指定的投影,则在创建GeoSeries或GeoDataFrame期间,它不能被另一个不一致的投影覆盖:

>>> array.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
...
>>> GeoSeries(array, crs=4326)  # crs=4326 is okay, as it matches the existing CRS
>>> GeoSeries(array, crs=3395)  # crs=3395 is forbidden as array already has CRS
ValueError: CRS mismatch between CRS of the passed geometries and 'crs'. Use 'GeoSeries.set_crs(crs, allow_override=True)' to overwrite CRS or 'GeoSeries.to_crs(crs)' to reproject geometries.
    GeoSeries(array, crs=3395).crs

如果要覆盖投影,则可以手动将其分配给GeoSeries,或使用以下任一方法将几何重新投影到目标投影 GeoSeries.set_crs(epsg=3395, allow_override=True)GeoSeries.to_crs(epsg=3395)

所有基于GeometryArray的操作都会保留投影;但是,如果循环遍历包含几何图形的列,则此信息可能会丢失。

使用pyproj>2.2和proj>6升级到GeoPandas 0.7#

从GeoPandas 0.7开始, .crs 属性将CRS信息存储为 pyproj.CRS ,而不再作为项目字符串或词典。

以前,你可能见过这样的场景:

>>> gdf.crs
{'init': 'epsg:4326'}

而现在,您将看到类似以下内容:

>>> gdf.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
...
>>> type(gdf.crs)
pyproj.crs.CRS

这提供了更好的用户界面,并集成了对pyproj和proj 6的改进,但也可能需要对代码进行一些更改。看见 this blogpost 了解更多背景信息,下面的小节介绍了不同可能的迁移问题。

请参阅 pyproj docs 了解更多有关 pyproj.CRS 对象。

从文件导入数据#

使用读取地理空间文件时 geopandas.read_file() ,事情基本上应该是开箱即用的。例如,阅读示例国家/地区数据集将生成正确的CRS:

In [9]: df = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))

In [10]: df.crs
Out[10]: 
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

但是,在某些情况下(使用较旧的CRS格式),生成的CRS对象可能不完全像预期的那样。请参阅 section below 可能的原因以及如何解决它。

手动指定CRS#

在代码中手动指定CRS时(例如,因为您的数据还没有CRS,或者在转换为另一个CRS时),这可能需要更改您的代码。

"init" proj4 strings/dicts

目前,很多人(之前的GeoPandas文档也显示了这一点)使用“init”proj4字符串指定EPSG代码:

## OLD
GeoDataFrame(..., crs={'init': 'epsg:4326'})
# or
gdf.crs = {'init': 'epsg:4326'}
# or
gdf.to_crs({'init': 'epsg:4326'})

上面的代码现在将从pyproj中引发一个弃用警告,您应该只使用EPSG代码本身,而不是“init”proj4字符串,如下所示:

## NEW
GeoDataFrame(..., crs="EPSG:4326")
# or
gdf.crs = "EPSG:4326"
# or
gdf.to_crs("EPSG:4326")

proj4 strings/dicts

尽管不建议使用完整的proj4字符串(与上面的“init”字符串相反),但如果可能,仍建议使用EPSG代码对其进行更改。

例如,不是:

gdf.crs = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

我们建议执行以下操作:

gdf.crs = "EPSG:2163"

if 您知道您正在使用的投影的EPSG代码。

找出EPSG代码的一种可能方法是使用pyproj来执行此操作:

>>> import pyproj
>>> crs = pyproj.CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")
>>> crs.to_epsg()
2163

(您可能需要设置 min_confidence 的关键字 to_epsg 如果匹配不完美,则设置为较低的值)

此外,在像这样的网站上 spatialreference.orgepsg.io 可以找到许多CRS的描述,包括它们的EPSG代码和Proj4字符串定义。

其他格式

除了上面提到的EPSG代码之外,还有其他指定CRS的方法:实际的 pyproj.CRS 对象、WKT字符串、项目JSON字符串等。 pyproj.CRS.from_user_input() 可以通过指定给 crs GeoPandas中的关键字/属性。

也兼容CRS对象,例如来自 rasterio 包,可以直接传递给GeoPandas。

CRS的轴序#

从PROJ6/PYPROJ2开始,遵循正式EPSG定义的轴顺序。例如,在标准EPSG:4326中使用地理坐标(经度和纬度)时,CRS将如下所示:

>>> pyproj.CRS(3EPSG:4326")
<Geographic 2D CRS: EPSG:4326>
...
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
...

这提到的顺序为(经度,经度),因为这是EPSG:4326中的官方坐标顺序。然而,在GeoPandas中,坐标始终存储为(x,y),因此存储为(经度、经度)顺序,而不考虑CRS(即,在地理信息系统中使用的“传统”顺序)。在重新投影时,GeoPandas和pyproj将在幕后处理轴顺序的差异,因此用户不需要关心这一点。

为什么它不能正确识别我的CRS?#

有许多文件源和CRS定义可能不完全符合Proj>6的新标准(Proj4字符串、较旧的WKT格式,等等)。在这种情况下,您将获得一个 pyproj.CRS 对象可能不完全符合您的预期(例如,不等于预期的EPSG代码)。下面我们列出几个可能的案例。

我得到了一个“绑定CRS”?#

某些CRS定义包括 “Towgs84”条款 ,这可能会给识别实际的CRS带来问题。

例如,EPSG:31370(比利时使用的本地投影)的Proj4和WKT表示可在 https://spatialreference.org/ref/epsg/31370/ 包括这个。从该站点获取这些定义之一并创建CRS对象时:

>>> import pyproj
>>> crs = pyproj.CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m +no_defs")
>>> crs
<Bound CRS: +proj=lcc +lat_1=51.16666723333333 +lat_2=49.83333 ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: Transformation from unknown to WGS84
- method: Position Vector transformation (geog2D domain)
Datum: Unknown based on International 1909 (Hayford) ellipsoid
- Ellipsoid: International 1909 (Hayford)
- Prime Meridian: Greenwich
Source CRS: unknown

您会注意到,上面并不是一个预期的“投影CRS”,而是一个“绑定CRS”。这是因为它“绑定”到WGS84的转换,并且在重新投影时将始终使用它,而不是让Proj确定最佳转换。

若要获取实际的基础投影CRS,可以使用 .source_crs 属性:

>>> crs.source_crs
<Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unk ...>
Name: unknown
...

现在我们有了“Projected CRS”,现在它还将识别正确的EPSG编号:

>>> crs.to_epsg()

>>> crs.source_crs.to_epsg()
31370

我有一个不同的轴顺序?#

如上所述,pyproj现在遵循EPSG定义的轴顺序。然而,proj4字符串或较早的WKT版本没有正确指定这一点,这可能是CRS对象不等于预期的EPSG代码的一个原因。

请考虑以下加拿大计划CRS“EPSG:2953”的示例。从上提供的WKT字符串构造CRS对象时 https://epsg.io/2953

>>> crs = pyproj.CRS("""PROJCS["NAD83(CSRS) / New Brunswick Stereographic",
...     GEOGCS["NAD83(CSRS)",
...         DATUM["NAD83_Canadian_Spatial_Reference_System",
...             SPHEROID["GRS 1980",6378137,298.257222101,
...                 AUTHORITY["EPSG","7019"]],
...             AUTHORITY["EPSG","6140"]],
...         PRIMEM["Greenwich",0,
...             AUTHORITY["EPSG","8901"]],
...         UNIT["degree",0.0174532925199433,
...             AUTHORITY["EPSG","9122"]],
...         AUTHORITY["EPSG","4617"]],
...     PROJECTION["Oblique_Stereographic"],
...     PARAMETER["latitude_of_origin",46.5],
...     PARAMETER["central_meridian",-66.5],
...     PARAMETER["scale_factor",0.999912],
...     PARAMETER["false_easting",2500000],
...     PARAMETER["false_northing",7500000],
...     UNIT["metre",1,
...         AUTHORITY["EPSG","9001"]],
...     AUTHORITY["EPSG","2953"]]""")

>>> crs
<Projected CRS: PROJCS["NAD83(CSRS) / New Brunswick Stereographic" ...>
Name: NAD83(CSRS) / New Brunswick Stereographic
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
...

尽管这是在线找到的“EPSG:2953”的WKT字符串,但此CRS对象的计算结果不等于此EPSG代码:

>>> crs == "EPSG:2953"
False

如果我们从EPSG代码(截断输出)构造CRS对象:

>>> pyproj.CRS("EPSG:2953")
<Projected CRS: EPSG:2953>
Name: NAD83(CSRS) / New Brunswick Stereographic
Axis Info [cartesian]:
- N[north]: Northing (metre)
- E[east]: Easting (metre)
...

您可以看到,从WKT字符串构造的CRS对象的轴顺序为“Easting,Northing”(即x,y),而从EPSG代码构造的CRS对象的轴顺序为(Northing,Eating)。

在GeoPandas中使用CRS时,只有轴顺序上的这种差异是没有问题的,因为GeoPandas总是使用(x,y)顺序来存储数据,而不管CRS定义如何。但是,您可能仍然希望验证它是否等同于预期的EPSG代码。通过降低 min_confidence ,则将忽略轴顺序:

>>> crs.to_epsg()

>>> crs.to_epsg(min_confidence=20)
2953

这个 .crs 属性不再是词典或字符串#

如果您依赖于 .crs 对象是字典或字符串,这样的代码可以被破解,因为它现在是 pyproj.CRS 对象。但该对象实际上提供了一个更健壮的接口来获取有关CRS的信息。

例如,如果您使用以下代码获取EPSG代码:

gdf.crs['init']

这将不再起作用。获取EPSG代码 crs 对象,则可以使用 to_epsg() 方法。

或检查CRS是否为某个UTM区域:

'+proj=utm ' in gdf.crs

可替换为更强大的检查(需要pyproj 2.6+):

gdf.crs.utm_zone is not None

上还提供了许多其他方法。 pyproj.CRS 类以获取有关CRS的信息。