使用 Cartopy 合并国家

Question

我正在使用以下代码制作瑞典、挪威和芬兰作为一个区域的地图。但是，我正在为此苦苦挣扎。我正在关注这个例子，Python Matplotlib Cartopy Color One Country 中的映射。

from shapely.geometry import Polygon
from cartopy.io import shapereader
import cartopy.io.img_tiles as cimgt
import cartopy.crs as ccrs
import geopandas
import matplotlib.pyplot as plt


def rect_from_bound(xmin, xmax, ymin, ymax):
    """Returns list of (x,y)'s for a rectangle"""
    xs = [xmax, xmin, xmin, xmax, xmax]
    ys = [ymax, ymax, ymin, ymin, ymax]
    return [(x, y) for x, y in zip(xs, ys)]

# request data for use by geopandas
resolution          = '10m'
category            = 'cultural'
name                = 'admin_0_countries'
countries           = ['Norway', 'Sweden', 'Finland']
shpfilename         = shapereader.natural_earth(resolution, category, name)
df                  = geopandas.read_file(shpfilename)
extent              = [2, 32, 55, 72]
# get geometry of a country
for country in (countries):
    poly = [df.loc[df['ADMIN'] == country]['geometry'].values[0]]
    stamen_terrain              = cimgt.StamenTerrain()
    # projections that involved
    st_proj = stamen_terrain.crs  #projection used by Stamen images
    ll_proj = ccrs.PlateCarree()  #CRS for raw long/lat
    # create fig and axes using intended projection
    fig = plt.figure(figsize=(8,9))
    ax = fig.add_subplot(122, projection=st_proj)
    ax.add_geometries(poly, crs=ll_proj, facecolor='none', edgecolor='black')
    pad1 = 0.5  #padding, degrees unit
    exts = [poly[0].bounds[0] - pad1, poly[0].bounds[2] + pad1, poly[0].bounds[1] - pad1, poly[0].bounds[3] + pad1];
    ax.set_extent(exts, crs=ll_proj)
    # make a mask polygon by polygon's difference operation
    # base polygon is a rectangle, another polygon is simplified switzerland
    msk = Polygon(rect_from_bound(*exts)).difference( poly[0].simplify(0.01) )
    msk_stm  = st_proj.project_geometry (msk, ll_proj)  # project geometry to the projection used by stamen
    # get and plot Stamen images
    ax.add_image(stamen_terrain, 8) # this requests image, and plot
    # plot the mask using semi-transparency (alpha=0.65) on the masked-out portion
    ax.add_geometries( msk_stm, st_proj, zorder=12, facecolor='white', edgecolor='none', alpha=0.65)
    ax.gridlines(draw_labels=True)
    plt.show()

我有的是分开的地图。 THoguh 我只需要他们的一张地图。你能帮忙吗？谢谢。

Answer 1

你在工作中改编的代码对一个国家来说是好的。如果多个相邻的国家是新目标，则需要 select 所有这些国家并分解成一个单一的几何图形。只需要修改几行代码。

示例：新目标国家/地区：['Norway'、'Sweden'、'Finland']

需要替换的代码行：

poly = [df.loc[df['ADMIN'] == 'Switzerland']['geometry'].values[0]]

将其替换为以下代码行：

scan3 = df[ df['ADMIN'].isin(['Norway','Sweden', 'Finland']) ]
scan3_dissolved = scan3.dissolve(by='LEVEL')
poly = [scan3_dissolved['geometry'].values[0]]

你应该得到类似这样的情节：

使用 Cartopy 合并国家

Merge countries using Cartopy

matplotlib

cartopy

geopandas