Version of half_plume.py used to generate the figures currently in the Thesis TeX file

metrics/half_plume.py

@@ -2,6 +2,8 @@ import numpy as np
 import netCDF4 as nc
 import matplotlib.pyplot as plt
 import argparse
+from mpl_toolkits.basemap import Basemap
+import pickle as pkl
 
 CHIMERE_days = ['04', '05', '06', '07', '08', '09', '10', '11']
 radius, g, p0 = 6.371e6, 9.81, 1e5 # Earth radius, gravity, reference pressure for hybrid coordinate
@@ -13,19 +15,21 @@ native_area      = '/data/WORKSPACE/dubos/Transport/NetCDF/Ai.nc'
 wrf_nc           = '/home/smailler/PUY60/geog_USGS_PUY60.nc'
 chimere_nc       = '/data/PLT8/pub/smailler/CHIMOUT-PUY60-tier2/out.201106%s00_24_PUY60-tier2.nc'
 
+
 def compute_pressure(ap, bp, p0, ps):
     lev, lon, lat = len(ap), ps.shape[1], ps.shape[0]
     p = np.zeros((lev, lat, lon)) # lat changes with row, lon with column
-    for i in range (0, lev): # 0 is ground level and not corresponding to a value of ilev?
+    for i in range (0, lev):
         p[i, :, :] = p0*ap[i] + ps*bp[i] # should return a pressure matrix(lat, lon) at each level
     return p
 
+
 def check_compute_pressure(ap, bp, p0, ps, p, rel_tol=1e-6):
     pp = compute_pressure(ap, bp, p0, ps)
     dp = np.absolute(pp-p)#.abs()
-#    print("Checking formula for hybrid pressure coordinate : max(p)=%f, max(pp)=%f, max(dp)%f" % (p.max(), pp.max(), dp.max()) )
     return dp.max()>p0*rel_tol  # true if there is a discrepancy
 
+
 def compute_volume(cell_area, hlay):
     lev = hlay.shape[0]
     volume = np.zeros(hlay.shape)
@@ -38,6 +42,7 @@ def compute_volume(cell_area, hlay):
         volume[l,:,:] = vol
     return volume
 
+
 def compute_mass(cell_area, p, g):
     ilev, nlat, nlon = p.shape
     mass = np.zeros((ilev-1, nlat, nlon))  # array of air mass per cell
@@ -45,6 +50,7 @@ def compute_mass(cell_area, p, g):
         mass[i, :, :] = cell_area * (p[i, :, :] - p[i+1, :, :]) / g
     return mass
 
+
 def max_volume_plot(variable_chim, variable_dyn, day_list):
     GTRC1_array = []
     for day in day_list:
@@ -67,6 +73,7 @@ def max_volume_plot(variable_chim, variable_dyn, day_list):
     plt.legend()#([q_for_plot, GTRC1_array],['q', 'gtrc'])
     plt.savefig('metric_figures/max_concentrations_versus_time.png')
 
+
 def to_radian(angle):
     return angle*(np.pi/180)
 
@@ -96,8 +103,6 @@ def read_DYNAMICO_common(hybrids, hour, area, ps, pmid, phi, q):
 def read_DYNAMICO(dynamico_nc, hybrids, day=7, hour=22):
     hour = day*24+hour # assume output every hour, first day is day 0  
 
-#    print('Reading DYNAMICO data from %s ...'%dynamico_nc)
-
     ds = nc.Dataset(dynamico_nc)
     lon, lat = ds['lon'][:], ds['lat'][:]
     nlon, nlat = len(lon), len(lat)
@@ -136,17 +141,16 @@ def read_DYNAMICO_native(dynamico_nc, area_nc, hybrids, day=7, hour=22):
 #    print(area.shape)
 
     lon, lat, area, ps, pmid, phi, q = map(extend, (lon, lat, area, ps, pmid, phi, q))
-
     volume, SO2_density = read_DYNAMICO_common(hybrids, hour, area, ps, pmid, phi, q)
     return lon, lat, area, volume, SO2_density
 
+
 def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
 
     hour = day*24+hour # assume output every hour
     day, hour = hour//24, hour%24
 
-    # NOTE :we work with molecule number rather than mass (constant factor = molar mass)
-#    print('Reading WRF data from %s ...'%wrf_nc)
+    # NOTE: we work with molecule number rather than mass (constant factor = molar mass)
     wrf = nc.Dataset(wrf_nc)
     cell_area = wrf.DX*wrf.DY / wrf['MAPFAC_MX'][0,:,:]/ wrf['MAPFAC_MY'][0,:,:] # m^2
 
