Create maps using Etopo1 and matplotlib

I often create maps for visualizing scientific data for a specific geographic region, or to simply show the bathymetry of an area. I find that the Etopo1 Global Relief Model in combination with Python, matplotlib, and Basemap provides the essential tools and data necessary to create high-resolution, nice looking maps, suitable for publication. You can download the from the NOAA website and according to NOAA: “ETOPO1 is a 1 arc-minute global relief model of Earth’s surface that integrates land topography and ocean bathymetry. It was built from numerous global and regional data sets, and is available in “Ice Surface” (top of Antarctic and Greenland ice sheets) and “Bedrock” (base of the ice sheets) versions.” So, you actually got both bathymetry, topography, and ice in one file.

As an example on how to use these data I created a script that creates a map for a pre-defined region. The example creates a map for the Bay of Biscay (latitude:40-81N, longitude:-10-75E). Since the ETOPO1 file is quite large (445MB) I added a function to the script that cuts out only the area of interest from the ETOPO1 file. This saves space and time (see function: findSubsetIndices). Note that this script require a number of Python modules including the Python NetCDF4 interface, Basemap, Numpy, and Matplotlib. Details on how to install all these modules are found here.  I also smooth the ETOPO1 data with a Laplace filter that you can download below. The Laplace filter was written by Bjørn Aadlandsvik at IMR in Bergen. The Laplace filter is called from the main python script (createMapsEtopo1.py). Finally, the mpl_util.py function helps generate smooth color contours automatically.
download
createMapsEtopo1.py (1660)

import os, sys, datetime, string
import numpy as np
from netCDF4 import Dataset
import numpy.ma as ma
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, shiftgrid, NetCDFFile
from pylab import *
import laplaceFilter
import mpl_util

__author__   = 'Trond Kristiansen'
__email__    = 'trond.kristiansen (at) imr.no'
__created__  = datetime.datetime(2008, 8, 15)
__modified__ = datetime.datetime(2009, 7, 21)
__version__  = "1.0"
__status__   = "Development"

def findSubsetIndices(min_lat,max_lat,min_lon,max_lon,lats,lons):
"""Array to store the results returned from the function"""
res=np.zeros((4),dtype=np.float64)
minLon=min_lon; maxLon=max_lon
distances1 = []; distances2 = []
indices=[]; index=1
for point in lats:
s1 = max_lat-point # (vector subtract)
s2 = min_lat-point # (vector subtract)
distances1.append((np.dot(s1, s1), point, index))
distances2.append((np.dot(s2, s2), point, index-1))
index=index+1
distances1.sort()
distances2.sort()
indices.append(distances1[0])
indices.append(distances2[0])
distances1 = []; distances2 = []; index=1

for point in lons:
s1 = maxLon-point # (vector subtract)
s2 = minLon-point # (vector subtract)
distances1.append((np.dot(s1, s1), point, index))
distances2.append((np.dot(s2, s2), point, index-1))
index=index+1
distances1.sort()
distances2.sort()
indices.append(distances1[0])
indices.append(distances2[0])

""" Save final product: max_lat_indices,min_lat_indices,max_lon_indices,min_lon_indices"""
minJ=indices[1][2]
maxJ=indices[0][2]
minI=indices[3][2]
maxI=indices[2][2]
res[0]=minI; res[1]=maxI; res[2]=minJ; res[3]=maxJ;
return res

def makeMap(lonStart,lonEnd,latStart,latEnd,name,stLon,stLat):
plt.figure(figsize=(8,8))
"""Get the etopo2 data"""
etopo1name='ETOPO1_Ice_g_gmt4.grd'
etopo1 = Dataset(etopo1name,'r')
lons = etopo1.variables["lon"][:]
lats = etopo1.variables["lat"][:]

res = findSubsetIndices(latStart-5,latEnd+5,lonStart-40,lonEnd+10,lats,lons)
lon,lat=np.meshgrid(lons[res[0]:res[1]],lats[res[2]:res[3]])
print "Extracted data for area %s : (%s,%s) to (%s,%s)"%(name,lon.min(),lat.min(),lon.max(),lat.max())
bathy = etopo1.variables["z"][int(res[2]):int(res[3]),int(res[0]):int(res[1])]
bathySmoothed = laplaceFilter.laplace_filter(bathy,M=None)
levels=[-6000,-5000,-3000, -2000, -1500, -1000,-500, -400, -300, -250, -200, -150, -100, -75, -65, -50, -35, -25, -15, -10, -5, 0]
if lonStart< 0 and lonEnd < 0:
lon_0= - (abs(lonEnd)+abs(lonStart))/2.0
else:
lon_0=(abs(lonEnd)+abs(lonStart))/2.0
map = Basemap(llcrnrlat=latStart,urcrnrlat=latEnd,\
llcrnrlon=lonStart,urcrnrlon=lonEnd,\
rsphere=(6378137.00,6356752.3142),\
resolution='l',area_thresh=1000.,projection='lcc',\
lat_1=latStart,lon_0=lon_0)
x, y = map(lon,lat)
map.drawcoastlines()
map.drawcountries()
map.fillcontinents(color='grey')
map.drawmeridians(np.arange(lons.min(),lons.max(),10),labels=[0,0,0,1])
map.drawparallels(np.arange(lats.min(),lats.max(),4),labels=[1,0,0,0])
#map.bluemarble()

