使用 geopandas 和 matplotlib 绘制地图
Plotting a map using geopandas and matplotlib
我有一个小的 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()
尝试df.unary_union。该函数会将点聚合到一个几何图形中。
Jupyter 笔记本 can plot it
GeoPandas 文档包含有关如何向地图添加背景的示例 (https://geopandas.readthedocs.io/en/latest/gallery/plotting_basemap_background.html)，下面将对此进行更详细的说明。
您将不得不处理 tiles，即通过网络服务器提供的 (png) 图像，URL 喜欢
http://.../Z/X/Y.png, where Z is the zoom level, and X and Y identify the tile
geopandas 的文档展示了如何将图块设置为地块的背景、获取正确的图块以及完成所有其他困难的空间同步工作等...
安装
假设 GeoPandas 已经安装，您需要 contextily package in addition. If you are under windows, you may want to pick a look at
用例
创建一个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 来为地图找到合适的分辨率。 例如/即。 :
... 此类决议暗中要求更改 x/y 限制。
只想添加有关缩放的用例，根据新的 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
我有一个小的 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()
尝试df.unary_union。该函数会将点聚合到一个几何图形中。 Jupyter 笔记本 can plot it
GeoPandas 文档包含有关如何向地图添加背景的示例 (https://geopandas.readthedocs.io/en/latest/gallery/plotting_basemap_background.html)，下面将对此进行更详细的说明。
您将不得不处理 tiles，即通过网络服务器提供的 (png) 图像，URL 喜欢
http://.../Z/X/Y.png, where Z is the zoom level, and X and Y identify the tile
geopandas 的文档展示了如何将图块设置为地块的背景、获取正确的图块以及完成所有其他困难的空间同步工作等...
安装
假设 GeoPandas 已经安装，您需要 contextily package in addition. If you are under windows, you may want to pick a look at
用例
创建一个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 来为地图找到合适的分辨率。 例如/即。 :
... 此类决议暗中要求更改 x/y 限制。
只想添加有关缩放的用例，根据新的 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