@@ -154,7 +158,6 @@ def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
         return ds[var][hour,:,:,:]
 
     filename = chimere_nc%CHIMERE_days[day]
-#    print('Reading CHIMERE data from %s ...'%filename)
     ds = nc.Dataset(filename)
     hlay        = getvar('hlay')    # altitude above ground of upper interface (m)
     GTRC1       = getvar('GTRC1')   # SO2 molecular mixing ratio (ppb)
@@ -166,10 +169,11 @@ def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
 
     air_density = air_density*1e6                     # air molecule number per unit volume (m^3)
     SO2_density = 1e-9*GTRC1*air_density              # SO2 molecule number per unit volume (m^3)
-    SO2_density = SO2_density*SO2_molar_mass/avogadro # SO2 mass (kg) per unit volume (m^3)
+    SO2_density = SO2_density*SO2_molar_mass/avogadro # SO2 mass (kg) per unit volume (m^3) dimensions: lev=0-98, lat/y=0-398, lon/x=0-398
     
     return ds['lon'][:,:], ds['lat'][:,:], cell_area, cell_volume, SO2_density
 
+
 def volume_half_mass(volume, density, label):
     volume, density  = volume.flatten(), density.flatten()
     total_volume = volume.sum()
@@ -182,13 +186,13 @@ def volume_half_mass(volume, density, label):
     half_amount   = 0.5*cum_amount[-1]
 
     j = bisection_method(half_amount, cum_amount)
-    volume = sorted_volume[j:-1].sum()
+    half_mass_volume = sorted_volume[j:-1].sum()
     print('(%s) half-mass        = %f '%(label, half_amount))
-    print('(%s) half-mass volume = %f '%(label, volume/total_volume))
+    print('(%s) half-mass volume = %f '%(label, half_mass_volume/total_volume))
+    return half_mass_volume # half_mass_volume/total_volume DIVIDE BY DYNAMICO VOLUME
 
-def bisection_method(half_mass, cum_mass):
 
-#    print('Dichotomy for half-mass index ... ')
+def bisection_method(half_mass, cum_mass):
 
     N = len(cum_mass)
     a = 0   # lower_bound
@@ -197,21 +201,16 @@ def bisection_method(half_mass, cum_mass):
 
     if (cum_mass[m] <= half_mass) and (half_mass <= cum_mass[m+1]):
         return m
-
     else:
         while (cum_mass[m] > half_mass) or (half_mass > cum_mass[m+1]):
-#            print('a, b, m = ', a, b, m)
-#            print('m, m+1, half_mass = ', cum_mass[m], cum_mass[m+1], half_mass)
             if (cum_mass[m]>half_mass):
                 b = m # new upper bound
             else:
                 a = m # new lower bound
             m = int((a+b)/2)
-
-#        print('a, b, m = ', a, b, m)
-#        print('m, m+1, half_mass = ', cum_mass[m], cum_mass[m+1], half_mass)
         return m    
 
+
 def check_area(area, volume, radius, label):
     area = area.sum()
     print('Total %s domain area= %f = %d percent of Earth surface'%(label, area, 100*area/(4*np.pi*radius**2) )) 
@@ -227,7 +226,8 @@ def amount(lon, lat, area, vol, dens):
     index = np.unravel_index(np.argmax(mass_per_column), mass_per_column.shape)
     return total, lon[index], lat[index]
 
-def plot_amount_emitted(read_DYN, read_CHI, ndays=4, hour_step=3):
+
+def plot_amount_emitted(read_DYN, read_CHI, ndays=8, hour_step=3):
     n=24*ndays//hour_step
     hours = [ i*hour_step for i in range(n)]
 
@@ -243,9 +243,9 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=4, hour_step=3):
 
     amount_DYN, lon_DYN, lat_DYN = zip(*DYN)
     amount_CHI, lon_CHI, lat_CHI = zip(*CHI)
-    
-#    print(amount_DYN)
-#    print(amount_CHI)
+
+    # --------------------------------
+    # Amount of SO2 versus time figure
 
     plt.figure(); 
     plt.plot(hours, amount_DYN, label='DYNAMICO')
@@ -258,26 +258,120 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=4, hour_step=3):
     plt.savefig('amount_emitted.png')
     plt.close()
         
