Skip to content

Commit

Permalink
Merge pull request #3 from tsherwen/dev
Browse files Browse the repository at this point in the history 
Dev
  • Loading branch information
@tsherwen
tsherwen authored Mar 1, 2022
2 parents 3dbbe35 + 3d33165 commit 3f2b9a3
Show file tree
Hide file tree
Showing 9 changed files with 9,445 additions and 115 deletions.
6,876 changes: 6,876 additions & 0 deletions arna/#plotting.py#
View file

Large diffs are not rendered by default.

10 changes: 10 additions & 0 deletions arna/GEOS.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2037,6 +2037,16 @@ def add_extra_derived_species(ds, debug=False):
# attrs['units'] = 'unitless'
# ds[NewVar].attrs = attrs

# Try adding NOx if not present
NewVar = 'NOx'
try:
ds['NOx']
except KeyError:
try:
ds[NewVar] = ds['NO'] + ds['NO2']
except KeyError:
pass

# Add NOy divided by PM2.5 concentration
NewVar = 'NOy*PM2.5(dust)'
ds[NewVar] = ds['NOy'] * ds['PM2.5(dust)']
Expand Down
80 changes: 71 additions & 9 deletions arna/GEOS_Chem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,11 +34,18 @@ def PF_TRAXXX_2TracerName(TRA_XX, folder=None, RTN_dict=False):
return dict(zip(nums, TRAs))[TRA_XX]


