投影#

坐标参考系统#

坐标参考系统 (CRS) 非常重要，因为 GeoSeries 或 GeoDataFrame 对象中的几何形状仅仅是一组在任意空间中的坐标。CRS 告诉 Python 这些坐标与地球上的位置之间的关系。

有关最常用投影的参考代码，请参见 spatialreference.org。

相同的坐标参考系统（CRS）通常可以用多种方式来表示。例如，最常用的CRS之一是WGS84纬度-经度投影。这可以通过权限代码 "EPSG:4326" 来表示。

GeoPandas可以接受pyproj.CRS.from_user_input()接受的任何内容：

CRS WKT 字符串
一个权威字符串（即“epsg:4326”）
一个 EPSG 整数代码（即 4326）
A pyproj.CRS
一个具有 to_wkt 方法的对象。
PROJ 字符串
PROJ参数字典
参数的PROJ关键字参数
带有PROJ参数的JSON字符串

供参考，一些非常常见的投影及其EPSG代码：

WGS84 纬度/经度: EPSG:4326
UTM区域（北）：EPSG:32633
UTM区域（南部）：EPSG:32733

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

通常，WKT或SRID比PROJ字符串更受欢迎，因为它们可以包含有关给定坐标参考系统（CRS）的更多信息。WKT和PROJ字符串之间的转换在大多数情况下会导致信息丢失，可能会导致错误的转换。如果可能，应该使用WKT2。

有关更多细节，请参阅描述坐标参考系统的最佳格式是什么。

设置投影#

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

设置投影可能是必要的，当某些原因导致GeoPandas拥有坐标数据（x-y值），但没有关于这些坐标如何对应于现实世界位置的信息时。设置投影是告诉GeoPandas如何解释坐标的方法。如果没有设置CRS，GeoPandas几何操作仍然可以工作，但坐标转换将无法进行，并且导出的文件可能无法被其他软件正确解释。

请注意，大多数情况下您不需要设置投影。从可靠来源加载的数据（使用 geopandas.read_file() 命令）应该始终包含投影信息。您可以通过 GeoSeries.crs 属性查看对象当前的 CRS。

然而，偶尔你可能会收到不包含投影的数据。在这种情况下，你必须设置CRS，以便GeoPandas知道如何解释坐标。

例如，如果您手动将经纬度的电子表格转换为GeoSeries，您可以通过将WGS84纬度-经度坐标参考系传递给GeoSeries.set_crs()方法（或通过设置GeoSeries.crs属性）来设置投影：

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

重投影#

重投影是将位置的表现从一个坐标系更改为另一个坐标系的过程。地球上位置的所有投影到二维平面的过程中都有失真。有关更多信息，请参见 Which projection is best。最适合您应用的投影可能与您导入的数据相关联的投影不同。在这些情况下，可以使用 GeoDataFrame.to_crs() 命令对数据进行重投影：

In [1]: import geodatasets

# load example data
In [2]: usa = geopandas.read_file(geodatasets.get_path('geoda.natregimes'))

# Check original projection
# (it's Plate Carrée! x-y are long and lat)
In [3]: usa.crs
Out[3]: 
<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 [4]: ax = usa.plot()

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

# Reproject to Albers contiguous USA
In [6]: usa = usa.to_crs("ESRI:102003")

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

In [8]: ax.set_title("NAD 1983 Albers contiguous USA");

多个几何列的投影#

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的操作都保留投影；然而，如果您遍历包含几何的列，这些信息可能会丢失。

升级到 GeoPandas 0.7 与 pyproj > 2.2 和 PROJ > 6#

从 GeoPandas 0.7 开始，.crs 属性的 GeoSeries 或 GeoDataFrame 存储 CRS 信息为一个 pyproj.CRS，不再是 proj4 字符串或字典。

之前，你可能见过这样的内容：

>>> 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的改进，但可能也需要对您的代码进行一些更改。有关更多信息，请参见这篇博客文章。下面的各个小节涵盖了不同的迁移问题。

有关更多信息，请参见pyproj 文档中关于pyproj.CRS对象的内容。

从文件导入数据#

当使用 geopandas.read_file() 读取地理空间文件时，大多数情况下应该能够开箱即用。例如，读取示例国家数据集会产生一个正确的 CRS：

