0

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)
Snj
  • 33
  • 1
  • 9
  • You should share your code that does it in Basemap, and ask for modifications to migrate to Cartopy. – swatchai May 22 '21 at 06:42

1 Answers1

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

contour1

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

swatchai
  • 17,400
  • 3
  • 39
  • 58
  • Thank you very much, this is helpful.. What I still don't understand is how to transform the topographic data such that they get a shape that can be plotted with contourf. In Basemap I used the map.transform_scalar argument to generate that shape of the topodat. I'm unsure how I can do that now. Would you have an idea? – Snj May 27 '21 at 16:25
  • @Snj You better post that as a new question, more people will see and help. And the correct way of thanking is marking the `accept` button or voting up (when you have enough reputation). – swatchai May 29 '21 at 00:48
  • @Snj You need to get 2 more reputation points to enable voting up/down questions and answers, and you will be able to comment in other people's Q/A also. FYI, accepting an answer to your question worths 2 points. – swatchai Jun 03 '21 at 20:27