Merge pull request #14 from E3SM-Project/hsiangheleellnl/chemdyg_o3hole

O3 hole plot x-axis fix

chemdyg/templates/chemdyg_index.bash

@@ -111,7 +111,7 @@ import pandas as pd
 pathout = './'
 
 indexfile = open(pathout+'index.html',"w")
-index = '<h> E3SM chem diagnostics package v1.1 (September 2022) </h>'
+index = '<h> E3SM chemistry diagnostics package (ChemDyg) </h>'
 index = index + '<pre> Test: ${short} </pre>'
 index = index + '<pre> Reference: Observations and Reanalysis </pre>'
 index = index + '<pre>   </pre>'

chemdyg/templates/chemdyg_noaa_co_comparison.bash

@@ -184,38 +184,51 @@ for n in range(len(Sta)):
     slope_noaa, intercept, r_value, p_value, std_err = linregress(nmonth_noaa[mask_noaa], CO_noaa[mask_noaa])
     lin_noaa = nmonth_noaa*slope_noaa+intercept
     lin_noaa_xa = xr.DataArray(lin_noaa, coords=[time_range_noaa], dims=["time"])
+    diff_noaa = CO_noaa - lin_noaa_xa
 
     nmonth = np.arange(0,len(CO_sel_1D),1)
     mask = ~np.isnan(CO_sel_1D)
-    if (len(nmonth[mask]) == 0):
-        continue
-    slope_e3sm, intercept, r_value, p_value, std_err = linregress(nmonth[mask], CO_sel_1D[mask])
-    lin_e3sm = nmonth*slope_e3sm+intercept
-    lin_e3sm_xa = xr.DataArray(lin_e3sm, coords=[time_range_noaa], dims=["time"])
 