In [9]: df = geopandas.read_file(geodatasets.get_path('naturalearth.land'))

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对象可能与预期不完全相符。请参阅下面的部分以了解可能的原因以及如何解决。

手动指定坐标参考系统#

在代码中手动指定CRS时（例如，因为您的数据尚未有CRS，或者在转换为另一个CRS时），这可能需要更改您的代码。

“初始化” proj4 字符串/字典

目前，有很多人（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 字符串/字典

尽管完整的 proj4 字符串没有被弃用（与上面的“init”字符串相反），但如果可能，仍然建议用 EPSG 代码进行更改。

例如，如果您知道您正在使用的投影的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"

找出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

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

此外，在像 Spatial Reference 和 epsg.org 这样的网站上可以找到许多CRS的描述，包括它们的EPSG代码和proj4字符串定义。

其他格式

除了上面提到的EPSG代码，还有其他方式来指定坐标参考系（CRS）：一个实际的 pyproj.CRS 对象，一个WKT字符串，一个PROJ JSON字符串，等等。任何被 pyproj.CRS.from_user_input() 接受的内容都可以指定给 crs 关键字/属性在GeoPandas中。

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

CRS的轴顺序#

从PROJ 6 / pyproj 2开始，官方EPSG定义的轴顺序得到了尊重。例如，在使用地理坐标（经度和纬度的度数）时，标准EPSG:4326的坐标参考系统将如下所示：

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

这提到的顺序是 (lat, lon)，因为这是 EPSG:4326 中坐标的官方顺序。然而，在 GeoPandas 中，坐标始终以 (x, y) 的形式存储，因此顺序为 (lon, lat)，无论 CRS 如何（即在 GIS 中使用的“传统”顺序）。在重新投影时，GeoPandas 和 pyproj 会在底层处理这一轴顺序差异，因此用户无需担心这个。

为什么我的CRS没有正确识别？#

在“野外”中，有许多文件源和CRS定义，可能有一个CRS描述并不完全符合新的PROJ > 6标准（proj4字符串，旧的WKT格式等）。在这种情况下，您将获得一个pyproj.CRS对象，可能并不完全符合您的预期（例如，与预期的EPSG代码不相等）。以下是一些可能情况的列表。

我得到一个“边界坐标参考系统”？#

一些 CRS 定义包括一个 “towgs84” 条款，这可能会在识别实际 CRS 时出现问题。

例如，EPSG:31370（比利时使用的本地投影）的proj4和WKT表示形式都包含此条款。将该网站上的其中一个定义提取出来，并创建一个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

你会注意到上面并不是预期中的“投影坐标参考系统”，而是“绑定坐标参考系统”。这是因为它“绑定”了转换到WGS84，并且在重新投影时将始终使用这个，而不是让PROJ决定最佳的转换。

要获取实际的基础投影坐标参考系统（CRS），您可以使用 .source_crs 属性：

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

现在您有了一个“投影坐标参考系”，现在它也会识别正确的EPSG编号：

>>> crs.to_epsg()

>>> crs.source_crs.to_epsg()
31370

我有不同的轴顺序？#

如上所述，pyproj 现在尊重 EPSG 定义的坐标轴顺序。然而，proj4 字符串或旧版 WKT 并没有正确指定这一点，这可能是 CRS 对象与预期的 EPSG 代码不相等的原因。

考虑以下一个加拿大投影的参考坐标系 EPSG:2953 的例子。当从 WKT 字符串构建 CRS 对象时，正如在 EPSG: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 对象具有“东坐标，北坐标”（即 x，y）轴顺序，而从 EPSG 代码构建的 CRS 对象具有（北坐标，东坐标）轴顺序。

仅仅拥有这种轴顺序的差异在使用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']

这将不再有效。要从一个 crs 对象获取 EPSG 代码，您可以使用 to_epsg() 方法。

或者检查一个CRS是否为某个UTM区域：

'+proj=utm ' in gdf.crs

可以用更强大的检查替代（需要 pyproj 2.6+）：

gdf.crs.utm_zone is not None

此外，在pyproj.CRS类上还有许多其他方法可用于获取有关CRS的信息。

投影#