合并使用Cartopy的国家

问题描述

我正在使用以下代码为一个地区同时绘制瑞典，挪威和芬兰的地图.但是，我正在为此而苦苦挣扎.我正在关注这个示例，即Matplotlib Cartopy Color One国家中的Python映射.

  from shapely.geometry导入多边形从cartopy.io导入shapereader导入cartopy.io.img_tiles作为cimgt将cartopy.crs导入为ccrs进口大熊猫导入matplotlib.pyplot作为pltdef rect_from_bound(xmin，xmax，ymin，ymax):"返回矩形的(x，y)"列表.xs = [xmax，xmin，xmin，xmax，xmax]ys = [ymax，ymax，ymin，ymin，ymax]返回[(x，y)表示zip(xs，ys)中的x，y]#请求供Geopandas使用的数据分辨率='10m'类别=文化"名称='admin_0_countries'国家= ['挪威'，'瑞典'，'芬兰']shpfilename = shapereader.natural_earth(分辨率，类别，名称)df = geopandas.read_file(shpfilename)范围= [2、32、55、72]#获取一个国家的几何针对(国家/地区)的国家/地区:poly = [df.loc [df ['ADMIN'] == country] ['geometry'].values [0]]stamen_terrain = cimgt.StamenTerrain()#涉及的预测st_proj = stamen_terrain.crs#雄蕊图像使用的投影ll_proj = ccrs.PlateCarree()#CRS用于原始多头/空头#使用预期的投影创建无花果和轴无花果= 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，度单位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)#通过多边形的差值运算制作遮罩多边形#基本多边形是一个矩形，另一个多边形是简化的瑞士msk =多边形(rect_from_bound(* exts)).difference(poly [0] .simplify(0.01))msk_stm = st_proj.project_geometry(msk，ll_proj)#雄蕊所用投影的项目几何#获取并绘制雄蕊图像ax.add_image(stamen_terrain，8)#这将请求图像并绘制#在遮罩部分上使用半透明(alpha = 0.65)绘制遮罩ax.add_geometries(msk_stm，st_proj，zorder = 12，facecolor ='white'，edgecolor ='none'，alpha = 0.65)ax.gridlines(draw_labels = True)plt.show() 

我所拥有的是分开的地图.我只需要一张地图.你能帮忙吗?谢谢.

解决方案

代码

I am using the following code to make a map for Sweden, Norway and Finland together as one area. however, I am struggling with it. I'm following this example, Python Mapping in 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()

What I have is separated maps. THoguh I need only one map of them.Can you please help?Thank you.

解决方案

The code here that you adapted to your work is good for a single country. If multiple contiguous countries are new target, one need to select all of them and dissolve into a single geometry. Only a few lines of code need to be modified.

Example: new target countries: ['Norway','Sweden', 'Finland']

The line of code that need to be replaced:

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

Replace it with these lines of code:

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

And you should get a plot similar to this:

这篇关于合并使用Cartopy的国家的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！