Software‎ > ‎Tutorials‎ > ‎

Visualizing Air Quality

Horizontal (Map) Visualization

Below is sample code that will produce a "tile-plot" of Odd-Oxygen with country and state maps where available. 

  1. Download the sample input data from the resources section of the site (sample_surface_geos_chem_output.nc) and save it to your scratch folder on HPC
  2. Copy the text below into a script file (e.g., plotting.py) and save it in the same folder with the sample output.
  3. Activate the aqc python installation (source /hpc/scratch/barronh/aqc/bin/activate)
  4. Run `python plotting.py`
  5. Download the test.png output from HPC via FileZilla
  6. View it.

Code: Sample code to visualize surface level Odd-Oxygen from a GEOS-Chem simulation.

#!/usr/bin/env python -i
from matplotlib import use
use('Agg')

from netCDF4 import Dataset;
import numpy as np;
from pylab import *
from mpl_toolkits.basemap import Basemap
import textwrap

# read netcdf file and create file object (f)
f = Dataset('sample_surface_geos_chem_output.nc')

# get the data variable from the file's variables
C = f.variables['IJ-AVG-$_Ox'];

# Disable tex for faster processing
# Enable tex for prettier fonts and latex support
rcParams['text.usetex'] = False

# Normalization is set for 10 deciles
# comment the next line to use standard matplotlib
# normalization or set your own normalization
norm = None

def dates_from(time):
    """
    Returns date from time variable
     time - variable with units property in the form "unit since YYYY-MM-DD HH:MM:SS UTC"
    """
    unit, sep, sincetxt = time.units.partition(' ')
    dummy, sep, date = sincetxt.partition(' ')
    date = date[:-4]
    refdate = datetime.datetime.strptime(date, '%Y-%m-%d %H:%M:%S')
    f = vectorize(lambda tv: refdate + datetime.timedelta(**{unit: float(tv)}))
    return f(time)

# Get coordinate variables (bounds are start and stop of intervals)
time_edge = f.variables['time_bounds'];
lat_edge = f.variables['latitude_bounds'][:];
lon_edge = f.variables['longitude_bounds'][:];

# Change lat/lon from (point, vertices) to (points,) shape
lat_edge = np.append(lat_edge[:, 0], lat_edge[-1, -1])
lon_edge = np.append(lon_edge[:, 0], lon_edge[-1, -1])

# create matrix from vectors
LON, LAT = np.meshgrid(lon_edge, lat_edge);

# Open a new figure for plotting
fig = figure(figsize = (8,5.5))

# Create a map object using the edges in the file
map = Basemap(llcrnrlon = lon_edge.min(), llcrnrlat = lat_edge.min(), urcrnrlon = lon_edge.max(), urcrnrlat = lat_edge.max());

if norm is None:
    nbins = 10
    boundaries = np.percentile(C[:].ravel(), np.linspace(0, 100, nbins + 1).tolist())
    norm = matplotlib.colors.BoundaryNorm(boundaries, 256)

ax = fig.add_axes([0.05, 0.175, 0.9, 0.675])
# Plot the data variable using LAT and LON matrices as corners
map.pcolor(LON, LAT, C[:], norm = norm)

# Add coastlines, countries, and states
map.drawcoastlines(); map.drawcountries(); map.drawstates();

# Add a super title (i.e., fig title)
suptitle('Odd-Oxygen from GEOS-Chem', size = 24)

# Convert time variable to date objects
start_date, end_date = dates_from(time_edge)[0]

# Add sub title (i.e., axis title) using date range
title('%s to %s' % (start_date.strftime('%Y-%m-%d %H:%M:%S'), end_date.strftime('%Y-%m-%d %H:%M:%S')), size = 20)
histstr = 'Processing: %s' % '\n'.join(textwrap.wrap(f.history.strip(), 140))
if rcParams['text.usetex']:
    histstr = histstr.replace('_', '\_')

fig.text(0.005, 0.005, histstr, horizontalalignment = 'left', verticalalignment = 'bottom', size = 8)
cax = fig.add_axes([0.125, 0.12, 0.75, 0.03])
cbar = colorbar(cax = cax, orientation = 'horizontal')
cax.yaxis.set_label_position("right")
xl = cax.set_ylabel(C.units.strip(), rotation = 0)