-    diff = CO_sel_1D - lin_e3sm_xa
-    diff_noaa = CO_noaa - lin_noaa_xa
-    # plotting
-    fig, (ax1,ax2) = plt.subplots(2, 1,figsize=(10, 5))
-    ax1.plot(time_range_noaa,CO_sel_1D,'k')
-    ax1.plot(time_range_noaa[mask],lin_e3sm_xa[mask],'k--')
-    ax1.plot(time_range_noaa,CO_noaa,'r')
-    ax1.plot(time_range_noaa[mask_noaa],lin_noaa_xa[mask_noaa],'r--')
-    #ax1.set_ylim(0,200)
-    ax1.set_title('Surface CO at '+Sta[n]+' (Lat '+str(lat_noaa[0].values)+', Lon '+str(lon_noaa[0].values)+')')
-    line1 = 'E3SM mean:'+str(np.round(CO_sel_1D[mask].mean(),2))
-    line2 = 'E3SM trend:'+str(np.round(slope_e3sm*12,2))+' ppb/yr'
-    line3 = 'NOAA mean:'+str(np.round(CO_noaa[mask_noaa].mean().values,2))
-    line4 = 'NOAA trend:'+str(np.round(slope_noaa*12,2))+' ppb/yr'
-    ax1.legend([line1,line2, line3,line4])
-    #ax1.grid(True)
-
-    ax2.plot(time_range_noaa[mask],diff[mask],'k')
-    ax2.plot(time_range_noaa[mask_noaa],diff_noaa[mask_noaa],'r')
-    ax2.set_xlabel('time')
-    #ax2.set_ylim(-60,60)
-    ax2.set_ylabel('Anomalies')
-    pylab.savefig(pathout+'NOAA_CO_'+Sta[n]+'.png', dpi=600)
+    if (len(nmonth[mask]) != 0):
+        slope_e3sm, intercept, r_value, p_value, std_err = linregress(nmonth[mask], CO_sel_1D[mask])
+        lin_e3sm = nmonth*slope_e3sm+intercept
+        lin_e3sm_xa = xr.DataArray(lin_e3sm, coords=[time_range_noaa], dims=["time"])
+        diff = CO_sel_1D - lin_e3sm_xa
+
+        # plotting
+        fig, (ax1,ax2) = plt.subplots(2, 1,figsize=(10, 5))
+        ax1.plot(CO_noaa['time'],CO_sel_1D,'k')
+        ax1.plot(CO_noaa['time'][mask],lin_e3sm_xa[mask],'k--')
+        ax1.plot(CO_noaa['time'],CO_noaa,'r')
+        ax1.plot(CO_noaa['time'][mask_noaa],lin_noaa_xa[mask_noaa],'r--')
+        ax1.set_title('Surface CO at '+Sta[n]+' (Lat '+str(lat_noaa[0].values)+', Lon '+str(lon_noaa[0].values)+')')
+        line1 = 'E3SM mean:'+str(np.round(CO_sel_1D[mask].mean(),2))
+        line2 = 'E3SM trend:'+str(np.round(slope_e3sm*12,2))+' ppb/yr'
+        line3 = 'NOAA mean:'+str(np.round(CO_noaa[mask_noaa].mean().values,2))
+        line4 = 'NOAA trend:'+str(np.round(slope_noaa*12,2))+' ppb/yr'
+        ax1.legend([line1,line2, line3,line4])
+
+        ax2.plot(CO_noaa['time'][mask],diff[mask],'k')
+        ax2.plot(CO_noaa['time'][mask_noaa],diff_noaa[mask_noaa],'r')
+        ax2.set_xlabel('time')
+        #ax2.set_ylim(-60,60)
+        ax2.set_ylabel('Anomalies')
+        pylab.savefig(pathout+'NOAA_CO_'+Sta[n]+'.png', dpi=600)
+    else:
+        # plotting
+        fig, (ax1,ax2) = plt.subplots(2, 1,figsize=(10, 5))
+        ax1.plot(CO_noaa['time'],CO_noaa,'r')
+        ax1.plot(CO_noaa['time'][mask_noaa],lin_noaa_xa[mask_noaa],'r--')
+        ax1.set_title('Surface CO at '+Sta[n]+' (Lat '+str(lat_noaa[0].values)+', Lon '+str(lon_noaa[0].values)+')')
+        line3 = 'NOAA mean:'+str(np.round(CO_noaa[mask_noaa].mean().values,2))
+        line4 = 'NOAA trend:'+str(np.round(slope_noaa*12,2))+' ppb/yr'
+        ax1.legend([line3,line4])
+        print('E3SM data is not avaiable within this time period')
+        ax2.plot(CO_noaa['time'][mask_noaa],diff_noaa[mask_noaa],'r')
+        ax2.set_xlabel('time')
+        #ax2.set_ylim(-60,60)
+        ax2.set_ylabel('Anomalies')
+        pylab.savefig(pathout+'NOAA_CO_'+Sta[n]+'.png', dpi=600)
 
 EOF
 

chemdyg/templates/chemdyg_o3_hole_diags.bash

@@ -177,25 +177,30 @@ TOZ_array = TOZ_min_time.values.reshape((years,184))
 TOZ_mean = TOZ_array.mean(axis=0)
 TOZ_std = TOZ_array.std(axis=0)
 
-fig = plt.figure(figsize=(10,5))
-plt.plot(time_range[181:365],oz_avg, label ='Obs.')
-plt.fill_between(time_range[181:365],oz_avg+oz_std,oz_avg-oz_std, alpha=.5, linewidth=0)
-plt.plot(time_range[181:365],TOZ_mean, label ='E3SM')
-plt.fill_between(time_range[181:365],TOZ_mean+TOZ_std,TOZ_mean-TOZ_std, alpha=.5, linewidth=0)
+npdate = np.array(time_range[181:365])
+fig, ax = plt.subplots(figsize=(10, 5))
+plt.plot(npdate,oz_avg, label ='Obs.')
+plt.fill_between(npdate,oz_avg+oz_std,oz_avg-oz_std, alpha=.5, linewidth=0)
+plt.plot(npdate,TOZ_mean, label ='E3SM')
+plt.fill_between(npdate,TOZ_mean+TOZ_std,TOZ_mean-TOZ_std, alpha=.5, linewidth=0)
+date_form = DateFormatter("%b-%d")
+ax.xaxis.set_major_formatter(date_form)
 plt.legend(loc = 'upper left')
-plt.title('SH minimum total O3')
-plt.xlabel('Time')
+plt.title('SH minimum total O3 ('+str(startyear)+' - '+str(endyear)+')')
+plt.xlabel('Date')
 plt.ylabel('O3 conc. (DU)',fontsize='large')
 pylab.savefig(pathout+'to3mins_O3hole.png', dpi=600)
 
