PODAAC

Using Python Packages Matplotlib and Basemap to Create Plots Using Your Data

The following two packages are incredibly versatile in generating various scientific plots of geo-gridded Earth science data: Matplotlib and the extension for Basemap for Python. In this example, we will will create a time series plot using netCDF files and also how to create some simple contours.

Note: we will need the netCDF4 library so that we can easily access and manipulate netCDF data. We can download the netCDF4 library here:
http://netcdf4-python.googlecode.com/sv ... odule.html

Here is a simple Python program that reads in data from a netCDF file and uses it to create a contour.

###################################################################
# import necessary libraries and rename them so they are easier to use
import numpy as N
import netCDF4 as S
import matplotlib.pyplot as plt
import mpl_toolkits.basemap as bm
from netCDF4 import Dataset

# create a Dataset object called Dset using the Dataset constructor from the
# netCDF4 library
# To use this program simply replace <netCDF_File> with the name of your file.
Dset = Dataset("<netCDF_File>.nc", "r", format ="NETCDF4")

# Get information from our Dataset object, Dset and store it in varialbes
# my netCDF file contains information on global sea surface temperature.
# Dset.variables calls my Dataset object and goes to the variables section of
# the information that it holds. We find the names of the variables that we need
# and assign them to variables in python. The array that stores all of the SST
# information is called “analysed_sst” and our first line gets that data and
# stores it in the variable sst.
sst = Dset.variables["analysed_sst"][:]
sst_units = Dset.variables["analysed_sst"].units

# Here we do the same thing except now we are getting information about the lon and lat variables.
lon_scale = Dset.variables["lon"][:]
lon_units = Dset.variables["lon"].units
lon_max = Dset.variables["lon"].valid_max
lon_min = Dset.variables["lon"].valid_min

lat_scale = Dset.variables["lat"][:]
lat_units = Dset.variables["lat"].units
lat_max = Dset.variables["lat"].valid_max
lat_min = Dset.variables["lat"].valid_min

#____________________
# These lines create the basemap for the data.
Mbase = bm.Basemap(projection = 'cyl', llcrnrlat = -90.0, llcrnrlon = -180.0, urcrnrlat = 90.0, urcrnrlon = 180.0)

Mbase.drawparallels(N.array([-20, -10, 0, 10 ,20]), labels=[0,0,0,1])
Mbase.drawmeridians(N.array([-180, -90, 0, 90, 180]), labels = [0,0,0,1])

Mbase.drawcoastlines()
Mbase.fillcontinents(color = 'blue') # blue can be changed to any color that you would like
#____________________

mymap = plt.contour(lon_scale, lat_scale, sst, 15, colors = 'k')
mymap2 = plt.contourf(lon_scale, lat_scale, sst, 15, cmap=plt.cm.Greens)
plt.colorbar(mymap2, orientation = 'horizontal')

plt.axis([lon_min,lon_max,lat_min,lat_max])
plt.xlabel(lon_units)
plt.ylabel(lat_units)
plt.show()


###################################################################

If you have some sort of two dimensional data you may want to create a simple time series. Here is an example on how to create time series.

###################################################################
import numpy as N

import matplotlib.pyplot as plt
from netCDF4 import Dataset

Dset = Dataset("netCDF_File.nc", "r", format ="NETCDF4")

sst = Dset.variables["analysed_sst"][:]
sst_units = Dset.variables["analysed_sst"].units
time_units = Dset.variables["time"].units


plt.plot(sst)
plt.xlabel("Units of 12HRS into January, 2008")
plt.ylabel(sst_units)
plt.show()
###################################################################