问题描述
我正在根据某些值在cartopy地图上为国家着色。我正在使用来自https://www.naturalearthdata.com/
的geopandas和shapefile在遍历数据框df以获取某些国家/地区的几何图形时,我遇到了一个问题。我可以得到具有Multipolygon几何体的国家/地区的几何图形,但是我无法获得具有多边形几何体的国家/地区的几何图形,例如比利时或奥地利。
这是我的代码:
#imports
import matplotlib.pyplot as plt
import matplotlib
import cartopy
from cartopy.io import shapereader
import cartopy.crs as ccrs
import geopandas
import numpy as np
# get natural earth data (http://www.naturalearthdata.com/)
# get country borders
resolution = '10m'
category = 'cultural'
name = 'admin_0_countries'
shpfilename = shapereader.natural_earth(resolution,category,name)
# read the shapefile using geopandas
df = geopandas.read_file(shpfilename)
# Set up the canvas
fig = plt.figure(figsize=(20,20))
central_lon,central_lat = 0,45
extent = [-10,28,35,65]
ax = plt.axes(projection=cartopy.crs.Orthographic(central_lon,central_lat))
ax.set_extent(extent)
ax.gridlines()
# Add natural earth features and borders
ax.add_feature(cartopy.feature.BORDERS,linestyle='-',alpha=0.8)
ax.add_feature(cartopy.feature.OCEAN,facecolor=("lightblue"))
ax.add_feature(cartopy.feature.LAND,facecolor=("lightgreen"),alpha=0.35)
ax.coastlines(resolution='10m')
# Countries and value
countries = ['Sweden','Netherlands','Ireland','Denmark','Germany','Greece','France','Spain','Portugal','Italy','United Kingdom','Austria']
value = [47.44,32.75,27.53,23.21,20.08,18.08,17.23,13.59,12.13,5.66,22.43,7]
# normalise values
value_norm = (value-np.nanmin(value))/(np.nanmax(value) - np.nanmin(value))
# Colourmap
cmap = matplotlib.cm.get_cmap('YlOrBr')
for country,value_norm in zip(countries,value_norm):
# read the borders of the country in this loop
poly = df.loc[df['ADMIN'] == country]['geometry'].values[0]
# get the color for this country
rgba = cmap(value_norm)
# plot the country on a map
ax.add_geometries(poly,crs=ccrs.PlateCarree(),facecolor=rgba,edgecolor='none',zorder=1)
# Add a scatter plot of the original data so the colorbar has the correct numbers
dummy_scat = ax.scatter(value,value,c=value,cmap=cmap,zorder=0)
fig.colorbar(mappable=dummy_scat,label='Percentage of Do and Dont`s [%]',orientation='horizontal',shrink=0.8)
plt.show()
fig.savefig("Länderübersicht.jpg")
我该如何遍历这些国家或为其着色,还是必须获取另一个shapefile? 谢谢!
解决方法
命令ax.add_geometries()
要求提供几何列表,因此,单个几何将导致错误。要修复您的代码,您可以替换以下行:
ax.add_geometries(poly,crs=ccrs.PlateCarree(),facecolor=rgba,edgecolor='none',zorder=1)
使用以下代码行:
# plot the country on a map
if poly.geom_type=='MultiPolygon':
# `poly` is a list of geometries
ax.add_geometries(poly,zorder=1)
elif poly.geom_type=='Polygon':
# `poly` is a geometry
# Austria,Belgium
# Plot it `green` for checking purposes
ax.add_geometries([poly],facecolor="green",zorder=1)
else:
pass #do not plot the geometry
请注意,如果poly.geom_type
是'Polygon',我只用[poly]代替poly
。
从错误代码TypeError: 'Polygon' object is not iterable
中汲取灵感,我从一个假设开始,即我们需要某种可迭代的对象,例如多边形列表。从this answer进行绘制,我发现函数shapely.geometry.MultiPolygon
可以完成这项工作。您只需将多边形列表传递给它即可。添加一些逻辑以仅在检测到Polygon
而不是MultiPolygon
并且我们有以下情况时才执行此操作:
poly = df.loc[df['ADMIN'] == country]['geometry'].values[0]
if type(poly) == shapely.geometry.polygon.Polygon:
simple_poly = df.loc[df['ADMIN'] == country]['geometry'].values[0]
list_polys = [poly,poly]
poly = shapely.geometry.MultiPolygon(list_polygons)
这是一个非常棘手的解决方案,它将打印两次多边形,因此请注意,如果您以后决定使其透明或类似。或者,您可以使用[poly,poly]
代替[poly,some_other_poly_outside_map_area]
。