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

如何解决将高分辨率底部地形添加到 Cartopy 地图

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

更新：底图中的代码

map = Basemap(projection = 'merc',llcrnrlat = 67.2,urcrnrlat = 69.5,\
llcrnrlon = 8,urcrnrlon = 16.5,lat_ts = 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)

解决方法

首先，这是一个演示代码及其使用 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,zorder=12)

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

要修改代码以绘制数据，只需将 (xs,zs) 更改为 (ttxi,ttyi,topodat)。