This code makes plots from the NASA Goddard Space Flight Center two-dimensional atmospheric chemistry and dynamics model, specifically the version by Brian Thomas. Three primary runs are used:
If you would like to run this yourself, you will need the specout*.nc files. You will also need to change the directory locations within this code.
Cell 1: import statements Cell 2: definitions of chemicals, latitude, and altitude arrays Cell 3: defining some functions used in plotting Cell 4: main plotting cell. Several different plots are defined, and these can all be turned on or off with the comments at the bottom of the cell. Cells 5-7: the exact code used to create the figures. Most of this is copy/pasted from the plotting routines in Cell 4, but placed here for easy replication.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
from matplotlib import ticker, cm
from matplotlib import colors
import matplotlib.patheffects as pe
import subprocess
from scipy.signal import argrelextrema
from scipy import interpolate
from scipy.optimize import curve_fit
from scipy.io import netcdf
# import netCDF4 as nc
%matplotlib inline
# Plot GSFC data from the specout files
# Similar to ncl or Panoply, but now I have control over it
# Why? Cause I want to plot column densities of anything, not just what's in the coldens file
CONTROLRUN = 'controlRun_PMCs_norestart'
# all molecule units are mol/cm^3
# Age has unknown units
# OzoneDiag is in DU/km (DU is Dobson units)
# Temp is temperature in K
# UbarD is mean zonal wind speed in m/s
# HeatOut is heating rate in K/day
# CoolOut is cooling rate in K/day
# W is vertical velocity in cm/s
# V is horizontal velocity (north-south) in cm/s
cname = [
'O(3p)', 'O(1d)', 'O2', 'O3', 'NO', 'NO2', 'NO3', 'N2O5', 'N', 'HNO3', 'N2O', 'H', 'OH', 'HO2',\
'H2O', 'H2O2', 'H2', 'CH4', 'CO', 'CO2', 'ClOx', 'CH3O2', 'CH2O', 'CH3OOH', 'HOCl',\
'ClO3', 'Cl', 'ClO', 'HCl', 'ClONO2', 'NOy', 'NOx', 'Cly', 'CFCl3', 'CF2Cl2', 'CCl4',\
'CH3Cl', 'HO2NO2', 'Ox', 'CH3CCl3', 'O2(1Delta)', 'HONO', 'Br2', 'BrO', 'Br', 'HBr',\
'BrONO2', 'Bry', 'CH3Br', 'CBrF3', 'CBrClF2', 'CHClF2', 'C2Cl3F3', 'C2Cl2F4', 'C2ClF5',\
'HF', 'CF2O', 'BrOx', 'H2O_rainout', 'BrCl', 'Cl2O2', 'OClO', 'ClNO2', 'Cl2', 'ClOO',\
'HNO3_solid', 'LiqH2O', 'HOBr', 'CH3CCl2F', 'CH3CClF2', 'CHCl2CF3', 'C2Br2F4', 'CHx',\
'HOx', 'CBr2F2', 'CH2Br2', 'CHBr3', 'Age', 'BrONO', 'OzoneDiag',\
'Temp', 'UbarD', 'HeatOut', 'CoolOut', 'W', 'V'
]
# define latitude and altitude arrays
# note that zalt90arr is approximate height
# for pressure-based height, need to calculate it with calc_z_specoutnc() defined later
lat4arr = np.linspace(-88, 88, 45, endpoint=True)
zalt90arr = np.concatenate((np.linspace(0.5, 59.5, 60, endpoint=True), np.linspace(61, 91, 16, endpoint=True)))
# Rather than load the entire freaking 20 yr dataset into RAM,
# (too much data to load all 24+ years, so can't include all)
# make a function that reads file data when needed for plots.
# This will be slightly more time-consuming when making plots,
# but will be able to do all years without maxing out the RAM.
def check_variables(species, day, alt, lat):
# Do some basic checks to make sure input data is okay
# If variables are outside range, return a string
# 'species' is a string corresponding to 'cname' above
if species in cname:
si = cname.index(species)
else:
errorstring = '#####\n##### Error: species not found in cname\n#####'
return errorstring
# 'day' is a list of ints from 0 to (720 * number of runs - 1)
if isinstance(day, int) or isinstance(day, np.int64):
day = np.array([day])
if isinstance(day, list):
day = np.array(day)
if not all(t >= 0 for t in day):
errorstring = '#####\n##### Error: values for day less than 0.\n#####'
return errorstring
# 'alt' is a list of ints from 0 to 75
if isinstance(alt, int):
alt = np.array([alt])
if isinstance(alt, list):
alt = np.array(alt)
if not all(0 <= a <= 75 for a in alt):
errorstring = '#####\n##### Error: values for alt outside 0 - 75.\n#####'
return errorstring
# 'lat' is a list of ints from 0 to 44
if isinstance(lat, int):
lat = np.array([lat])
if isinstance(lat, list):
lat = np.array(lat)
if not all(0 <= l <= 44 for l in lat):
errorstring = '#####\n##### Error: values for lat outside 0 - 44.\n#####'
return errorstring
return si, day, alt, lat
def get_filenames(runname, day, headdir):
# determine which files are needed based on runname and days
dirname = headdir+runname+'/'
run = day // 720 + 1
modr = day % 720
# split modr into sequence of arrays (super janky but it works)
# e.g., run = [1, 1, 2, 2, 2] and modr = [718, 719, 0, 1, 2]
# yields modrsplit = [[718, 719], [0, 1, 2]]
modrsplit = []
for i in set(run):
added = []
for j in range(len(run)):
if i == run[j]:
added.append(modr[j])
modrsplit.append(added)
# make array of filenames to use
filenames = []
for r in set(run):
fileyear = 1975 + 2*(r-1)
filenames.append('specout_daily-2yr_'+runname+'_r'+str(r)+'_startYr'+str(fileyear)+'.out')
# full filename includes directory, 'U200' means it can be up to 200 characters long
filenames_long = np.empty_like(filenames, dtype='U200')
for i in range(len(filenames)):
filenames_long[i] = dirname + str(filenames[i])
return filenames_long, modrsplit
def get_specout_data(runname, species, day, alt, lat, headdir='/home/jesse/storage/atmoRuns/'):
# Get data to plot from specout.out files
variablesokay = check_variables(species, day, alt, lat)
if isinstance(variablesokay, str):
print(variablesokay)
return None
else:
si, day, alt, lat = variablesokay
filenames_long, modrsplit = get_filenames(runname, day, headdir)
# loop over input files and read in data to 'outdata' array
hlines = 3617 # number of lines in header of specout
daylines = 297618 # number of lines of data per day in specout
outdataflat = []
for ifile in range(len(filenames_long)):
# get all necessary lines in file
# use modrsplit[ifile] instead of day
lns = []
for d in range(len(modrsplit[ifile])):
for a in range(len(alt)):
for l in range(len(lat)):
#calculate line number
ln = hlines
ln += daylines * modrsplit[ifile][d] + 2
ln += 3916 * alt[a] + 1
ln += 87 * lat[l] + 1
ln += si
lns.append(ln)
# go through file, extract all relevant values
# store into outdataflat, then reshape the array afterwards
# print('extracting', len(lns), 'values from', filenames[ifile])
with open(filenames_long[ifile]) as f:
for i, line in enumerate(f):
if i < lns[0]:
# if i is low, don't check if it's in lns
continue
else:
# add data to outdataflat and remove element from lns
outdataflat.append(float(line))
lns.remove(i)
if len(lns) == 0:
# no more values, so stop looping through file
break
outdatashape = (len(day), len(alt), len(lat))
outdata = np.array(outdataflat).reshape(outdatashape)
# print('ending function and returning data')
return outdata
def get_specoutnc_data(runname, species, day, headdir='/home/jesse/storage/atmoRuns/'):
# Get data to plot from specout.out.nc files
# Check input variables
# since we don't do alt or lat, put in dummy arrays
variablesokay = check_variables(species, day, [0], [0])
if isinstance(variablesokay, str):
print(variablesokay)
return None
else:
si, day, alt, lat = variablesokay
# Get filenames
filenames_long, modrsplit = get_filenames(runname, day, headdir)
filenames_long = [i+'.nc' for i in filenames_long]
# Extract data from nc file(s)
outdata = []
for ifile in range(len(filenames_long)):
file2read = netcdf.NetCDFFile(filenames_long[ifile],'r')
temp = file2read.variables[species]
data = temp[:]*1
file2read.close()
for d in modrsplit[ifile]:
outdata.append(data[d])
# need to swap last to axes, e.g., (6, 45, 76) -> (6, 76, 45)
outdata = np.swapaxes(outdata, 1, 2)
return outdata
def calc_z_specoutnc(runname, day, headdir='/home/jesse/storage/atmoRuns/'):
# from 'input_9hb.f' and 'setdaily_9hb.f'
# make sure 'day' is an okay value, ignore the other parts of check_variables
_, day, _, _ = check_variables('Temp', day, [0], [0])
# get temperature values for given days
mytemp = get_specoutnc_data(runname, 'Temp', day, headdir)
mytemp = np.swapaxes(mytemp, 1, 2)
# calculate basic cell grid
zalte = np.concatenate((np.linspace(0, 60, 60, endpoint=False), np.linspace(60, 116, 29, endpoint=True)))
zalt = np.zeros(len(zalte)-1)
for i in range(len(zalte)-1):
zalt[i] = (zalte[i] + zalte[i+1]) / 2.0
press = 1013.0 * np.exp(-zalt/7.0)
presse = 1013.0 * np.exp(-zalte/7.0)
# dp is change in pressure
dp = np.zeros(len(zalte)-1)
for i in range(len(dp)):
dp[i] = np.log(presse[i+1] / presse[i])
# now for stuff from setdaily
# zgp means geopotential height
# NOTE: commented out if/else populates zgp, which is the height
# It's commented out because it increases runtime by ~2-3x
# and we pretty much never need actual height more precise than zkm90arr
zgp = np.zeros((len(day), 45, 76))
mydeltaz = np.zeros_like(zgp)
RE = 6.371e3
# do the actual calculation of deltaz and z
for d in range(len(zgp)):
for ij in range(len(zgp[0])):
for ik in range(len(zgp[0][0])):
H = 29.3 * mytemp[d][ij][ik] * 1.0e-3
# if ik == 0:
# zgp[d][ij][ik] = H * np.log(1013.0/press[0])
# else:
# zgp[d][ij][ik] = 29.3 * (mytemp[d][ij][ik-1]+mytemp[d][ij][ik])/2.0 * 1.0e-3 * np.log(press[ik-1]/press[ik]) + zgp[d][ij][ik-1]
mydeltaz[d][ij][ik] = -H * dp[ik]
zkm = RE * (1.0 / (1.0 - zgp/RE)-1.0)
for ij in range(len(zgp)):
for ik in range(len(zgp[0])):
# print('ij, ik', ij, ik, mydeltaz[0][ij][ik], deltaz[0][ij][ik])
continue
return zkm, mydeltaz
# Plot characteristics of a single run
var = 'LiqH2O'
time = [8640-360, 8640-270, 8640-180, 8640-90]
# time = np.arange(8640-360, 8640, 10)
# For PMCs, do days 6685 and/or 7043
# For mesopause temp min/max, do days 6882, 7052
# For ozone, do 7200-270 (6930)
cmap = 'seismic_r' # 'PuBuGn' 'RdBu_r' 'gray_r' 'seismic'
plt.rcParams.update({'font.size': 14})
runname = 'H2Otop_coldcloud3'
runname2 = CONTROLRUN
# 2 yr increments: 0, 720, 1440, 2160, 2880, 3600, 4320, 5040, 5760, 6480, 7200, 7920, 8640
def unitlabel(var):
if var == 'Temp':
name, unit = 'Temperature', ' [K]'
elif var == 'UbarD':
name, unit = 'Mean zonal wind speed', ' [m/s]'
elif var == 'HeatOut':
name, unit = 'Heating rate', ' [K/day]'
elif var == 'CoolOut':
name, unit = 'Cooling rate', ' [K/day]'
elif var == 'W':
name, unit = 'Vertical velocity', ' [cm/s]'
elif var == 'V':
name, unit = 'Horizontal velocity (north-south)', ' [cm/s]'
elif var == 'H2O_rainout':
name, unit = 'H2O rainout', ' [molec/cm$^3$/s]'
elif var == 'Age':
name, unit = 'Clock tracer (age)', ' [???]'
elif var == 'OzoneDiag':
name, unit = 'Ozone diagnostic', ' [DU/km]'
else:
name, unit = var + ' concentration', ' [molec/cm$^3$]'
return name, unit
# returns name, unit
def find_minmax_temp(temparr):
# finds all local minima and maxima temperature values
# used to find boundaries between atmospheric layers (e.g., mesopause)
# get indices of min/max values
# localmax, localmin = [], []
# for i in range(len(lat4arr)):
# profile = temparr[i]
# a = argrelextrema(profile, np.greater)[0]
# b = argrelextrema(profile, np.less)[0]
# localmax.append([zalt90arr[x] for x in a])
# localmin.append([zalt90arr[x] for x in b])
#####################################
# try another method
localmax, localmin = [], []
for i in range(len(lat4arr)):
profile = temparr[i]
latmin, latmax = [], []
# # don't interpolate, stay true to your values!
# for j in range(len(profile)):
# if j in [0, 1, len(profile)-1, len(profile)-2]:
# continue
# if profile[j-2] < profile[j] > profile[j+2]:
# latmax.append(zalt90arr[j])
# if profile[j-2] > profile[j] < profile[j+2]:
# latmin.append(zalt90arr[j])
# do a spline!
altinterp = np.linspace(1,90,10000)
tck = interpolate.splrep(zalt90arr, profile, s=0.5, k=3)
spline = interpolate.BSpline(*tck)(altinterp)
profile = [x for x in spline]
for j in range(len(profile)):
if j in [0, len(profile)-1]:
continue
if profile[j-1] < profile[j] > profile[j+1]:
latmax.append(altinterp[j])
if profile[j-1] > profile[j] < profile[j+1]:
latmin.append(altinterp[j])
localmax.append(latmax)
localmin.append(latmin)
print('done')
return localmin, localmax
def make_altlat_plot(runname, var, time, operation=None, logColor=False, temp_minmax=False, save=False):
linColor_low, linColor_high = -100, 15
logColor_low, logColor_high = 2, 5
ylim_lower, ylim_upper = 0, 90
numbins = 24
cbarformat='%.1f'
contourlines=True
index = cname.index(var)
operationlist = ['data', 'ratio', 'percentage', 'difference']
if operation not in operationlist:
print('Unknown operation type, assuming \'data\'. ')
operation = 'data'
if operation == operationlist[0]:
# data
varplot = get_specoutnc_data(runname, var, time)[0]
titlestring = runname
elif operation == operationlist[1]:
# ratio
varplot1 = get_specoutnc_data(runname, var, time)[0]
varplot2 = get_specoutnc_data(CONTROLRUN, var, time)[0]
varplot = varplot1 / varplot2
titlestring = runname + ' vs ' + runname2 + ' ' + operation
elif operation == operationlist[2]:
# percentage change
varplot1 = get_specoutnc_data(runname, var, time)[0]
varplot2 = get_specoutnc_data(CONTROLRUN, var, time)[0]
varplot = (varplot1 - varplot2)/varplot2 * 100
titlestring = runname + ' vs ' + runname2 + ' ' + operation
elif operation == operationlist[3]:
# difference
varplot1 = get_specoutnc_data(runname, var, time)[0]
varplot2 = get_specoutnc_data(CONTROLRUN, var, time)[0]
varplot = varplot1 - varplot2
titlestring = runname + ' vs ' + runname2 + ' ' + operation
# print min and max values of plot (ignore location)
print('day {:d}, min ='.format(time), np.amin(varplot), ' max = {:.3e}'.format(np.amax(varplot)))
f, ax = plt.subplots(figsize=(9,6))
# set colorbar, define over and under colors
cmapcopy = cm.get_cmap(cmap).copy().with_extremes(over='xkcd:cyan', under='xkcd:pink')#'#ffd8d8', under='k')
# if 'difference', then center the colorbar where difference is 0
if linColor_low < 0 and linColor_high > 0:
split_across_0 = True
else:
split_across_0 = False
if (operation == 'difference' or operation == 'percentage') and split_across_0:
divnorm = colors.TwoSlopeNorm(vcenter=0.0, vmin=linColor_low, vmax=linColor_high)
else:
divnorm = None
# Make the actual plot
if logColor:
levels = np.logspace(logColor_low, logColor_high, numbins)
cs = plt.contourf(lat4arr, zalt90arr, varplot, cmap=cmapcopy,
norm=colors.LogNorm(), locator=ticker.LogLocator(), levels=levels, extend='both')
else:
levels = np.linspace(linColor_low, linColor_high, numbins, endpoint=True)
cs = plt.contourf(lat4arr, zalt90arr, varplot, cmap=cmapcopy,
norm=divnorm, levels=levels, extend='both')
if contourlines:
cs2 = plt.contour(lat4arr, zalt90arr, varplot, colors='k', linewidths=0.3, linestyles='-',
levels=levels)
# Muck around with the colorbar label
# cbarticks = [0, -2, -4, -6, -8, -10, -12]#np.linspace(0, 12, 7).tolist()
cbar = plt.colorbar(cs, format=cbarformat)#, ticks=cbarticks)
# cbar.ax.set_yticklabels([str(c) for c in cbarticks])
cbarname, cbarunit = unitlabel(var)
if operation == 'ratio':
cbarunit = ' ratio'
elif operation == 'difference':
cbarunit = ' difference'+cbarunit
elif operation == 'percentage':
cbarname = cbarname[0:len(cbarname)-13] + 'change'
cbarunit = ' [%]'
cbar.set_label(cbarname+cbarunit)
# cbar.set_label('Temperature change [K]')
if contourlines:
cbar.add_lines(cs2)
# call find_minmax_temp and plot points, if desired
# this will plot locations of minmax values, used for highlighting stratopause, mesopause, etc.
if temp_minmax:
t_index = cname.index('Temp')
temparr = get_specoutnc_data(runname, var, time)
localmin, localmax = find_minmax_temp(temparr)
# plot minmax points
for i in range(len(lat4arr)):
for j in range(len(localmax[i])):
plt.scatter(lat4arr[i], localmax[i][j], c='r', edgecolors='k')
for j in range(len(localmin[i])):
plt.scatter(lat4arr[i], localmin[i][j], c='b', edgecolors='k')
# Final plotting steps
# plt.text(0.02, 0.98, 'Day {0:d}'.format(int(time)-7200+360), \
# size=20, ha='left', va='top', transform=ax.transAxes)
plt.xlim([-90, 90])
plt.ylim(bottom=ylim_lower, top=ylim_upper)
# plt.title(titlestring+'\n'+var+' concentration on day '+str(time))
plt.xlabel('Latitude')
plt.ylabel('Altitude [km]')
xticklocs = ax.get_xticks()
xlabels = [str(int(abs(i)))+'S' if i < 0 else 0 if i == 0 else str(int(i))+'N' for i in xticklocs]
plt.xticks(xticklocs, xlabels)
plt.xlim([-90, 90])
plt.text(-80, 42, 'day '+str(time-8640+360), color='k', fontweight='bold', fontsize=16)
if save:
plt.savefig(runname+'_'+var+'_altlat_day'+str(time)+str(operation)+'.png', bbox_inches='tight', dpi=200)
# plt.close()
plt.show()
return
def make_colden_time_plot(runname, var, starttime, endtime, runname2=CONTROLRUN, operation=None, \
zlow=0.0, zhigh=100.0, save=False):
linColor_low, linColor_high = -4, 4
numbins = 17
cbarformat = '%.1f'
vlines_time = False
ylim_lower, ylim_upper = -88, 88
contourlines = True
index = cname.index(var)
timearr = np.arange(starttime, endtime)
operationlist = ['data', 'ratio', 'percentage', 'difference']
cbarlevels = np.linspace(linColor_low, linColor_high, numbins)
# get data and deltaz values
wholedata = get_specoutnc_data(runname, var, timearr)
varplot = np.zeros((len(wholedata), len(wholedata[0][0])))
_, mydeltaz = calc_z_specoutnc(runname, timearr)
# get relevant concentrations and deltaz values
# multiply them together using the dot product, 1e5 converts cm to km
for t in range(len(wholedata)):
for l in range(len(wholedata[0][0])):
concs = wholedata[t].T[l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
dzvals = mydeltaz[t][l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
varplot[t][l] = np.dot(concs, dzvals) *1.0e5
# do the same calculation for the control run
wholedata2 = get_specoutnc_data(runname2, var, timearr)
varplot2 = np.zeros((len(wholedata2), len(wholedata2[0][0])))
_, mydeltaz2 = calc_z_specoutnc(runname2, timearr)
for t in range(len(wholedata2)):
for l in range(len(wholedata2[0][0])):
concs2 = wholedata2[t].T[l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
dzvals2 = mydeltaz2[t][l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
varplot2[t][l] = np.dot(concs2, dzvals2) *1.0e5
if operation not in operationlist:
print('Unknown operation type, assuming \'data\'. ')
operation = 'data'
if operation == operationlist[0]:
# data
varplot = varplot
titlestring = runname
elif operation == operationlist[1]:
# ratio
varplot = varplot / varplot2
titlestring = runname + ' vs ' + runname2 + ' ' + operation
elif operation == operationlist[2]:
# percentage change
varplot = (varplot - varplot2) / varplot2 * 100.0
titlestring = runname + ' vs ' + runname2 + ' ' + operation
elif operation == operationlist[3]:
# difference
varplot = varplot - varplot2
titlestring = runname + ' vs ' + runname2 + ' ' + operation
# print min and max values of plot (ignore location)
print('min = {:.2e} max = {:.2e}'.format(np.amin(varplot), np.amax(varplot)))
# make plot
cmapcopy = cm.get_cmap(cmap).copy().with_extremes(over='k', under='w')#'0.8')
# if 'difference', then center the colorbar where difference is 0
if linColor_low < 0 and linColor_high > 0:
split_across_0 = True
else:
split_across_0 = False
if (operation == 'difference' or operation == 'percentage') and split_across_0:
divnorm = colors.TwoSlopeNorm(vcenter=0.0, vmin=linColor_low, vmax=linColor_high)
else:
divnorm = None
f, ax = plt.subplots(figsize=(14,8))
cs = plt.contourf(timearr, lat4arr, varplot.T,
levels=cbarlevels, norm=divnorm,
cmap=cmapcopy, extend='both')
if contourlines:
cs2 = ax.contour(timearr, lat4arr, varplot.T, colors='k', linewidths=0.3, linestyles='-',
levels=cbarlevels)
cbar = plt.colorbar(cs, format=cbarformat, pad=0.02)
if var == 'O3':
# NOTE: conversion to DU not actually implemented!
cbarname = var
cbarunit = ' concentration'# [DU]'
else:
cbarname, cbarunit = unitlabel(var)
cbarunit = list(cbarunit)
if cbarunit[-3] == '3':
cbarunit[-3] = '2'
else:
cbarunit.insert(-1, ' * cm')
cbarunit = ''.join(cbarunit)
cbar.set_label(cbarname+cbarunit)
if operation == 'ratio':
cbarunit = ' ratio'
elif operation == 'percentage':
cbarunit = ' difference [%]'
cbar.set_label(cbarname+cbarunit)
# cbar.set_label('H$_2$O ice column density [molec/cm$^2$]')
cbar.set_label('Change in ozone column density [%]')
if contourlines:
cbar.add_lines(cs2)
if vlines_time:
vlinecolors = ['C0', 'C1']
for i in range(len(time)):
plt.axvline(time[i], c=vlinecolors[i])
# plt.title(titlestring+'\ncolumn density'+', from '+str(zlow)+' km to '+str(zhigh)+' km')
plt.xlabel('Day of year')
plt.ylabel('Latitude')
plt.ylim(bottom=ylim_lower, top=ylim_upper)
yticklocs = ax.get_yticks()
ylabels = [str(int(abs(i)))+'S' if i < 0 else 0 if i == 0 else str(int(i))+'N' for i in yticklocs]
plt.yticks(yticklocs, ylabels)
plt.ylim(bottom=ylim_lower, top=ylim_upper)
# Option to change the x-axis labels to something more sensible
# xlocs = np.linspace(8640-360, 8640, 13)
# xlabels = ['January', 'April', 'July', 'October', ' ']
# xlabels = ['J', 'F', 'M', 'A', 'M', 'J', 'J', 'A', 'S', 'O', 'N', 'D', ' ']
# xlabels = np.arange(0, 361, 30)
# plt.xticks(xlocs, xlabels)
if save:
plt.savefig(runname+'_'+var+'_colden'+'.png', bbox_inches='tight', dpi=200)
plt.show()
return
def make_colden_globalavg_plot(runname, var, starttime, endtime, zlow=0.0, zhigh=100.0, \
runname2=CONTROLRUN, operation=None, save=False):
# start very similar to make_colden_time_plot(),
# then do a (weighted) average over latitude
index = cname.index(var)
timearr = np.arange(starttime, endtime)
operationlist = ['data', 'ratio', 'percentage', 'difference']
# get data and deltaz values
wholedata = get_specoutnc_data(runname, var, timearr)
varplot = np.zeros((len(wholedata), len(wholedata[0][0])))
_, mydeltaz = calc_z_specoutnc(runname, timearr)
# get relevant concentrations and deltaz values
# multiply them together using the dot product, 1e5 converts cm to km
# Note that varplot is the 2D column density
for t in range(len(wholedata)):
for l in range(len(wholedata[0][0])):
concs = wholedata[t].T[l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
dzvals = mydeltaz[t][l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
varplot[t][l] = np.dot(concs, dzvals) * 1.0e5
# do the same calculation for the control run
wholedata2 = get_specoutnc_data(CONTROLRUN, var, timearr)
varplot2 = np.zeros((len(wholedata2), len(wholedata2[0][0])))
_, mydeltaz2 = calc_z_specoutnc(CONTROLRUN, timearr)
for t in range(len(wholedata2)):
for l in range(len(wholedata2[0][0])):
concs2 = wholedata2[t].T[l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
dzvals2 = mydeltaz2[t][l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
varplot2[t][l] = np.dot(concs2, dzvals2) * 1.0e5
# if looking at ozone (O3), convert to Dobson units
if var == 'O3':
varplot = varplot/2.687e16
varplot2 = varplot2/2.687e16
if operation not in operationlist:
print('Unknown operation type, assuming \'data\'. ')
operation = 'data'
if operation == operationlist[0]:
# data
varplot = varplot
titlestring = runname
elif operation == operationlist[1]:
# ratio
varplot = varplot / varplot2
titlestring = runname + ' vs ' + runname2 + ' (' + operation + ')'
elif operation == operationlist[2]:
# percentage change
varplot = (varplot - varplot2) / varplot2 * 100.0
titlestring = runname + ' vs ' + runname2 + ' (' + operation + ')'
elif operation == operationlist[3]:
# difference
varplot = varplot - varplot2
titlestring = runname + ' vs ' + runname2 + ' (' + operation + ')'
# Instead of a simple average of latitude values,
# we need to weight by their surface area they cover.
# The area of a band on a sphere is A = 2*pi*R^2*(sin(phi1)-sin(phi2))
# Then to normalize, divide by total surface area, 4*pi*R^2
# So in short, the weighting factor is 0.5*(sin(phi1)-sin(phi2))
latweight = np.zeros_like(lat4arr)
for i in range(len(latweight)):
radlow = (lat4arr[i] - 2.0) * np.pi / 180.
radhigh = (lat4arr[i] + 2.0) * np.pi / 180.
latweight[i] = 0.5 * (np.sin(radhigh) - np.sin(radlow))
globalavg = np.zeros(len(varplot))
for i in range(len(varplot)):
weighted = np.dot(varplot[i], latweight)
globalavg[i] = weighted
if var == 'O3':
yname, yunit = var+' concentration', ' [DU]'
else:
yname, yunit = unitlabel(var)
yunit = list(yunit)
if yunit[-3] == '3':
yunit[-3] = '2'
else:
yunit.insert(-1, ' * cm')
yunit = ''.join(yunit)
plt.figure(figsize=(8,8))
plt.plot(timearr, globalavg)
plt.xlim(left=timearr[0], right=timearr[-1])
plt.xlabel('Time [days]')
if operation == 'ratio':
yunit = ' ratio'
elif operation == 'percentage':
yunit = ' [%]'
plt.ylabel(yname+yunit)
plt.title(titlestring+'\n'+var+', globally-averaged, from '+str(zlow)+' km to '+str(zhigh)+' km')
plt.axvline(360, c='k', linestyle='--') # start of cloud
plt.axvline(5760, c='k', linestyle='--') # end of cloud
if save:
plt.savefig(runname+'_'+var+'_colden_globalavg'+'.png', bbox_inches='tight', dpi=200)
plt.show()
return
def make_altprofile_plot(runname, var, time, lat, runname2=CONTROLRUN, plotcontrol=False, \
showlegend=True, logx=True, save=False):
# lat can be either a single number or a list of values
xbottom, xtop = 150, 300
ybottom=0
index = cname.index(var)
if not hasattr(lat, '__len__'):
lat = [lat]
degreevals = [int(lat4arr[i]) for i in lat]
# get data
varplot2d = get_specoutnc_data(runname, var, time)[0].T
varplot = np.zeros((len(lat), len(zalt90arr)))
for i in range(len(lat)):
for j in range(len(zalt90arr)):
varplot[i][j] = varplot2d[lat[i]][j]
# get data from second run
varplot22d = get_specoutnc_data(runname2, var, time)[0].T
varplot2 = np.zeros((len(lat), len(zalt90arr)))
for i in range(len(lat)):
for j in range(len(zalt90arr)):
varplot2[i][j] = varplot22d[lat[i]][j]
colors = ['C0', 'C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7', 'C8', 'C9']
for i in range(len(varplot)):
labelname = runname + ', ' + str(degreevals[i]) + '$\degree$'
plt.plot(varplot[i], zalt90arr, '-', c=colors[i], label=labelname)#'{:n}$\degree$'.format(lat[i]))
if plotcontrol:
for i in range(len(varplot)):
labelname = runname2 + ', ' + str(degreevals[i]) + '$\degree$'
plt.plot(varplot2[i], zalt90arr, '--', c=colors[i], label=labelname)#'C {:n}$\degree$'.format(lat[i]))
if logx:
plt.xscale('log')
plt.title(runname+'\n'+var+' profile plot on day '+str(time)+', lat {}$\degree$'.format(degreevals))
plt.xlabel(var+unitlabel(var)[1])
plt.ylabel('Altitude [km]')
plt.xlim(left=xbottom, right=xtop)
plt.ylim(bottom=ybottom)
if showlegend:
plt.legend()
if save:
plt.savefig(runname+'_'+var+'_altprofile_day_'+str(time)+'.png', bbox_inches='tight', dpi=200)
plt.show()
return
if isinstance(time, int):
time = [time]
# for t in time:
# make_altlat_plot(runname, 'O3', t, operation='percentage', logColor=False, \
# temp_minmax=False, save=False)
# make_colden_time_plot(runname, 'LiqH2O', 8640-360, 8640, zlow=60.0, zhigh=100.0, \
# operation='data', save=False)
make_colden_time_plot(runname, 'O3', 8640-360, 8640, zlow=0.0, zhigh=100.0, \
operation='percentage', save=False)
# make_colden_time_plot(chemarr, 'LiqH2O', 0, 5760, zlow=60.0, save=False)
# make_colden_globalavg_plot(chemarr, 'O3', 0, 2160, zlow=0.0, chemarr2=chemarr2, operation='percentage', \
# save=False)
# for t in time:
# make_altprofile_plot(chemarr, 'H2O', t, [84], \
# chemarr2=chemarr2, plotcontrol=True, logvar=True, differences=True, save=False)
/home/jesse/anaconda3/lib/python3.9/site-packages/scipy/io/netcdf.py:306: RuntimeWarning: Cannot close a netcdf_file opened with mmap=True, when netcdf_variables or arrays referring to its data still exist. All data arrays obtained from such files refer directly to data on disk, and must be copied before the file can be cleanly closed. (See netcdf_file docstring for more information on mmap.) warnings.warn((
min = -4.40e+00 max = 3.84e+00
# Make 5-panel plot of ozone depletion for paper.
# Copy-paste "make_altlat_plot()" and modify slightly to produce all plots at once
# Copy-paste "make_colden_plot()" for bottom panel
time = [8640-360+42, 8640-360+212]
runnames = ['H2Otop_coldcloud4', 'H2Otop_coldcloud3']
plt.rcParams.update({'font.size': 14})
logColor_low, logColor_high = -1, 2
linColor_low, linColor_high = -100, 0
ylim_lower, ylim_upper = 40, 90
numbins = 21
cbarformat='%.1f'
cmap = 'jet_r'
var='O3'
index = cname.index(var)
# fig, axs = plt.subplots(2, 2, figsize=(12,8), sharex=False, sharey=True)
fig = plt.figure(figsize=(16,12))
ax1 = plt.subplot2grid(shape=(3, 3), loc=(0, 0), colspan=1)
ax2 = plt.subplot2grid(shape=(3, 3), loc=(0, 1), colspan=1)
ax3 = plt.subplot2grid(shape=(3, 3), loc=(1, 0), colspan=1)
ax4 = plt.subplot2grid(shape=(3, 3), loc=(1, 1), colspan=1)
ax5 = plt.subplot2grid(shape=(3, 3), loc=(2, 0), colspan=2)
axs = [ax1, ax2, ax3, ax4, ax5]
# set colorbar, define over and under colors
cmapcopy = cm.get_cmap(cmap).copy().with_extremes(over='w', under='k')
# Make the actual plot
levels = np.linspace(linColor_low, linColor_high, numbins, endpoint=True)
mytext = [['a','b'],['c','d']]
myruns = ['LB', 'LxCC']
mydays = ['Day '+str(time[0]-8640+360), 'Day '+str(time[1]-8640+360)]
for i in range(2):
for j in range(2):
varplot1 = get_specoutnc_data(runnames[i], var, time[j])[0]
varplot2 = get_specoutnc_data(CONTROLRUN, var, time[j])[0]
varplot = (varplot1 - varplot2)/varplot2 * 100
axindex = 2*i+j
cs = axs[axindex].contourf(lat4arr, zalt90arr, varplot, cmap=cmapcopy,
levels=levels, extend='max')
cs2 = axs[axindex].contour(lat4arr, zalt90arr, varplot, colors='k', linewidths=0.3, linestyles='-',
levels=levels)
tcolor = 'w'
axs[axindex].text(0.05, 0.07, mytext[i][j], transform=axs[axindex].transAxes, c=tcolor, weight='bold', fontsize=20,
path_effects=[pe.withStroke(linewidth=2, foreground='k')])
axs[axindex].text(0.7, 0.12, myruns[i], transform=axs[axindex].transAxes, c=tcolor, fontsize=16,
path_effects=[pe.withStroke(linewidth=2, foreground='k')])
axs[axindex].text(0.7, 0.04, mydays[j], transform=axs[axindex].transAxes, c=tcolor, fontsize=16,
path_effects=[pe.withStroke(linewidth=2, foreground='k')])
print('min', mytext[i][j], np.amin(varplot))
plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0.1, hspace=0.05)
# Colorbar
fig.subplots_adjust(right=0.88)
cbar_ax = fig.add_axes([0.63, 0.37, 0.02, 0.53])
cbar = fig.colorbar(cs, cax=cbar_ax)
cbar.add_lines(cs2)
cbar.set_label('Change in ozone concentration [%]')
# Final plotting steps
xticklocs = np.arange(-90, 91, 30)
xlabels = [str(int(abs(i)))+'S' if i < 0 else 0 if i == 0 else str(int(i))+'N' for i in xticklocs]
for ax in range(len(np.array(axs).flat[0:4])):
axs[ax].axis(xmin=-90, xmax=90, ymin=ylim_lower, ymax=ylim_upper)
# label y axes
if ax in [0,2]:
axs[ax].set(ylabel='Altitude [km]')
else:
axs[ax].yaxis.set_major_locator(ticker.NullLocator())
# label x axes
if ax in [2,3]:
axs[ax].set(xlabel='Latitude')
axs[ax].set_xticks(xticklocs, xlabels)
else:
axs[ax].xaxis.set_major_locator(ticker.NullLocator())
######################################################
# Now make a panel with the total ozone column density
box = axs[4].get_position()
box.y0 = box.y0 - 0.1
box.y1 = box.y1 - 0.05
axs[4].set_position(box)
runname = 'H2Otop_coldcloud3'
var = 'O3'
starttime, endtime = 8640-360, 8640
zlow, zhigh = 0.0, 100.0
linColor_low, linColor_high = -4, 4
numbins = 17
cbarformat = '%.1f'
vlines_time = False
ylim_lower, ylim_upper = -90, 90
contourlines = True
cmap = 'seismic_r'
index = cname.index(var)
timearr = np.arange(starttime, endtime)
cbarlevels2 = np.linspace(linColor_low, linColor_high, numbins)
# get data and deltaz values
wholedata = get_specoutnc_data(runname, var, timearr)
varplot = np.zeros((len(wholedata), len(wholedata[0][0])))
_, mydeltaz = calc_z_specoutnc(runname, timearr)
# get relevant concentrations and deltaz values
# multiply them together using the dot product, 1e5 converts cm to km
for t in range(len(wholedata)):
for l in range(len(wholedata[0][0])):
concs = wholedata[t].T[l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
dzvals = mydeltaz[t][l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
varplot[t][l] = np.dot(concs, dzvals) *1.0e5
# do the same calculation for the control run
wholedata2 = get_specoutnc_data(runname2, var, timearr)
varplot2 = np.zeros((len(wholedata2), len(wholedata2[0][0])))
_, mydeltaz2 = calc_z_specoutnc(runname2, timearr)
for t in range(len(wholedata2)):
for l in range(len(wholedata2[0][0])):
concs2 = wholedata2[t].T[l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
dzvals2 = mydeltaz2[t][l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
varplot2[t][l] = np.dot(concs2, dzvals2) *1.0e5
# calculate percentage
varplot = (varplot - varplot2) / varplot2 * 100.0
# print min and max values of plot (ignore location)
print('min = {:.2e} max = {:.2e}'.format(np.amin(varplot), np.amax(varplot)))
# make plot
cmapcopy2 = cm.get_cmap(cmap).copy().with_extremes(over='k', under='w')
divnorm2 = colors.TwoSlopeNorm(vcenter=0.0, vmin=linColor_low, vmax=linColor_high)
# f, ax = plt.subplots(figsize=(14,8))
cs2 = axs[4].contourf(timearr, lat4arr, varplot.T,
levels=cbarlevels2, norm=divnorm2,
cmap=cmapcopy2, extend='both')
cs2x = axs[4].contour(timearr, lat4arr, varplot.T, colors='k', linewidths=0.3, linestyles='-',
levels=cbarlevels2)
# Colorbar
cbar_ax2 = fig.add_axes([0.63, -0.02, 0.02, 0.34])
cbar2 = fig.colorbar(cs2, cax=cbar_ax2)
cbar2.add_lines(cs2x)
cbar2.set_label('Change in ozone column [%]')
# Text, subfigure label
tcolor = 'w'
axs[4].text(0.025, 0.07, 'e', transform=axs[4].transAxes, c=tcolor, weight='bold', fontsize=20,
path_effects=[pe.withStroke(linewidth=2, foreground='k')])
# Axes ticks and labels
axs[4].set_xlabel('Day of year')
axs[4].set_ylabel('Latitude')
yticklocs = np.arange(-90, 91, 30)
ylabels = [str(int(abs(i)))+'S' if i < 0 else 0 if i == 0 else str(int(i))+'N' for i in yticklocs]
axs[4].set_yticks(yticklocs, ylabels)
axs[4].set_ylim(bottom=ylim_lower, top=ylim_upper)
# Option to change the x-axis labels to something more sensible
xlocs = np.linspace(8640-360, 8640, 13)
# xlabels = ['January', 'April', 'July', 'October', ' ']
# xlabels = ['J', 'F', 'M', 'A', 'M', 'J', 'J', 'A', 'S', 'O', 'N', 'D', ' ']
xlabels = np.arange(0, 361, 30)
axs[4].set_xticks(xlocs, xlabels)
# plt.savefig('OzoneDepletion_5panel.png', bbox_inches='tight', dpi=200, transparent=False)
plt.show()
min a -99.52013 min b -97.83116 min c -99.98737 min d -99.88817 min = -4.40e+00 max = 3.84e+00
# Make 2-panel temperature difference plot
time = [8640-360+30, 8640-360+210]#, 8640-360+120, 8640-360+180, 8640-360+240, 8640-360+300]
# time = np.arange(8640-360, 8640, 60)
print(len(time), 'plots')
var = 'Temp'
plt.rcParams.update({'font.size': 14})
logColor_low, logColor_high = -1, 2
linColor_low, linColor_high = -12, 12
ylim_lower, ylim_upper = 40, 90
numbins = 25
cbarformat='%.1f'
cmap = 'seismic'
contourlines=True
index = cname.index(var)
fig, axs = plt.subplots(1, 2, figsize=(10,4))
# set colorbar, define over and under colors
cmapcopy = cm.get_cmap(cmap).copy().with_extremes(over='xkcd:pink', under='xkcd:cyan')#'#ffd8d8', under='k')
# center the colorbar where difference is 0
if linColor_low < 0 and linColor_high > 0:
split_across_0 = True
else:
split_across_0 = False
if split_across_0:
divnorm = colors.TwoSlopeNorm(vcenter=0.0, vmin=linColor_low, vmax=linColor_high)
else:
divnorm = None
# array of the alphabet
mytext = [chr(i) for i in range(ord('a'),ord('z')+1)][0:len(time)]
for i in range(len(axs)):
# for j in range(len(axs[0])):
t_index = i#*(len(axs[0])) + j
# Make the actual plot
varplot1 = get_specoutnc_data('H2Otop_coldcloud3', var, time[t_index])[0]
varplot2 = get_specoutnc_data(CONTROLRUN, var, time[t_index])[0]
varplot = varplot1 - varplot2
levels = np.linspace(linColor_low, linColor_high, numbins, endpoint=True)
cs = axs[i].contourf(lat4arr, zalt90arr, varplot, cmap=cmapcopy,
norm=divnorm, levels=levels, extend='both')
if contourlines:
cs2 = axs[i].contour(lat4arr, zalt90arr, varplot, colors='k', linewidths=0.3, linestyles='-',
levels=levels)
tcolor='k'
axs[i].text(0.05, 0.05, mytext[t_index], transform=axs[i].transAxes, c='w', weight='bold', fontsize=20,
path_effects=[pe.withStroke(linewidth=2, foreground='k')])
axs[i].text(0.7, 0.05, 'day '+str(int(time[t_index]-8640+360)), transform=axs[i].transAxes, c='w', fontsize=16,
path_effects=[pe.withStroke(linewidth=2, foreground='k')])
plt.subplots_adjust(left=0.1, right=0.9, wspace=0.05, hspace=0.05)
# Colorbar
fig.subplots_adjust(right=0.88)
cbar_ax = fig.add_axes([0.90, 0.1, 0.02, 0.8])
cbar = fig.colorbar(cs, cax=cbar_ax)
cbar.add_lines(cs2)
cbar.set_label('Temperature difference [K]')
# Final plotting steps
# xticklocs = axs[0].get_xticks()
xticklocs = np.arange(-60, 91, 30)
xlabels = [str(int(abs(i)))+'S' if i < 0 else 0 if i == 0 else str(int(i))+'N' for i in xticklocs]
for ax in axs.flat:
ax.axis(xmin=-90, xmax=90, ymin=ylim_lower, ymax=ylim_upper)
ax.set(xlabel='Latitude', ylabel='Altitude [km]')
ax.set_xticks(xticklocs, xlabels)
ax.label_outer()
# plt.savefig('TempDiff_LxCCv3.png', bbox_inches='tight', dpi=200)
plt.show()
2 plots
# Make 2-panel ice column density plot
plt.rcParams.update({'font.size': 14})
linColor_low, linColor_high = 0, 3e15 # will get overwritten in loop
numbins = 16
cbarformat = '%.1e'
ylim_lower, ylim_upper = 50, 88
zlow, zhigh = 60, 100
contourlines = True
var = 'LiqH2O'
index = cname.index(var)
timearr = np.arange(8640-360, 8640)
operationlist = ['data', 'ratio', 'percentage', 'difference']
cmap = 'jet'
cbarlevels = np.linspace(linColor_low, linColor_high, numbins) # also gets overwritten
myruns = [CONTROLRUN, 'H2Otop_coldcloud4', 'H2Otop_coldcloud3']
mytext1 = ['a', 'b', 'c']
mytext2 = ['control', 'LB', 'LxCC']
cmapcopy = cm.get_cmap(cmap).copy().with_extremes(over='k', under='w')#'0.8')
fig, axs = plt.subplots(len(myruns), 1, figsize=(12,12), sharex=True)
for i in range(len(myruns)):
# reset colorbars since they're on two different scales
if i == 0:
linColor_low, linColor_high = 0.0, 6e13
cbarlevels = np.linspace(linColor_low, linColor_high, numbins)
elif i == 1:
linColor_low, linColor_high = 0.0, 3e14
cbarlevels = np.linspace(linColor_low, linColor_high, numbins)
elif i == 2:
linColor_low, linColor_high = 0.0, 3e15
cbarlevels = np.linspace(linColor_low, linColor_high, numbins)
# calculate column density, put it into varplot
wholedata = get_specoutnc_data(myruns[i], var, timearr)
varplot = np.zeros((len(wholedata), len(wholedata[0][0])))
_, mydeltaz = calc_z_specoutnc(myruns[i], timearr)
for t in range(len(wholedata)):
for l in range(len(wholedata[0][0])):
concs = wholedata[t].T[l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
dzvals = mydeltaz[t][l][np.where((zalt90arr>zlow) & (zalt90arr<zhigh))]
varplot[t][l] = np.dot(concs, dzvals) *1.0e5
cs = axs[i].contourf(timearr, lat4arr, varplot.T,
levels=cbarlevels,
cmap=cmapcopy, extend='both')
if contourlines:
cs2 = axs[i].contour(timearr, lat4arr, varplot.T, colors='k', linewidths=0.3, linestyles='-',
levels=cbarlevels)
cbar = fig.colorbar(cs, ax=axs[i], format='%.1e', pad=0.02)
if contourlines:
cbar.add_lines(cs2)
cbar.set_label('NLC ice column [molec/cm$^2$]')
axs[i].text(0.03, 0.07, mytext1[i], transform=axs[i].transAxes, c='w', fontsize=20, weight='bold',
path_effects=[pe.withStroke(linewidth=2, foreground='k')])
axs[i].text(0.88, 0.07, mytext2[i], transform=axs[i].transAxes, c='w', fontsize=16,
path_effects=[pe.withStroke(linewidth=2, foreground='k')])
print('max {:.3e}'.format(np.amax(varplot)))
plt.subplots_adjust(left=0.1, right=0.9, wspace=0.05, hspace=0.1)
# Option to change the y-axis labels to something more sensible
# yticklocs = axs[0].get_yticks()
yticklocs = np.arange(50, 88, 5)
ylabels = [str(int(abs(i)))+'S' if i < 0 else 0 if i == 0 else str(int(i))+'N' for i in yticklocs]
for ax in axs.flat:
ax.axis(xmin=8640-360, xmax=8640, ymin=ylim_lower, ymax=ylim_upper)
ax.set(ylabel='Latitude')
ax.set_yticks(yticklocs, ylabels)
ax.set_xticks(np.arange(8640-360, 8641, 30), [str(i) for i in np.arange(0, 361, 30)])
ax.label_outer()
ax.axis(xmin=8640-360, xmax=8640, ymin=ylim_lower, ymax=ylim_upper)
axs[len(myruns)-1].set(xlabel='Day of year')
# plt.savefig('IceColumn_2panel.png', bbox_inches='tight', dpi=200)
plt.show()
max 4.646e+13 max 2.891e+14 max 2.688e+15