使用 geopandas 和 matplotlib 绘制地图

Posted 2023-05-07

【中文标题】使用 geopandas 和 matplotlib 绘制地图

【英文标题】：Plotting a map using geopandas and matplotlib

【发布时间】：2019-01-08 09:33:24

【问题描述】：

我有一个小型 csv，其中包含来自英格兰伯明翰的 6 个坐标。我用 pandas 读取 csv，然后将其转换为 GeoPandas DataFrame，用 Shapely Points 更改我的纬度和经度列。我现在正在尝试绘制我的 GeoDataframe，我能看到的只是点。我如何获得伯明翰地图的代表？非常感谢 GeoPandas 上的良好文档来源。

from shapely.geometry import Point
import geopandas as gpd
import pandas as pd

df = pd.read_csv('SiteLocation.csv')
df['Coordinates'] = list(zip(df.LONG, df.LAT))
df['Coordinates'] = df['Coordinates'].apply(Point)
# Building the GeoDataframe 
geo_df = gpd.GeoDataFrame(df, geometry='Coordinates')
geo_df.plot()  

【问题讨论】：

@Kanak 你能把它添加为答案吗？

@joris 添加为答案。

【参考方案1】：

GeoPandas 文档包含一个关于如何为地图添加背景的示例 (https://geopandas.readthedocs.io/en/latest/gallery/plotting_basemap_background.html)，下面将对此进行更详细的说明。

---

您必须处理tiles，即通过网络服务器提供的 (png) 图像，其 URL 类似于

http://.../Z/X/Y.png，其中 Z 是缩放级别，X 和 Y 标识图块

geopandas 的文档展示了如何将图块设置为绘图的背景、获取正确的图块以及完成空间同步等所有其他困难的工作......

---

安装

假设 GeoPandas 已经安装，你还需要 contextily 包。如果你是windows下的，不妨看看How to install Contextily?

用例

创建一个python脚本并定义contextily helper function

import contextily as ctx

def add_basemap(ax, zoom, url='http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'):
    xmin, xmax, ymin, ymax = ax.axis()
    basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=zoom, url=url)
    ax.imshow(basemap, extent=extent, interpolation='bilinear')
    # restore original x/y limits
    ax.axis((xmin, xmax, ymin, ymax))

然后玩

import matplotlib.pyplot as plt
from shapely.geometry import Point
import geopandas as gpd
import pandas as pd

# Let's define our raw data, whose epsg is 4326
df = pd.DataFrame(
    'LAT'  :[-22.266415, -20.684157],
    'LONG' :[166.452764, 164.956089],
)
df['coords'] = list(zip(df.LONG, df.LAT))

# ... turn them into geodataframe, and convert our
# epsg into 3857, since web map tiles are typically
# provided as such.
geo_df = gpd.GeoDataFrame(
    df, crs  ='init': 'epsg:4326',
    geometry = df['coords'].apply(Point)
).to_crs(epsg=3857)

# ... and make the plot
ax = geo_df.plot(
    figsize= (5, 5),
    alpha  = 1
)
add_basemap(ax, zoom=10)
ax.set_axis_off()
plt.title('Kaledonia : From Hienghène to Nouméa')
plt.show()

---

注意：您可以使用zoom 来找到适合地图的分辨率。 例如/I.e. ：

...此类决议隐含地要求更改 x/y 限制。

【讨论】：

有任何问题@Herc01 吗？

感谢您的回答。关于缩放级别的最后一句话是什么意思？如果更改缩放级别，通常不需要手动编辑 x/y 限制。

@joris，我的意思是使用为参数设置的值zoom 在某种程度上等同于更改地图的（图像）分辨率。你同意吗？

是的，这正是“缩放”的意思。这可能有点令人困惑，因为地图本身没有放大（查看较小的区域），但较低/较高的 zoom 会放大整个区域，从而导致较低/较高的分辨率。

对不起，我打算在编辑中添加评论，为什么我删除它，但忘记了：我知道它肯定很难安装，但我认为这里会分散实际答案的注意力（而且它也不是唯一可能的安装方式，例如对于使用 conda 的人来说，这将是错误的方式）。您想提出一个新问题“如何根据上下文进行安装？”，然后在该部分添加您的答案。然后我们可以从这里链接到那个问题（并将这个答案集中在如何）。

【参考方案2】：

只想添加有关缩放的用例，根据新的xlim 和ylim 坐标更新底图。我想出的解决方案是：

首先在ax 上设置可以检测xlim_changed 和ylim_changed 的回调 一旦检测到两者都已更改，则获取新的plot_area，调用ax.get_xlim() 和ax.get_ylim() 然后清除ax 并重新绘制底图和任何其他数据

