Add high-resolution bottom topography to Cartopy map

Question

I'm migrating from basemap to Cartopy and would like to plot ocean bottom topography with high resolution for a limited area. In basemap I was using ETOPO1_Ice_g_gmt4.grd and transforming it to map coordinates based on documentation I had found somewhere. I don't know how to do that for Cartopy. Can anyone help? Cheers, Sünnje

Update: Code in Basemap

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)

Answer 1

First, here is a demo code and its output map using cartopy. It uses the projection and map's extent that you specified.

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()

To modify the code to plot your data, just change (xs,ys,zs) to (ttxi, ttyi, topodat).