# Save the figure
savefig('test_surface.png')

Example output:

Zonal Visualization

Below is sample code that will produce a "zonal-plot" of Odd-Oxygen. 

  1. Download the sample input data from the resources section of the site (sample_zonal_geos_chem_output.nc) and save it to your scratch folder on HPC
  2. Copy the text below into a script file (e.g., zonal_plotting.py) and save it in the same folder with the sample output.
  3. Activate the aqc python installation (source /hpc/scratch/barronh/aqc/bin/activate)
  4. Run `python zonal_plotting.py`
  5. Download the test_zonal.png output from HPC via FileZilla
  6. View it.

Code: Sample code to visualize zonal Odd-Oxygen from a GEOS-Chem simulation.

#!/usr/bin/env python -i
from matplotlib import use
use('Agg')

from netCDF4 import Dataset;
import numpy as np;
from pylab import *
import textwrap

# read netcdf file and create file object (f)
f = Dataset('sample_zonal_geos_chem_output.nc')

# get the data variable from the file's variables
C = f.variables['IJ-AVG-$_Ox'];

# Disable tex for faster processing
# Enable tex for prettier fonts and latex support
rcParams['text.usetex'] = False

# Normalization is set for 10 deciles
# comment the next line to use standard matplotlib
# normalization or set your own normalization
norm = None


def dates_from(time):
    """
    Returns date from time variable
     time - variable with units property in the form "unit since YYYY-MM-DD HH:MM:SS UTC"
    """
    unit, sep, sincetxt = time.units.partition(' ')
    dummy, sep, date = sincetxt.partition(' ')
    date = date[:-4]
    refdate = datetime.datetime.strptime(date, '%Y-%m-%d %H:%M:%S')
    f = vectorize(lambda tv: refdate + datetime.timedelta(**{unit: float(tv)}))
    return f(time)

# Get coordinate variables (bounds are start and stop of intervals)
time_edge = f.variables['time_bounds'];
lat_edge = f.variables['latitude_bounds'][:];
lay_edge = f.variables['layer48'][:];

# Change lat from (point, vertices) to (points,) shape
lat_edge = np.append(lat_edge[:, 0], lat_edge[-1, -1])

# create matrix from vectors
LAT, LAY = np.meshgrid(lat_edge, lay_edge);

# Open a new figure for plotting
fig = figure(figsize = (8,7))
ax = fig.add_axes([0.1, 0.15, 0.75, 0.725])

# Create evenly distributed boundaries
if norm is None:
    nbins = 10
    boundaries = np.percentile(C[:].ravel(), np.linspace(0, 100, nbins + 1).tolist())
    norm = matplotlib.colors.BoundaryNorm(boundaries, 256)

# Plot C
pcolor(LAT, LAY, C[:].squeeze(), norm = norm)

# Create an axis for the colorbar
cax = fig.add_axes([0.875, 0.15, 0.025, 0.725])

# Add the colorbar
cbar = colorbar(ax = ax, cax = cax)

# Put units on the x axis
xl = cbar.ax.set_xlabel(C.units.strip(), rotation = 0, size = 14)

# Reverse y limits because the y-axis is in pressure units
ylim(LAY.max(), LAY.min())

# Fit axis to data
xlim(LAT.min(), LAT.max())


ylabel('pressure (hPa)')
xlabel('latitude (degrees N)')

# Add a super title (i.e., fig title)
suptitle('Odd-Oxygen from GEOS-Chem', size = 24)

# Convert time variable to date objects
start_date, end_date = dates_from(time_edge)[0]

# Add sub title (i.e., axis title) using date range
title('%s to %s' % (start_date.strftime('%Y-%m-%d %H:%M:%S'), end_date.strftime('%Y-%m-%d %H:%M:%S')), size = 20)

# Add history to figure
histstr = 'Processing: %s' % '\n'.join(textwrap.wrap(f.history.strip(), 140))
if rcParams['text.usetex']:
    histstr = histstr.replace('_', '\_')

fig.text(0.005, 0.005, histstr, horizontalalignment = 'left', verticalalignment = 'bottom', size = 8)

# Save the figure
savefig('test_zonal.png')
 

 

Example Output:


Comments