CS1 = map.contourf(x,y,bathySmoothed,levels,
cmap=mpl_util.LevelColormap(levels,cmap=cm.Blues_r),
extend='upper',
alpha=1.0,
origin='lower')
CS1.axis='tight'
"""Plot the station as a position dot on the map"""
pt,ypt = map(stLon,stLat)
map.plot([xpt],[ypt],'ro', markersize=10)
plt.text(xpt+100000,ypt+100000,name)
plt.title('Area %s'%(name))
plotfile='figures/map_'+str(name)+'.pdf'
plt.savefig(plotfile,dpi=150,orientation='portrait')
plt.show()

def main():
names=['bayOfBiscay']
lat_start=[40]
lat_end  =[81]
lon_start=[-10]
lon_end  =[75]
"""List of stations for each area"""
stationlonlist=[ 2.4301, -22.6001, -47.0801,  13.3801, -67.2001]
stationlatlist=[54.5601, 63.7010,  60.4201,  67.5001,  41.6423]
for i in range(len(lat_start)):
if i==0:
makeMap(lon_start[i],lon_end[i],lat_start[i],lat_end[i],names[i],stationlonlist[i],stationlatlist[i])
if __name__ == "__main__":
main()

download
Laplace Filter (1125)

import numpy as np
def laplace_X(F,M):
"""1D Laplace Filter in X-direction (axis=1)"""
jmax, imax = F.shape
# Add strips of land
F2 = np.zeros((jmax, imax+2), dtype=F.dtype)
F2[:, 1:-1] = F
M2 = np.zeros((jmax, imax+2), dtype=M.dtype)
M2[:, 1:-1] = M
MS = M2[:, 2:] + M2[:, :-2]
FS = F2[:, 2:]*M2[:, 2:] + F2[:, :-2]*M2[:, :-2]
return np.where(M > 0.5, (1-0.25*MS)*F + 0.25*FS, F)

def laplace_Y(F,M):
"""1D Laplace Filter in Y-direction (axis=1)"""
jmax, imax = F.shape
# Add strips of land
F2 = np.zeros((jmax+2, imax), dtype=F.dtype)
F2[1:-1, :] = F
M2 = np.zeros((jmax+2, imax), dtype=M.dtype)
M2[1:-1, :] = M
MS = M2[2:, :] + M2[:-2, :]
FS = F2[2:, :]*M2[2:, :] + F2[:-2, :]*M2[:-2, :]
return np.where(M > 0.5, (1-0.25*MS)*F + 0.25*FS, F)

def laplace_filter(F, M=None):
if M == None:
M = np.ones_like(F)

return 0.5*(laplace_X(laplace_Y(F, M), M) +
laplace_Y(laplace_X(F, M), M))

download
mpl_util.py (1053)

import matplotlib
from numpy import ma
import pylab as pl
"""A set of utility functions for matplotlib"""

# ------------------
# Plot land mask
# ------------------

def landmask(M, color='0.8'):

   # Make a constant colormap, default = grey
   constmap = pl.matplotlib.colors.ListedColormap([color])

   jmax, imax = M.shape
   # X and Y give the grid cell boundaries,
   # one more than number of grid cells + 1
   # half integers (grid cell centers are integers)
   X = -0.5 + pl.arange(imax+1)
   Y = -0.5 + pl.arange(jmax+1)

   # Draw the mask by pcolor
   M = ma.masked_where(M > 0, M)
   pl.pcolor(X, Y, M, shading='flat', cmap=constmap)

# -------------
# Colormap
# -------------

def LevelColormap(levels, cmap=None):
    """Make a colormap based on an increasing sequence of levels"""
   
    # Start with an existing colormap
    if cmap == None:
        cmap = pl.get_cmap()

    # Spread the colours maximally
    nlev = len(levels)
    S = pl.arange(nlev, dtype='float')/(nlev-1)
    A = cmap(S)

    # Normalize the levels to interval [0,1]
    levels = pl.array(levels, dtype='float')
    L = (levels-levels[0])/(levels[-1]-levels[0])

    # Make the colour dictionary
    R = [(L[i], A[i,0], A[i,0]) for i in xrange(nlev)]
    G = [(L[i], A[i,1], A[i,1]) for i in xrange(nlev)]
    B = [(L[i], A[i,2], A[i,2]) for i in xrange(nlev)]
    cdict = dict(red=tuple(R),green=tuple(G),blue=tuple(B))

    # Use
    return matplotlib.colors.LinearSegmentedColormap(
        '%s_levels' % cmap.name, cdict, 256)

4 thoughts on “Create maps using Etopo1 and matplotlib”

  1. Thanks !!!

    The LevelColormap was superuseful !

    I could not find anywhere but here how to change colormap in Basemap
    (of course except not very tempting cmap=cm.s3pcpn and cm.GMT_haxby from oficial matplotlib basemap) and you have a very neat and flexible solution.

Leave a Reply

Your email address will not be published. Required fields are marked *