def get_dict_of_GEOSChem_model_output(res='0.5x0.625',
def get_dict_of_GEOSChem_model_output(res='0.5x0.625', folder4netCDF=False,
RunSet='MERRA2-0.5-initial',
CoreRunsOnly=False):
"""
Retrieve dictionary of model run names and their full location paths
Parameters
-------
folder4netCDF (bool): append subfolder for NetCDFs (OutputDir) to folder
Notes
-----
"""
RunRoot = get_local_folder('RunRoot')
if res == '4x5':
Expand Down Expand Up @@ -73,9 +80,41 @@ def get_dict_of_GEOSChem_model_output(res='0.5x0.625',
# Isotherm
# RunStr = 'DustUptake.JNIT.Isotherm.BCs'
# RunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II'
RunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags.v2'
# IsoRunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags.'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, IsoRunStr+'v2')
# d['Acid-4x5-Isotherm.v2.0'] = folder
# ... Other Isotherm versions
# RunStr = IsoRunStr + 'v2.very_low'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-Isotherm.v2.3'] = folder
# RunStr = IsoRunStr + 'v2.low'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-Isotherm.v2.2'] = folder
# RunStr = IsoRunStr + 'v2.medium'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-Isotherm.v2.1'] = folder
# RunStr = IsoRunStr + 'v2.deli'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-Isotherm.v2.0.deli'] = folder
# RunStr = IsoRunStr + 'v2.4'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-Isotherm.v2.4'] = folder
# RunStr = IsoRunStr + 'v2.4.deli'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-Isotherm.v2.4.deli'] = folder
# RunStr = IsoRunStr + 'v2.4.deli.H2O'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-Isotherm.v2.4.deli.H2O'] = folder
# Version 3
RunStr = IsoRunStr + 'v3.0.H2O'
folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
d['Acid-4x5-Isotherm.v3.0.H2O'] = folder
# Version 3 plus v3.0
RunStr = IsoRunStr + 'v3.0.H2O.Acid'
folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
d['Acid-4x5-Isotherm.v2'] = folder
d['Acid-4x5-Isotherm.v3.0.H2O.Acid'] = folder
# Re-run acid exploration?

# Isotherm + HONO 100%
# RunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags.HONO100'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
Expand All @@ -88,25 +127,44 @@ def get_dict_of_GEOSChem_model_output(res='0.5x0.625',
RunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags.J50'
folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
d['Acid-4x5-J50'] = folder
# ACID + J50 (no BB)
# RunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags.J50.BBx0'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-J50-BBx0'] = folder
# Acid plus v3.0


# Isotherm + BBx3 + NH3x3
RunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags'
RunStr += '.J50.BBx3.NH3x3'
folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
d['Acid-4x5-J50-BBx3-NH3x3'] = folder
# RunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags'
# RunStr += '.J50.BBx3.NH3x3'
# folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
# d['Acid-4x5-J50-BBx3-NH3x3'] = folder
# Isotherm + African BBx3 + African NH3x3
RunStr = 'DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags'
RunStr += '.J50.BBx3AFRICA.NH3x3/'
folder = '{}/{}{}/'.format(RunRoot, AcidRunStr, RunStr)
d['Acid-4x5-J50-AfBBx3-NH3x3'] = folder
# Temporary runs - MERRA-2 met
# RunStr = 'merra2_4x5_standard.v12.9.0.BASE.2019.2020.PF'
# folder = '{}/{}{}/'.format(RunRoot, RunStr, '')
# d['BASE-4x5-MERRA-2'] = folder

# Only return the core runs
if CoreRunsOnly:
runs2use = [
# 'Acid-4x5',
# 'Acid-4x5-J25',
'Acid-4x5-J00',
'Acid-4x5-J50',
'Acid-4x5-Isotherm.v2',
'Acid-4x5-J50-AfBBx3-NH3x3',
# 'Acid-4x5-Isotherm.v2.4',
# 'Acid-4x5-J50-AfBBx3-NH3x3',
# 'Acid-4x5-Isotherm.v2.4',
# 'Acid-4x5-Isotherm.v2',
'Acid-4x5-Isotherm.v3.0.H2O',
# 'Acid-4x5-Isotherm.v3.0.H2O.Acid',
# Temp
# 'BASE-4x5-MERRA-2',
# 'Acid-4x5-J50-BBx0',
]
dNew = {}
for run in runs2use:
Expand Down Expand Up @@ -228,6 +286,10 @@ def get_dict_of_GEOSChem_model_output(res='0.5x0.625',
d['BC-BASE'] = folder
else:
pass
# Include NetCDF subfolder in run directory strings
if folder4netCDF:
for key in d.keys():
d[key] = d[key] + '/OutputDir/'
return d


Expand Down
8 changes: 8 additions & 0 deletions arna/observations.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1132,3 +1132,11 @@ def mk_planeflight_files4FAAM_campaigns(folder=None, testing_mode=False,
gc.collect()
del ds, ds_l
gc.collect()


def get_biomass_burning_flag_for_ARNA2(flight_ID='C225'):
"""
"""
folder = get_local_folder('DataRoot') + '/Misc/BB_flag/'
files = glob.glob( folder + '*_bb_flag.csv' )
dfs = [pd.read_csv(i)]
126 changes: 117 additions & 9 deletions arna/plotting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,6 @@ def get_vmin_value4var(input):
"""
Set a vmin value for a variable (below which values are set to NaN)
"""
#
if 'Dust' in input:
vmin = 15
elif 'NOy' in input:
Expand Down Expand Up @@ -357,8 +356,8 @@ def plot_up_longitudinally_sampled_locs(ds, var2use='noy', extents=None,
if isinstance(ds, type(None)):
folder = '/Users/tomassherwen/Google_Drive/Data/ARNA/GEOS_CF/'
folder += '/data_GEOSCF_2019_12_14/'
filename = 'ARNA_GEOSCF_chm_inst_1hr_g1440x721_p23_Cape_Verde_2019_353_noy_'
filename += 'lvls_1000_900_800_700_600_500.nc'
filename = 'ARNA_GEOSCF_chm_inst_1hr_g1440x721_p23_Cape_Verde'
filename += '_2019_353_noy_lvls_1000_900_800_700_600_500.nc'
ds = xr.open_dataset(folder + filename)
# Local area analysed as Cape Verde
d = get_analysis_region('local_CVAO_area')
Expand Down Expand Up @@ -518,8 +517,8 @@ def plt_ts4ds(ds, region='Cape_Verde', extr_str='',
max_vals = da.max(dim=['lat', 'lon'], skipna=True).values
time = AC.dt64_2_dt(da.time.values)
# Now plot up as a time series
prtstr = 'WARNING: Not plotting {} @ {:.0f}hPa because all NaNs!'
if mean_vals[np.isfinite(mean_vals)].shape[0] == 0:
prtstr = 'WARNING: Not plotting {} @ {:.0f}hPa because all NaNs!'
print(prtstr.format(var, lev))
pass
else:
Expand Down Expand Up @@ -2742,18 +2741,19 @@ def plt_comp_by_alt_4ARNA_together(dpi=320, just_SLR=True, show_plot=False,
context="paper", font_scale=0.75,
NOxAsLog=False,
CoreRunsOnly=False,
debug=False):
verbose=True, debug=False):
"""
Plot up altitude binned comparisons between core obs. and model data
"""
PrtStr = "NOTE: Plotting RunSet: '{}' @ res. of: '{}' (CoreRunsOnly?: {})"
print( PrtStr.format(RunSet, res, CoreRunsOnly) )
# Setup the pdf file to use
if isinstance(savetitle, type(None)):
savetitle = 'ARNA_altitude_binned_combined_file_{}'.format(res)
if just_SLR:
savetitle += '_JUST_SLR'
pdff = AC.plot2pdfmulti(title=savetitle, open=True, dpi=dpi)


# if isinstance(RunSet, type(None)):
# RunSet='FP-Nest'
# Call the standard species plotter
Expand Down Expand Up @@ -6749,7 +6749,6 @@ def mk_Andersen2021_figure_02(dpi=720, figsize=(7, 3), aspect=None, ms=None,
df3["Enhancement"] = df3["JpNO3"]/(7*10**-7)

# Plotting settings

if isinstance(fontsize, type(None)):
# fontsize = 20
fontsize = 8
Expand Down Expand Up @@ -6794,7 +6793,7 @@ def mk_Andersen2021_figure_02(dpi=720, figsize=(7, 3), aspect=None, ms=None,
Line2D([0], [0], color='black', lw=4,
label='Langmuir fit (this study)')]

# Creating colourscheme for scatterplots - should probably consider changing the colours...
# Creating colour scheme for scatterplots - should probably consider changing the colours...
c1_new = np.where(df1["Categories"] == 1, "blue", "orange")
c2_new = np.where(df1["Categories"] == 2, "lime", c1_new)
c3_new = np.where(df1["Categories"] == 4, "black", c2_new)
Expand Down Expand Up @@ -6857,7 +6856,7 @@ def mk_Andersen2021_figure_02(dpi=720, figsize=(7, 3), aspect=None, ms=None,
ax2.plot(x, y, color="black")
ax2.plot(x, y3, color="black", linestyle="dashed")
ax2.set_xlim(0.1, 10000)
# ax2.set_ylim(0, 20000) # Manually set the yaxis extents?
# ax2.set_ylim(0, 20000) # Manually set the y-axis extents?
ax2.set_xscale("log")
ax2.set_yscale("log")
ax2.tick_params(axis="both", labelsize=labelsize)
Expand All @@ -6876,3 +6875,112 @@ def mk_Andersen2021_figure_02(dpi=720, figsize=(7, 3), aspect=None, ms=None,
frameon=False)
# Save figure
AC.save_plot('ARNA_Andersen_figure_02', dpi=dpi, tight=tight)


def plt_seasonal_comparoisons_of_nitrate():
"""
Make a plot of seasonal nitrate at CVAO
"""
#
RunSet = 'ACID'
res = '4x5'
CoreRunsOnly = False
RunDict = get_dict_of_GEOSChem_model_output(res=res, RunSet=RunSet,
CoreRunsOnly=CoreRunsOnly,
folder4netCDF=True)
# Choose runs to use
d = {}
d['BC-BASE'] = '/users/ts551/scratch/GC/rundirs/P_ARNA//geosfp_4x5_standard.v12.9.0.BASE.2019.2020.ARNA.BCs.repeat//OutputDir/'
d['Acid-4x5-J00'] = '/users/ts551/scratch/GC/rundirs/P_ARNA//geosfp_4x5_aciduptake.v12.9.0.BASE.2019.2020.ARNA.DustUptake.JNIT.Isotherm.BCs.repeat.ON.II.diags.v2.J00//OutputDir/'
dates2use = [datetime.datetime(2019, 1+i, 1) for i in range(12)]
dsD = {}
for key in d.keys():
dsD[key] = AC.GetSpeciesConcDataset(wd=d[key],
dates2use=dates2use)

# Get model
from Prj_NITs_analysis import get_CVAO_NITs_data
# df_obs = get_CVAO_NITs_data()
# Update dates
dates = AC.dt64_2_dt( df_obs.index.values )
df_obs.index = [ update_year(i, year=2019) for i in df_obs.index ]

# Create a 'NIT-all'
prefix = 'SpeciesConc_'
for key in d.keys():
ds = dsD[key]
ds = AC.AddChemicalFamily2Dataset(ds, fam='NIT-all', prefix=prefix)
dsD[key] = ds

# Select for CVAO
from funcs4obs import gaw_2_loc
site = 'CVO'
lat, lon, alt, TZ = gaw_2_loc(site)
for key in d.keys():
ds = dsD[key]
ds = ds.sel(lat=lat, lon=lon, method='nearest')
ds = ds.sel(lev=ds.lev.values[0] )
dsD[key] = ds
# remove the redundent coordinates
for key in d.keys():
ds = dsD[key]
del ds['lat']
del ds['lon']
del ds['lev']
dsD[key] = ds

# plot up
import seaborn as sns
sns.set(color_codes=True)
var2plot = 'SpeciesConc_NIT-all'
fig, ax = plt.subplots()
colors = AC.get_CB_color_cycle()
for nKey, key in enumerate( d.keys() ):
ds = dsD[key]
AC.BASIC_seasonal_plot(dates=ds.time.values,
color=colors[nKey],
data=ds[var2plot].values*1E12,
ax=ax,
label=key)

# Add observations
AC.BASIC_seasonal_plot(dates=df_obs.index, color='k',
data=df_obs['Nitrate, pptv'].values, ax=ax,
label='Obs.')

plt.ylabel('All nitrate (inc. on dust)')
plt.legend()
plt.title('Seasonal cycle of ntirate at CVAO (pptv)')
AC.save_plot(title='NITs_comparison_CVAO_model_obs', dpi=720)
plt.close('all')

# rename obs
species = [i.split(prefix)[-1] for i in ds.data_vars]
for key in d.keys():
ds = dsD[key]
ds = ds.rename(name_dict=dict(zip(ds.data_vars, species)))
dsD[key] = ds

# plot up a stacked plot of nitrate by season.
key = 'Acid-4x5-J00'
ds = dsD[key]
del ds['ilev']
# vars2plot = [i for i in ds.data_vars if 'NIT' in i ]
vars2plot = ['NITD4', 'NITD3', 'NITD2', 'NITD1', 'NITs', 'NIT',]
df = ds[vars2plot+ ['NIT-all']].to_dataframe()
df = df * 1E12 # Conevrt to pptv
# for var in

plt.close('all')
fig, ax = plt.subplots()
df[vars2plot].plot.area()

# Add monthly mean for obs.
df_obs = df_obs.resample('1M').mean()
df_obs['Nitrate, pptv'].plot(label='obs', color='k', ls='--', zorder=99)
plt.autoscale(enable=True, axis='y')
plt.legend()

AC.save_plot(title='NITs_comparison_CVAO_by_species', dpi=720)
plt.close('all')

2 changes: 1 addition & 1 deletion arna/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -966,7 +966,7 @@ def get_local_folder(key, host=None, rtn_dict=False):
HEMCO_data = earth0_data_shelf + 'earth0_data/GEOS/ExtData/HEMCO/'
ARNA_data = '/users/ts551/scratch/data/ARNA/'
DataRoot = '/users/ts551/scratch/data/'
RunRoot = '/users/ts551/scratch/GC/rundirs/'
RunRoot = '/users/ts551/scratch/GC/rundirs/P_ARNA/'
elif ('earth0' in host):
NASA_data = '/work/data/NASA/'
ARNA_data = ''
Expand Down
Loading

0 comments on commit 3f2b9a3

Please sign in to comment.