>>> from env_helper import info; info()

页面更新时间： 2022-12-27 23:57:52
运行环境：
    Linux发行版本: Debian GNU/Linux 11 (bullseye)
    操作系统内核: Linux-5.10.0-20-amd64-x86_64-with-glibc2.31
    Python版本: 3.9.2

11.2. 地图工具¶

geopandas提供了一个用于制作地图的高级接口matplotlib库。其映射形状与GeoSeries或GeoDataFrame使用plot（）方法一样便于操作。

一些示例数据：

>>> %matplotlib inline
>>> import geopandas as gpd
>>> world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
>>>
>>> cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))

我们现在可以绘制这些GeoDataFrames：

>>> # Examine country GeoDataFrame
>>> world.head()

	pop_est	continent	name	iso_a3	gdp_md_est	geometry
0	920938	Oceania	Fiji	FJI	8374.0	MULTIPOLYGON (((180.00000 -16.06713, 180.00000...
1	53950935	Africa	Tanzania	TZA	150600.0	POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...
2	603253	Africa	W. Sahara	ESH	906.5	POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3	35623680	North America	Canada	CAN	1674000.0	MULTIPOLYGON (((-122.84000 49.00000, -122.9742...
4	326625791	North America	United States of America	USA	18560000.0	MULTIPOLYGON (((-122.84000 49.00000, -120.0000...

>>> world.plot();

注意，一般来说，任何可以输出到matplotlib中pyplot的选项（或者适用于行的样式选项）都可以输出给plot（）方法。

11.2.1. Chloropleth地图¶

geopandas创建Chloropleth地图比较简单（每个形状颜色对应于相关联变量值的地图）。只需使用plot命令，将列参数设置为要用于分配颜色值的列即可。

>>> # Plot by GDP per capta
>>> world = world[(world.pop_est>0) & (world.name!="Antarctica")]
>>>
>>> world['gdp_per_cap'] = world.gdp_md_est / world.pop_est
>>>
>>> world.plot(column='gdp_per_cap');

选择颜色。还可以使用cmap选项修改plot颜色（有关色彩图的完整列表，请参阅matplotlib网站）：

>>> world.plot(column='gdp_per_cap', cmap='OrRd');

颜色映射的缩放方式也可以使用scheme选项（如果你已经安装了pysal，可以通过conda install pysal来实现）。默认情况下，scheme设置为“equal_intervals”，但也可以调整为任何其他pysal选项，如“分位数”，“百分位数”等。

要运行下面代码，需要先安装 mapclassify 模块：

pip install mapclassify

>>> world.plot(column='gdp_per_cap', cmap='OrRd', scheme='quantiles');

11.2.2. 地图图层¶

在合并地图之前，请一定要确保它们的共享公用CRS（以便它们对齐）。

>>> # Look at capitals
>>> # Note use of standard `pyplot` line style options
>>> cities.plot(marker='*', color='green', markersize=5);
>>>
>>> # Check crs
>>> cities = cities.to_crs(world.crs)
>>>
>>> # Now we can overlay over country outlines
>>> # And yes, there are lots of island capitals
>>> # apparently in the middle of the ocean!

方法1

>>> base = world.plot(color='white')
>>> cities.plot(ax=base, marker='o', color='red', markersize=5);

使用matplotlib对象

>>> import matplotlib.pyplot as plt
>>>
>>> fig, ax = plt.subplots()
>>>
>>> # set aspect to equal. This is done automatically
>>> # when using *geopandas* plot on it's own, but not when
>>> # working with pyplot directly.
>>>
>>> ax.set_aspect('equal')
>>> world.plot(ax=ax, color='white')
>>> cities.plot(ax=ax, marker='o', color='red', markersize=5)

<AxesSubplot:>

>>> fig, ax = plt.subplots()
>>>
>>> # set aspect to equal. This is done automatically
>>> # when using *geopandas* plot on it's own, but not when
>>> # working with pyplot directly.
>>>
>>> ax.set_aspect('equal')
>>>
>>> world.plot(ax=ax, color='white')
>>>
>>> cities.plot(ax=ax, marker='o', color='red', markersize=5)
>>>
>>> plt.show()

11.2.3. 管理投影¶

1.坐标参考系统

CRS非常重要，其意义在于GeoSeries或GeoDataFrame对象中的几何形状是任意空间中的坐标集合。 CRS可以使Python的坐标与地球上的具体实物点相关联。

使用proj4字符串代码的方法引用CRS。你可以从www.spatialreference.org或remotesensing.org找到最常用的代码。

例如，最常用的CRS之一是WGS84纬度 - 经度投影。这个投影的proj4表示通常可以用多种不同的方式引用相同CRS：

"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

通用的投影可以通过EPSG代码来引用，因此这个通用的投影也可以使用proj4字符串调用

"+init=epsg:4326".

geopandas包含有CRS的表示方法，包括proj4字符串本身

("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

或在字典中分解的参数：

{'proj': 'latlong', 'ellps': 'WGS84', 'datum': 'WGS84', 'no_defs': True}

此外，也可直接采用EPSG代码完成一些功能。

以下为参考值，是一些非常常见的投影和它们的proj4字符串：

WGS84纬度/经度：“+ proj = longlat + ellps = WGS84 + datum = WGS84 + no_defs”或“+ init = epsg：4326”
UTM区域（北）：“+ proj = utm + zone = 33 + ellps = WGS84 + datum = WGS84 + units = m + no_defs”
UTM Zones（南）：“+ proj = utm + zone = 33 + ellps = WGS84 + datum = WGS84 + units = m + no_defs + south”

2.设置投影

投影有两个操作方法：设置投影和重新投影。

当地理坐标有坐标数据（x-y值），但没有关于这些坐标如何指向实际位置的信息时，需要设置投影。如果没有设置CRS，地理信息的几何操作虽可以工作，但不会变换坐标，导出的文件也可能无法被其他软件正常理解。

注意，大多数情况你不需要设置投影。使用from_file（）命令加载的数据会始终包括投影信息。您可以通过crs属性来查看当前CRS对象：my_geoseries.crs。

然而，您可能会获得不包含投影的数据。在这种情况下，您必须设置CRS，以便地理数据库作出解释坐标的合理操作。

例如，将纬度和经度的电子表格手动转换为GeoSeries，可以通过WGS84纬度 - 经度CRS分配给crs属性的方法来设置投影：

my_geoseries.crs = {'init' :'epsg:4326'}

3.重新投影

重新投影是从一个坐标系变到另一个坐标系并重新表示位置的过程。地球上的任意一点放置到二维平面上时，所有投影都是失真的，您应用最好的投影可能不同于与您导入的数据相关联的投影。在这些情况下，可以使用to_crs命令重新进行投影：

>>> # load example data
>>> world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
>>>
>>> # Check original projection
>>> # (it's Platte Carre! x-y are long and lat)
>>> world.crs

<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
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

>>> # Visualize
>>> world.plot()

<AxesSubplot:>

>>> # Reproject to Mercator (after dropping Antartica)
>>> world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]
>>> world = world.to_crs({'init': 'epsg:3395'}) # world.to_crs(epsg=3395) would also work
>>> world.plot();

/usr/lib/python3/dist-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))

目录

上一个主题

下一个主题

常见问题

11.2. 地图工具¶

11.2.1. Chloropleth地图¶

11.2.2. 地图图层¶

11.2.3. 管理投影¶