-fig = plt.figure(figsize=(10,5))
-plt.plot(time_range[181:365],area_avg, label ='Obs.')
-plt.fill_between(time_range[181:365],area_avg+area_std,area_avg-area_std, alpha=.5, linewidth=0)
-plt.plot(time_range[181:365],O3_mean, label ='E3SM')
-plt.fill_between(time_range[181:365],O3_mean+O3_std,O3_mean-O3_std, alpha=.5, linewidth=0)
+fig, ax = plt.subplots(figsize=(10, 5))
+plt.plot(npdate,area_avg, label ='Obs.')
+plt.fill_between(npdate,area_avg+area_std,area_avg-area_std, alpha=.5, linewidth=0)
+plt.plot(npdate,O3_mean, label ='E3SM')
+plt.fill_between(npdate,O3_mean+O3_std,O3_mean-O3_std, alpha=.5, linewidth=0)
+date_form = DateFormatter("%b-%d")
+ax.xaxis.set_major_formatter(date_form)
 plt.legend(loc = 'upper left')
-plt.title('O3 hole area (million of km2)')
-plt.xlabel('Time')
+plt.title('O3 hole area (million of km2) ('+str(startyear)+' - '+str(endyear)+')')
+plt.xlabel('Date')
 plt.ylabel('Area (million of km2)',fontsize='large')
 pylab.savefig(pathout+'to3areas_O3hole.png', dpi=600)
 

chemdyg/templates/chemdyg_summary_table.bash

@@ -218,8 +218,6 @@ for var in range(len(varname)):
    # write out annual chem tendency
             line_ann = line_ann + '<pre> '+ format('O3 burden (Tg)','28s')
             line_ann = line_ann + '     '+"{0:+.3e}".format(np.array(TOZ_total)) +'</pre>'
-            line_ann = line_ann + '<pre> '+ format('O3 emission (Tg/yr)','28s')
-            line_ann = line_ann + '     '+"{0:+.3e}".format(np.array(TDS_total)*unit_covet) +'</pre>'
             line_ann = line_ann + '<pre> '+ format('O3 deposition (Tg/yr)','28s')
             line_ann = line_ann + '     '+"{0:+.3e}".format(np.array(TDD_total)*unit_covet) +'</pre>'
             line_ann = line_ann + '<pre> '+ format('O3 production (Tg/yr)','28s')

chemdyg/templates/chemdyg_surf_o3_diags.bash

@@ -275,16 +275,18 @@ for i in range(len(O3EU_DJF)):
     O3_SEU_DJF[i] = (O3EU_DJF_SEU[i,:,:]*area_SEU[:,:]).sum()/area_SEU[:,:].sum()
 
 ndays = int(len(O3EU_JJA)/24)
-O3_ENA_JJA_24h = O3_ENA_JJA.reshape((ndays,24)).mean(axis=0)
-O3_WNA_JJA_24h = O3_WNA_JJA.reshape((ndays,24)).mean(axis=0)
-O3_NEU_JJA_24h = O3_NEU_JJA.reshape((ndays,24)).mean(axis=0)
-O3_SEU_JJA_24h = O3_SEU_JJA.reshape((ndays,24)).mean(axis=0)
+n_d = int(ndays*24)
+O3_ENA_JJA_24h = O3_ENA_JJA[0:n_d].reshape((ndays,24)).mean(axis=0)
+O3_WNA_JJA_24h = O3_WNA_JJA[0:n_d].reshape((ndays,24)).mean(axis=0)
+O3_NEU_JJA_24h = O3_NEU_JJA[0:n_d].reshape((ndays,24)).mean(axis=0)
+O3_SEU_JJA_24h = O3_SEU_JJA[0:n_d].reshape((ndays,24)).mean(axis=0)
 
 ndays = int(len(O3EU_DJF)/24)