显示首都的世界地图示例。当您放大地图时，您会注意到正在更新的分辨率。

import geopandas as gpd
import matplotlib.pyplot as plt
import contextily as ctx


figsize = (12, 10)
osm_url = 'http://tile.stamen.com/terrain/z/x/y.png'
EPSG_OSM = 3857
EPSG_WGS84 = 4326

class MapTools:
    def __init__(self):
        self.cities = gpd.read_file(
            gpd.datasets.get_path('naturalearth_cities'))
        self.cities.crs = EPSG_WGS84
        self.cities = self.convert_to_osm(self.cities)

        self.fig, self.ax = plt.subplots(nrows=1, ncols=1, figsize=figsize)
        self.callbacks_connect()

        # get extent of the map for all cities
        self.cities.plot(ax=self.ax)
        self.plot_area = self.ax.axis()

    def convert_to_osm(self, df):
        return df.to_crs(epsg=EPSG_OSM)

    def callbacks_connect(self):
        self.zoomcallx = self.ax.callbacks.connect(
            'xlim_changed', self.on_limx_change)
        self.zoomcally = self.ax.callbacks.connect(
            'ylim_changed', self.on_limy_change)

        self.x_called = False
        self.y_called = False

    def callbacks_disconnect(self):
        self.ax.callbacks.disconnect(self.zoomcallx)
        self.ax.callbacks.disconnect(self.zoomcally)

    def on_limx_change(self, _):
        self.x_called = True
        if self.y_called:
            self.on_lim_change()

    def on_limy_change(self, _):
        self.y_called = True
        if self.x_called:
            self.on_lim_change()

    def on_lim_change(self):
        xlim = self.ax.get_xlim()
        ylim = self.ax.get_ylim()
        self.plot_area = (*xlim, *ylim)
        self.blit_map()

    def add_base_map_osm(self):
        if abs(self.plot_area[1] - self.plot_area[0]) < 100:
            zoom = 13

        else:
            zoom = 'auto'

        try:
            basemap, extent = ctx.bounds2img(
                self.plot_area[0], self.plot_area[2],
                self.plot_area[1], self.plot_area[3],
                zoom=zoom,
                url=osm_url,)
            self.ax.imshow(basemap, extent=extent, interpolation='bilinear')

        except Exception as e:
            print(f'unable to load map: e')

    def blit_map(self):
        self.ax.cla()
        self.callbacks_disconnect()
        cities = self.cities.cx[
            self.plot_area[0]:self.plot_area[1],
            self.plot_area[2]:self.plot_area[3]]
        cities.plot(ax=self.ax, color='red', markersize=3)

        print('*'*80)
        print(self.plot_area)
        print(f'len(cities) cities in plot area')

        self.add_base_map_osm()
        self.callbacks_connect()

    @staticmethod
    def show():
        plt.show()


def main():
    map_tools = MapTools()
    map_tools.show()

if __name__ == '__main__':
    main()

在 Linux Python3.8 上运行并安装以下 pip

affine==2.3.0
attrs==19.3.0
autopep8==1.4.4
Cartopy==0.17.0
certifi==2019.11.28
chardet==3.0.4
Click==7.0
click-plugins==1.1.1
cligj==0.5.0
contextily==1.0rc2
cycler==0.10.0
descartes==1.1.0
Fiona==1.8.11
geographiclib==1.50
geopandas==0.6.2
geopy==1.20.0
idna==2.8
joblib==0.14.0
kiwisolver==1.1.0
matplotlib==3.1.2
mercantile==1.1.2
more-itertools==8.0.0
munch==2.5.0
numpy==1.17.4
packaging==19.2
pandas==0.25.3
Pillow==6.2.1
pluggy==0.13.1
py==1.8.0
pycodestyle==2.5.0
pyparsing==2.4.5
pyproj==2.4.1
pyshp==2.1.0
pytest==5.3.1
python-dateutil==2.8.1
pytz==2019.3
rasterio==1.1.1
requests==2.22.0
Rtree==0.9.1
Shapely==1.6.4.post2
six==1.13.0
snuggs==1.4.7
urllib3==1.25.7
wcwidth==0.1.7

特别注意contextily==1.0rc2的要求

在 Windows 上，我使用 Conda (P3.7.3) 并且不要忘记设置用户变量：

GDAL c:\Users\<username>\Anaconda3\envs\<your environment>\Library\share\gdal

PROJLIB c:\Users\<username>\Anaconda3\envs\<your environment>\Library\share

【参考方案3】：

试试df.unary_union。该函数会将点聚合成一个几何图形。 Jupyter笔记本can plot it

【讨论】：

对不起，伙计——那不是我想要的。不过谢谢。你看看其他一些未解决的问题