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.

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
Copy the text below into a script file (e.g., plotting.py) and save it in the same folder with the sample output.
Activate the aqc python installation (source /hpc/scratch/barronh/aqc/bin/activate)
Run `python plotting.py`
Download the test.png output from HPC via FileZilla
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.

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
Copy the text below into a script file (e.g., zonal_plotting.py) and save it in the same folder with the sample output.
Activate the aqc python installation (source /hpc/scratch/barronh/aqc/bin/activate)
Run `python zonal_plotting.py`
Download the test_zonal.png output from HPC via FileZilla
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:

Page updated

Report abuse