Compute lon-lat where column-integrated amount of SO2 is maximum

metrics/half_plume.py

@@ -68,32 +68,25 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
 #    print('Reading DYNAMICO data from %s ...'%dynamico_nc)
 
     ds = nc.Dataset(dynamico_nc)
-    lon, lat = to_radian(ds['lon'][:]), to_radian(ds['lat'][:])
+    lon, lat = ds['lon'][:], ds['lat'][:]
     nlon, nlat = len(lon), len(lat)
     q   = ds['q'][hour, 0, :, :, :]
     ps  = ds['ps'][hour, :, :]
     phi = ds['PHI'][hour, 1:, :, :] # geopotential'
 #    print('shape(phi) = ', phi.shape)
-#    print("Max of DYNAMICO a at hour=%d : %f"%(hour, q.max()))
+#    print("Max of DYNAMICO q at hour=%d : %f"%(hour, q.max()))
 
-    lon, lat = np.meshgrid(lon, lat)
-#    print('Min(cos(lat) = ', np.cos(lat).min() )
+    lon_deg, lat_deg = np.meshgrid(lon, lat)          # degrees
+    lon, lat = to_radian(lon_deg), to_radian(lat_deg) # radians
     area = (radius**2)*np.cos(lat)*(2*np.pi/nlon)*(np.pi/nlat)
 
     volume = compute_volume(area/g, phi) # m^3
-#    print('volume array length = ', len(volume), volume.size)
-#    print('volume array = ', volume)
-
 
 #    print('  geopotential last hour = ', phi[:,12,27])
 
     ds1 = nc.Dataset(dynamico_hybrids)
     ap, bp = ds1['hyai'][:], ds1['hybi'][:]    
-
-#    print('  ap shape, bp shape, ps shape = ',ap.shape, bp.shape, ps.shape)
-
     lev = len(ap) - 1
-
     p = np.zeros((lev, nlat, nlon)) # lat changes with row, lon with column
 
     if check_compute_pressure(ds1['hyam'][:], ds1['hybm'][:], p0, ps, ds['P'][hour,:,:,:] ):
@@ -106,7 +99,7 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
 
 #    print("Total SO2 mass in DYNAMICO = ", (q*air_mass).sum())
 #    print("Max SO2 mixing ratio (kg/kg) in DYNAMICO = ", q.max())
-    return radius, area, volume.flatten(), SO2_density.flatten()    
+    return radius, lon_deg, lat_deg, area, volume, SO2_density
 
 def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
 
@@ -134,12 +127,13 @@ 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 per unit volume (m^3)
+    SO2_density = SO2_density*SO2_molar_mass/avogadro # SO2 mass (kg) per unit volume (m^3)
     
-    return cell_area, cell_volume.flatten(), SO2_density.flatten() 
+    return ds['lon'][:,:], ds['lat'][:,:], cell_area, cell_volume, SO2_density
 
 def volume_half_mass(volume, density, label):
-    total_volume=volume.sum()
+    volume, density  = volume.flatten(), density.flatten()
+    total_volume = volume.sum()
     amount = volume*density
     # sort density per cell from small to large
     ind = np.argsort(density) # indices sorting density
@@ -195,21 +189,30 @@ getargs.add_argument("--amount-emitted", action='store_true', help='Total emitte
 getargs.add_argument("--half-mass-volume", action='store_true', help='Minimal volume containing half of emitted SO2 mass')
 args = getargs.parse_args()
 
+def amount(lon, lat, area, vol, dens):
+    mass = vol*dens # kg
+    total = mass.sum()  # kg
+    mass_per_column  = mass.sum(axis=0)/area # kg/m^2
+    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):
     n=24*ndays//hour_step
     hours = [ i*hour_step for i in range(n)]
 
-    amount_DYN=np.zeros(n)
-    amount_CHI=np.zeros(n)
-
-    def amount(vol, dens): return (vol*dens).sum()
-
+    DYN, CHI = [],[]
     for i in range(n):
-        radius, area_DYN, volume_DYN, density_DYN = read_DYN(3, hours[i]) # we shift the DYNAMICO data by 3 days
-        area_CHI,         volume_CHI, density_CHI = read_CHI(0, hours[i])
-        amount_DYN[i] = amount(volume_DYN, density_DYN)
-        amount_CHI[i] = amount(volume_CHI, density_CHI)
+        radius, lon_DYN, lat_DYN, area_DYN, volume_DYN, density_DYN = read_DYN(3, hours[i]) # we shift the DYNAMICO data by 3 days
+        lon_CHI, lat_CHI, area_CHI, volume_CHI, density_CHI = read_CHI(0, hours[i])
+        DYN.append(amount(lon_DYN, lat_DYN, area_DYN, volume_DYN, density_DYN)) # tuple (total, lon, lat)
+        CHI.append(amount(lon_CHI, lat_CHI, area_CHI, volume_CHI, density_CHI))
+
+        if DYN[i][0]>0.:
+            print('amount at %d'%i, DYN[i][0], CHI[i][0], CHI[i][0]/DYN[i][0])
 
+    amount_DYN, lon_DYN, lat_DYN = zip(*DYN)
+    amount_CHI, lon_CHI, lat_CHI = zip(*CHI)
+    
 #    print(amount_DYN)
 #    print(amount_CHI)
 
@@ -218,10 +221,29 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=4, hour_step=3):
     plt.plot(hours, amount_CHI, label='CHIMERE')    
     plt.yscale('log')
     plt.xlabel('hour')
-#    plt.legend()
+    plt.ylim(bottom=1e7)
+    plt.legend()
     plt.title('Total amount of SO2 (should be kg)')
     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')
+    plt.legend()
+    plt.savefig('max_lon.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')
+    plt.legend()
+    plt.savefig('max_lat.png')
+    plt.close()
         
 def main():
 #    dynamico_nc      = '/data/PLT8/pub/smailler/PUY60-dynamico/dynamico.nc'
@@ -239,9 +261,9 @@ def main():
     if args.amount_emitted:
         plot_amount_emitted(read_DYN, read_CHI)
     else:
-        radius, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
+        radius, _, _, 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)    
         check_area(area_CHIMERE, volume_CHIMERE, radius, 'CHIMERE')
 
         volume_half_mass_DYNAMICO = volume_half_mass(volume_DYNAMICO, density_DYNAMICO, 'DYNAMICO')