将高分辨率底部地形添加到Cartopy地图

Question

我正在从basemap迁移到Cartopy，并希望以高分辨率绘制有限区域的海底地形。 在basemap中，我使用ETOPO1_Ice_g_gmt4.grd并根据我在某处找到的文档将其转换为地图坐标。我不知道如何为Cartopy做到这一点。有人能帮忙吗？ 干杯，苏奈

更新：Basemap中的代码

地图=基本地图（投影='merc'，llcrnrlat=67.2，urcrnrlat=69.5，\
llcrnrlon=8，urcrnrlon=16.5，lat=67.5，）

topoFile = nc.NetCDFFile('/home/sunnje/data/ETOPO1_Ice_g_gmt4.grd','r')                                                       
topoLons = topoFile.variables['x'][:]                                                                                         
topoLats = topoFile.variables['y'][:]                                                                                         
topoZ = topoFile.variables['z'][:]                                                                                            
                                                                                                                              
# transform to nx x ny regularly spaced 1km native projection grid                                                            
nx = int((map.xmax - map.xmin)/1000.)+1                                                                                       
ny = int((map.ymax - map.ymin)/1000.)+1                                                                                       
                                                                                                                              
topodat = map.transform_scalar(topoZ,topoLons,topoLats,nx,ny)                                                                 
                                                                                                                              
tyi = np.linspace(map.ymin,map.ymax,topodat.shape[0])                                                                         
txi = np.linspace(map.xmin,map.xmax,topodat.shape[1])                                                                         
ttxi, ttyi = np.meshgrid(txi,tyi)                                                                                             
                                                                                                                              
cm = map.contour(ttxi, ttyi, topodat)

Answer 1

首先，这里是一个演示代码及其使用cartopy的输出映射。它使用您指定的投影和贴图范围

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from scipy.stats import multivariate_normal

# set extent necessary for data and plots
map_ext = [8, 16.5, 67.2, 69.5] #longmin,longmax,latmin,latmax

# for demo purposes
# set needed values for data creation
xmean, ymean = (map_ext[0]+map_ext[1])/2.0, (map_ext[2]+map_ext[3])/2.0
mu = np.array([xmean, ymean])
covar = np.array([[ 10. , -0.5], [-0.5,  10.5]])
lon_range = np.linspace(-180, 180, 200)
lat_range = np.linspace(-90, 90, 100)
xs,ys = np.meshgrid(lon_range, lat_range)
pos = np.empty(xs.shape + (2,))
pos[:, :, 0] = xs
pos[:, :, 1] = ys
# generate values as a function of (x,y) for contour genereation
zs = multivariate_normal(mu, covar).pdf(pos)

# setup projection for the map
projection = ccrs.Mercator(latitude_true_scale=67.5)
# create figure, axis for the map
fig, ax = plt.subplots(figsize=(8,6), subplot_kw={'projection': projection})

ax.set_extent(map_ext)
ax.coastlines()

ax.contourf(xs, ys, zs, transform=ccrs.PlateCarree(), zorder=10, alpha=0.6)
ax.contour(xs, ys, zs, transform=ccrs.PlateCarree(), zorder=12)

ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
plt.show()

要修改代码以绘制数据，只需将（xs、ys、zs）更改为（ttxi、ttyi、topodat）