-    plt.figure(); 
-    plt.plot(hours, lon_DYN, label='DYNAMICO')
-    plt.plot(hours, lon_CHI, label='CHIMERE')    
-    plt.xlabel('hour')
-    plt.title('Longitude of maximum column-integrated SO2')
+    # -------------------------------
+    # Trajectory of max column figure
+
+    m = Basemap(projection='spstere',boundinglat=-10,lon_0=90,resolution='i')
+    m.drawcoastlines()
+    m.fillcontinents(color='yellowgreen',lake_color='dodgerblue')
+    # draw parallels and meridians.
+    m.drawparallels(np.arange(-80.,81.,20.))
+    m.drawmeridians(np.arange(-180.,181.,20.))
+    m.drawmapboundary(fill_color='lightseagreen')
+
+    lon_DYN_plot = lon_DYN[7:-1]
+    lat_DYN_plot = lat_DYN[7:-1]
+    lon_CHI_plot = lon_CHI[7:-1]
+    lat_CHI_plot = lat_CHI[7:-1]
+
+    x_DYN, y_DYN = m(lon_DYN_plot, lat_DYN_plot)
+    x_CHI, y_CHI = m(lon_CHI_plot, lat_CHI_plot)
+
+    #m.scatter(x_DYN, y_DYN, marker = 'h', facecolors='none', edgecolors='r', zorder=3, label = 'DYNAMICO')
+    plt.plot(x_DYN, y_DYN, marker = 'h', markerfacecolor='none', markeredgecolor='r', markersize = 2, linewidth=0.5, color = 'r', label = 'DYNAMICO')
+    #m.scatter(x_CHI, y_CHI, marker = 's', facecolors='none', edgecolors='b', zorder=2, label = 'CHIMERE')
+    plt.plot(x_CHI, y_CHI, marker = 's', markerfacecolor='none', markeredgecolor='b', markersize = 2, linewidth=0.5, color = 'b', label = 'CHIMERE')
+
+    plt.legend(fontsize = 'xx-small')
+    plt.savefig('metric_figures/lonlat_plots.png', dpi=300)
+    plt.close()
+
+    # -------------------------
+    # Plume at level 75 figures
+
+    density_CHI_lat_lon = density_CHI[75, :, :] # lev, lat/y, lon/x
+    density_DYN_lat_lon = density_DYN[75, :, :]
+
+    plt.figure()
+    m2 = Basemap(projection='spstere',boundinglat=-10,lon_0=90,resolution='h')
+    m2.drawcoastlines()
+    m2.fillcontinents(color='yellow',lake_color='dodgerblue')
+    # draw parallels and meridians.
+    m2.drawparallels(np.arange(-80.,81.,20.))
+    m2.drawmeridians(np.arange(-180.,181.,20.))
+    m2.drawmapboundary(fill_color='azure')
+
+    x_CHI = np.linspace(0, m2.urcrnrx, density_CHI_lat_lon.shape[1])
+    y_CHI = np.linspace(0, m2.urcrnry, density_CHI_lat_lon.shape[0])
+    xx_CHI, yy_CHI = np.meshgrid(x_CHI, y_CHI)
+
+    m2.pcolormesh(xx_CHI, yy_CHI, density_CHI_lat_lon, zorder=2, alpha=0.5)
+
+    plt.colorbar(label='Sulfur dioxide density '+r'$(kg/m^3)$')
+    plt.savefig('metric_figures/CHI_plume75_plots.png', dpi=300)
+
+
+    plt.figure()
+    m3 = Basemap(projection='spstere',boundinglat=-10,lon_0=90,resolution='h')
+    m3.drawcoastlines()
+    m3.fillcontinents(color='yellow',lake_color='dodgerblue')
+    # draw parallels and meridians.
+    m3.drawparallels(np.arange(-80.,81.,20.))
+    m3.drawmeridians(np.arange(-180.,181.,20.))
+    m3.drawmapboundary(fill_color='azure')
+
+    x_DYN = np.linspace(0, m3.urcrnrx, density_DYN_lat_lon.shape[1])
+    y_DYN = np.linspace(0, m3.urcrnry, density_DYN_lat_lon.shape[0])
+    xx_DYN, yy_DYN = np.meshgrid(x_DYN, y_DYN)
+
+    m3.pcolormesh(xx_DYN, yy_DYN, density_DYN_lat_lon, zorder=2, alpha=0.5)
+
+    plt.colorbar(label='Sulfur dioxide density '+r'$(kg/m^3)$')
+    plt.savefig('metric_figures/DYN_plume75_plots.png', dpi=300)
+
+    # ----------------------------
+    # Mass of SO2 with time figure
+
+    x_axis = np.linspace(0, 24*8, num=(24*8)//3)
+    plt.figure()
+    plt.plot(x_axis, amount_DYN, marker='h', color='r', label = 'DYNAMICO')
+    plt.plot(x_axis, amount_CHI, marker='s', color='b', label = 'CHIMERE')
+    plt.yscale('log')
     plt.legend()
-    plt.savefig('max_lon.png')
+    plt.savefig('metric_figures/mass_vs_time.png')
     plt.close()
-        
-    plt.figure(); 
-    plt.plot(hours, lat_DYN, label='DYNAMICO')
-    plt.plot(hours, lat_CHI, label='CHIMERE')    
-    plt.xlabel('hour')
-    plt.title('Latitude of maximum column-integrated SO2')
+
+
+def half_volume_plot(half_volume_DYN, half_volume_CHI):
+
+    # Half-mass volume versus time figure
+    x_axis = [1, 2, 3, 4, 5, 6, 7]
+    plt.figure()
+    plt.plot(x_axis, half_volume_DYN, marker='h', color='r')
+    plt.plot(x_axis, half_volume_CHI, marker='s', color='b')
+    plt.ylabel(r'$\frac{V_{half ~ plume}}{V_{DYNAMICO ~ mesh}}$', labelpad=10, rotation=0, fontsize=16)
+    plt.xlabel('Days after eruption')
+    plt.ylim(top = 0.2)
+    plt.yscale('log')
+    plt.legend(('DYNAMICO','CHIMERE'), loc='upper left')
+    plt.tight_layout()
+    plt.savefig('metric_figures/half_volume_vs_time.png')
+    plt.close()
+
+
+def volume_vs_time_plot(volume_DYN, volume_CHI):
+
+    x_axis = np.linspace(0, 24*8, num=(24*8)//3) # taking 8 days starting from eruption
+    plt.figure()
+    plt.plot(x_axis, volume_DYN, marker='h', color='r', label = 'DYNAMICO')
+    plt.plot(x_axis, volume_CHI, marker='s', color='b', label = 'CHIMERE')
+    plt.yscale('log')
     plt.legend()
-    plt.savefig('max_lat.png')
+    plt.savefig('metric_figures/volume_vs_time.png')
     plt.close()
-        
-def main():
 
+      
+def main():
     getargs = argparse.ArgumentParser(description='Analyze 3D transport experiments.')
     getargs.add_argument("--native", action='store_true', help='Read DYNAMICO native data instead of lon-lat')
     getargs.add_argument("--amount-emitted", action='store_true', help='Total emitted SO2 vs time')
@@ -302,14 +396,19 @@ def main():
             _, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
 
         check_area(area_DYNAMICO, volume_DYNAMICO, radius, 'DYNAMICO')
-        _, _, area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHIMERE(wrf_nc, chimere_nc)    
+        _, _, area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHIMERE(wrf_nc, chimere_nc, day)    
         check_area(area_CHIMERE, volume_CHIMERE, radius, 'CHIMERE')
 
-        volume_half_mass_DYNAMICO = volume_half_mass(volume_DYNAMICO, density_DYNAMICO, 'DYNAMICO')
-        volume_half_mass_CHIMERE = volume_half_mass(volume_CHIMERE, density_CHIMERE, 'CHIMERE')
+        volume_half_mass_DYN.append(volume_half_mass(volume_DYNAMICO, density_DYNAMICO, 'DYNAMICO'))
+        volume_half_mass_CHI.append(volume_half_mass(volume_CHIMERE, density_CHIMERE, 'CHIMERE'))
 
-    #max_volume_plot('GTRC1', 'q', CHIMERE_days)
+    volume_CHI = volume_half_mass_CHI/(volume_DYNAMICO.sum())
+    volume_DYN = volume_half_mass_DYN/(volume_DYNAMICO.sum())
+             
+    half_volume_plot(volume_DYN, volume_CHI)
+#    volume_vs_time_plot(volume_DYN, volume_CHI) THESE ARE MESH VOLUMES NOT PLUME VOLUMES
 
+    #max_volume_plot('GTRC1', 'q', CHIMERE_days)
     #plots(metrics_CHIMERE, metrics_DYNAMICO)
 
 if __name__ == "__main__":