-O3_ENA_DJF_24h = O3_ENA_DJF.reshape((ndays,24)).mean(axis=0)
-O3_WNA_DJF_24h = O3_WNA_DJF.reshape((ndays,24)).mean(axis=0)
-O3_NEU_DJF_24h = O3_NEU_DJF.reshape((ndays,24)).mean(axis=0)
-O3_SEU_DJF_24h = O3_SEU_DJF.reshape((ndays,24)).mean(axis=0)
+n_d = int(ndays*24)
+O3_ENA_DJF_24h = O3_ENA_DJF[0:n_d].reshape((ndays,24)).mean(axis=0)
+O3_WNA_DJF_24h = O3_WNA_DJF[0:n_d].reshape((ndays,24)).mean(axis=0)
+O3_NEU_DJF_24h = O3_NEU_DJF[0:n_d].reshape((ndays,24)).mean(axis=0)
+O3_SEU_DJF_24h = O3_SEU_DJF[0:n_d].reshape((ndays,24)).mean(axis=0)
 
 
 # ## read observations
@@ -614,8 +616,8 @@ for d in range(ndays):
 
 O3US_MDA8 = O3US_navg.max(axis=1)
 O3EU_MDA8 = O3EU_navg.max(axis=1)
-O3EU_MDA8_xr = xr.DataArray(O3EU_MDA8, coords=[day_timearray['time'],lat_EU,lon_EU], dims=["time","lat","lon"])
-O3US_MDA8_xr = xr.DataArray(O3US_MDA8, coords=[day_timearray['time'],lat_US,lon_US], dims=["time","lat","lon"])
+O3EU_MDA8_xr = xr.DataArray(O3EU_MDA8, coords=[day_timearray['time'][0:ndays],lat_EU,lon_EU], dims=["time","lat","lon"])
+O3US_MDA8_xr = xr.DataArray(O3US_MDA8, coords=[day_timearray['time'][0:ndays],lat_US,lon_US], dims=["time","lat","lon"])
 
 # US subregion: Western (WNA) and eastern (ENA) North America is split by 264°E
 lon1 = 234
@@ -636,16 +638,16 @@ O3_MDA8_ENA_1D = np.zeros(ndays)
 O3_MDA8_WNA_1D = np.zeros(ndays)
 O3_MDA8_NEU_1D = np.zeros(ndays)
 O3_MDA8_SEU_1D = np.zeros(ndays)
-for i in range(len(day_timearray)):
+for i in range(ndays):
     O3_MDA8_ENA_1D[i] = (O3_MDA8_ENA[i,:,:]*area_ENA[:,:]).sum()/area_ENA[:,:].sum()
     O3_MDA8_WNA_1D[i] = (O3_MDA8_WNA[i,:,:]*area_WNA[:,:]).sum()/area_WNA[:,:].sum()
     O3_MDA8_NEU_1D[i] = (O3_MDA8_NEU[i,:,:]*area_NEU[:,:]).sum()/area_NEU[:,:].sum()
     O3_MDA8_SEU_1D[i] = (O3_MDA8_SEU[i,:,:]*area_SEU[:,:]).sum()/area_SEU[:,:].sum()
 
-O3_MDA8_ENA_1D = xr.DataArray(O3_MDA8_ENA_1D, coords=[day_timearray['time']], dims=["time"])
-O3_MDA8_WNA_1D = xr.DataArray(O3_MDA8_WNA_1D, coords=[day_timearray['time']], dims=["time"])
-O3_MDA8_NEU_1D = xr.DataArray(O3_MDA8_NEU_1D, coords=[day_timearray['time']], dims=["time"])
-O3_MDA8_SEU_1D = xr.DataArray(O3_MDA8_SEU_1D, coords=[day_timearray['time']], dims=["time"])
+O3_MDA8_ENA_1D = xr.DataArray(O3_MDA8_ENA_1D, coords=[day_timearray['time'][0:ndays]], dims=["time"])
+O3_MDA8_WNA_1D = xr.DataArray(O3_MDA8_WNA_1D, coords=[day_timearray['time'][0:ndays]], dims=["time"])
+O3_MDA8_NEU_1D = xr.DataArray(O3_MDA8_NEU_1D, coords=[day_timearray['time'][0:ndays]], dims=["time"])
+O3_MDA8_SEU_1D = xr.DataArray(O3_MDA8_SEU_1D, coords=[day_timearray['time'][0:ndays]], dims=["time"])
 
 # calculate monthly mean
 O3_MDA8_ENA_month = np.zeros(12)