plot gebco data in python basemap

Question

I have downloaded some gebco bathymetry data as a netCDF file. I would like to plot it with python-basemap. I have tried,

import netCDF4
from mpl_toolkits.basemap import Basemap


# Load data
dataset = netCDF4.Dataset('/home/david/Desktop/GEBCO/gebco_08_-30_45_5_65.nc')

# Extract variables
x = dataset.variables['x_range']
y = dataset.variables['y_range']
spacing = dataset.variables['spacing']

# Data limits
nx = (x[-1]-x[0])/spacing[0]   # num pts in x-dir
ny = (y[-1]-y[0])/spacing[1]   # num pts in y-dir

# Reshape data
zz = dataset.variables['z']
Z = zz[:].reshape(ny, nx)



# setup basemap.
m = Basemap(llcrnrlon=-30,llcrnrlat=45.0,urcrnrlon=5.0,urcrnrlat=65.0,
            resolution='i',projection='stere',lon_0=-15.0,lat_0=55.0)


# Set up grid
lons, lats = m.makegrid(nx, ny)
x, y = m(lons, lats)

m.contourf(x, y, flipud(Z))
m.fillcontinents(color='grey')
m.drawparallels(np.arange(10,70,10), labels=[1,0,0,0])
m.drawmeridians(np.arange(-80, 5, 10), labels=[0,0,0,1])

this gives the figure below, clearly not correct. The problem stems from how the areas are defined. For basemap area is defined by lower left corner lat,lon and upper right corner lat, lon. But the gebco data takes a maximum and minimum lon/lat defined along a center line. Anyone have any experience with gebco data or see a solution??

thanks D map

What's the motivation for creating X and Y (btw you do it twice) or for Z=zz[:].reshape(shape(X))? You also haven't defined nx or ny in the code provided. What do they represent? — Spencer Hill
– Spencer Hill, Commented Feb 19, 2014 at 19:43
Sorry, that was a bit unclear. Code has been updated. zz is reshaped as it is a single column vector. — Dave
– Dave, Commented Feb 19, 2014 at 20:25
What are x and y? If they are longitude and latitude arrays respectively, there is no need to use makegrid with nx and ny; just do x,y = m(*np.meshgrid(x, y)). Perhaps that also fixes the problem? — Spencer Hill
– Spencer Hill, Commented Feb 19, 2014 at 22:35
x and y contain the max and min values for lon and lat, spacing is the grid step size. I create the lon lat arrays from them. Using the arrays as you suggest instead of makegrid works perfectly, thanks. — Dave
– Dave, Commented Feb 20, 2014 at 14:28

Rich Signell · Accepted Answer · 2014-03-05 17:52:47Z

so just for the record, here's the answer that works, using the comments above:

import netCDF4
from mpl_toolkits.basemap import Basemap

# Load data
dataset = netCDF4.Dataset('/usgs/data1/rsignell/bathy/gebco_08_-30_-45_5_65.nc')

# Extract variables
x = dataset.variables['x_range']
y = dataset.variables['y_range']
spacing = dataset.variables['spacing']

# Compute Lat/Lon
nx = (x[-1]-x[0])/spacing[0]   # num pts in x-dir
ny = (y[-1]-y[0])/spacing[1]   # num pts in y-dir

lon = np.linspace(x[0],x[-1],nx)
lat = np.linspace(y[0],y[-1],ny)

# Reshape data
zz = dataset.variables['z']
Z = zz[:].reshape(ny, nx)

# setup basemap.
m = Basemap(llcrnrlon=-30,llcrnrlat=45.0,urcrnrlon=5.0,urcrnrlat=65.0,
            resolution='i',projection='stere',lon_0=-15.0,lat_0=55.0)

x,y = m(*np.meshgrid(lon,lat))

m.contourf(x, y, flipud(Z));
m.fillcontinents(color='grey');
m.drawparallels(np.arange(10,70,10), labels=[1,0,0,0]);
m.drawmeridians(np.arange(-80, 5, 10), labels=[0,0,0,1]);

which produces this plot. enter image description here

Collectives™ on Stack Overflow

plot gebco data in python basemap

1 Answer 1

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related