绘制正方形 Cartopy 地图

Question

我需要使用 Cartopy 绘制方形地图。我目前为我的地图使用以下代码：

plt.figure(figsize = (15, 15))

img = cimgt.GoogleTiles()

ax = plt.axes(projection = img.crs)
ax.set_extent((d['longitude'].min() - 0.05, d['longitude'].max() + 0.05,
               d['latitude'].min() - 0.05, d['latitude'].max() + 0.05))

ax.add_image(img, 10, interpolation = 'bicubic')

plt.scatter(d['longitude'], d['latitude'], transform = ccrs.PlateCarree(),
            c = '#E8175D', s = 14)

这很好用，除了地图不是正方形。相反，它只是适合情节的 (15, 15) 方块。

我想在左侧和右侧添加更多的地图，使绘图完美正方形而不扭曲。简单地将范围设置为相同的纬度和经度差异并不能完成这项工作，因为纬度和经度覆盖 Google（和大多数其他）地图投影中的不同距离。我也找到了 this post，但据我所知，这里的目的是扭曲地图。

我希望有人知道如何做到这一点。好像Cartopy在这方面不是很直观

Answer 1

要获得方形范围，您需要在地图投影坐标中指定它。这涉及到一些坐标变换。这是您需要的代码片段。

# crs of your choice
crg = cimgt.StamenTerrain().crs    # or cimgt.GoogleTiles().crs

# set map limits, in degrees
lonmin, lonmax = -22, -15
latmin, latmax = 63, 65

# do coordinate transformation
LL = crg.transform_point(lonmin, latmin, ccrs.Geodetic())
UR = crg.transform_point(lonmax, latmax, ccrs.Geodetic())
EW = UR[0] - LL[0]
SN = UR[1] - LL[1]

# get side of the square extent (in map units, usually meters)
side = max(EW, SN)    # larger value is in effect
mid_x, mid_y = LL[0]+EW/2.0, LL[1]+SN/2.0  # center location

# the extent preserves the center location
extent = [mid_x-side/2.0, mid_x+side/2.0, mid_y-side/2.0, mid_y+side/2.0]

# this sets square extent
# crs=crg signifies that projection coordinates is used in extent
ax.set_extent(extent, crs=crg)

希望对您有所帮助。

编辑

这是一个完整的工作代码及其生成的映射。

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.io.img_tiles as cimgt

def make_map(projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(10, 10),
                           subplot_kw=dict(projection=projection))
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = gl.ylabels_right = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax

request = cimgt.StamenTerrain()   # very responsive
crg = request.crs   #crs of the projection
fig, ax = make_map(projection = crg)

# specify map extent here
lonmin, lonmax = -22, -15
latmin, latmax = 63, 65

LL = crg.transform_point(lonmin, latmin, ccrs.Geodetic())
UR = crg.transform_point(lonmax, latmax, ccrs.Geodetic())
EW = UR[0] - LL[0]
SN = UR[1] - LL[1]
side = max(EW,SN)
mid_x, mid_y = LL[0]+EW/2.0, LL[1]+SN/2.0  #center location

extent = [mid_x-side/2.0, mid_x+side/2.0, mid_y-side/2.0, mid_y+side/2.0]   # map coordinates, meters

ax.set_extent(extent, crs=crg)
ax.add_image(request, 8)

# add a marker at center of the map
plt.plot(mid_x, mid_y, marker='o', \
         color='red', markersize=10, \
         alpha=0.7, transform = crg)

plt.show()

绘制正方形 Cartopy 地图

Plot square Cartopy map

python

matplotlib

cartopy