Module with pickles, without plots

metrics/half_plume_plus_pickles.py

+import numpy as np
+import netCDF4 as nc
+import matplotlib.pyplot as plt
+#import basemap
+import pickle as pkl
+
+
+CHIMERE_days = ['04', '05', '06', '07', '08', '09', '10', '11']
+
+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?
+        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)
+    for l in range(lev):
+        if l==0:
+            height = hlay[0,:,:]
+        else:
+            height = hlay[l,:,:]-hlay[l-1,:,:]
+        vol =  cell_area * height
+        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
+    for i in range(ilev-1):
+        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:
+        chimere_nc = '/data/PLT8/pub/smailler/CHIMOUT-PUY60-tier2/out.201106%s00_24_PUY60-tier2.nc' %day
+        dsf = nc.Dataset(chimere_nc)    
+        GTRC1_array.append(np.amax(dsf[variable_chim][:,:,:])) # each time is in a separate file
+
+    ds = nc.Dataset('/data/PLT8/pub/smailler/PUY60-dynamico/dynamico.nc')
+    q_array = ds[variable_dyn][:,0,:,:,:]
+    t_for_plot = []
+    q_for_plot = []
+    for t in range(0, q_array.shape[0], 48):
+        t_for_plot.append(t/48.)
+        q_for_plot.append(np.amax(ds[variable_dyn][t,0,:,:,:]))
+
+    plt.plot(t_for_plot, q_for_plot, '-r+')
+    plt.plot(day_list, GTRC1_array, '-b+')
+    plt.ylabel('Concentration')
+    plt.xlabel('t')
+    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)
+
+def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
+    hour = day*24+hour # assume output every hour, first day is day 0  
+    radius, g, p0 = 6.371e6, 9.81, 1e5 # Earth radius, gravity, reference pressure for hybrid coordinate
+
+#    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)
+    q   = ds['q'][hour, 0, :, :, :]
+    ps  = ds['ps'][hour, :, :]
+    phi = ds['PHI'][hour, 1:, :, :] # geopotential'
+#    print('shape(phi) = ', phi.shape)
+#    print("Max of DYNAMICO q at hour=%d : %f"%(hour, q.max()))
+
+    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('  geopotential last hour = ', phi[:,12,27])
+
+    ds1 = nc.Dataset(dynamico_hybrids)
+    ap, bp = ds1['hyai'][:], ds1['hybi'][:]    
+    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,:,:,:] ):
+        print("  Discrepancy in computed pressure. Stop.")
+        return True
+
+    p = compute_pressure(ap, bp, p0, ps)
+    air_mass = compute_mass(area, p, g)
+    SO2_density = q * air_mass / volume # (kg/kg) * (kg) / (m^3) = kg / m^3 dimensions: lev = 1-40, lat/y = -89 - 89, lon/x = 1 - 359
+
+#    print("Total SO2 mass in DYNAMICO = ", (q*air_mass).sum())
+#    print("Max SO2 mixing ratio (kg/kg) in DYNAMICO = ", q.max())
+    return radius, lon_deg, lat_deg, 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)
+    wrf = nc.Dataset(wrf_nc)
+    cell_area = wrf.DX*wrf.DY / wrf['MAPFAC_MX'][0,:,:]/ wrf['MAPFAC_MY'][0,:,:] # m^2
+
+    def getvar(var):
+        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)
+    air_density = getvar('airm')    # air molecule number per unit volume (cm^3)
+    cell_volume = compute_volume(cell_area, hlay)  # m^3
+
+    SO2_molar_mass = 64.066/1e3 # kg/mol
+    avogadro = 6.02214076e23    # Avogadro number
+
+    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) 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()
+    amount = volume*density
+    # sort density per cell from small to large
+    ind = np.argsort(density) # indices sorting density
+    sorted_amount = amount[ind]
+    sorted_volume = volume[ind]
+    cum_amount    = np.cumsum(amount)
+    half_amount   = 0.5*cum_amount[-1]
+
+    j = bisection_method(half_amount, cum_amount)
+    half_mass_volume = sorted_volume[j:-1].sum()
+    print('(%s) half-mass        = %f '%(label, half_amount))
+    print('(%s) half-mass volume = %f '%(label, half_mass_volume/total_volume))
+    return half_mass_volume/total_volume
+
+def bisection_method(half_mass, cum_mass):
+
+#    print('Dichotomy for half-mass index ... ')
+
+    N = len(cum_mass)
+    a = 0   # lower_bound
+    b = N-1 # upper_bound
+    m = int((a+b)/2)
+
+    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) )) 
+    volume = volume.sum()
+    print('Mean model top for %s is at %d meters'%(label, volume/area))
+
+# --------------------------------------------------------------
+
+import argparse
+import matplotlib.pyplot as plt
+
+getargs = argparse.ArgumentParser(description='Analyze 3D transport experiments.')
+getargs.add_argument("--amount-emitted", action='store_true', help='Total emitted SO2 vs time')
+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=8, hour_step=3):
+    n=24*ndays//hour_step
+    hours = [ i*hour_step for i in range(n)]
+
+    DYN, CHI = [],[]
+    for i in range(n):
+        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('lon_DYN = ', lon_DYN)
+    print('TYPE = ', type(lon_DYN))
+
+    with open('saved_amount_DYN.pkl', 'wb') as f:
+        pkl.dump(amount_DYN, f)
+    with open('saved_amount_CHI.pkl', 'wb') as f:
+        pkl.dump(amount_CHI, f)
+
+    with open('saved_density_DYN.pkl', 'wb') as f:
+        pkl.dump(density_DYN, f)
+    with open('saved_density_CHI.pkl', 'wb') as f:
+        pkl.dump(density_CHI, f)
+
+    with open('saved_lon_DYN.pkl', 'wb') as f:
+        pkl.dump(lon_DYN, f)
+    with open('saved_lat_DYN.pkl', 'wb') as f:
+        pkl.dump(lat_DYN, f)
+    with open('saved_lon_CHI.pkl', 'wb') as f:
+        pkl.dump(lon_CHI, f)
+    with open('saved_lat_CHI.pkl', 'wb') as f:
+        pkl.dump(lat_CHI, f)
+
+
+
+#    print(amount_DYN)
+#    print(amount_CHI)
+
+    plt.figure(); 
+    plt.plot(hours, amount_DYN, label='DYNAMICO')
+    plt.plot(hours, amount_CHI, label='CHIMERE')    
+    plt.yscale('log')
+    plt.xlabel('hour')
+    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'
+#    dynamico_hybrids = '/data/PLT8/pub/smailler/PUY60-dynamico/apbp.nc'
+    dynamico_nc      = '/data/PLT8/pub/smailler/multiscale/dynamico.nc'
+    dynamico_hybrids = '/data/PLT8/pub/smailler/multiscale/apbp.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 read_DYN(day,hour):
+        return read_DYNAMICO(dynamico_nc, dynamico_hybrids, day, hour)
+    def read_CHI(day,hour):
+        return read_CHIMERE(wrf_nc, chimere_nc, day, hour)
+
+    if args.amount_emitted:
+        plot_amount_emitted(read_DYN, read_CHI)
+    else:
+
+       volume_half_mass_DYNAMICO = []
+       volume_half_mass_CHIMERE = []
+
+       for day in (0, 1, 2, 3, 4, 5, 6):
+            radius, _, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids, day+3)
+            check_area(area_DYNAMICO, volume_DYNAMICO, radius, 'DYNAMICO')
+            _, _, 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.append(volume_half_mass(volume_DYNAMICO, density_DYNAMICO, 'DYNAMICO'))
+            volume_half_mass_CHIMERE.append(volume_half_mass(volume_CHIMERE, density_CHIMERE, 'CHIMERE'))
+            
+            with open('volume_half_mass_DYN.pkl', 'wb') as f:
+                pkl.dump(volume_half_mass_DYNAMICO, f)
+            with open('volume_half_mass_CHI.pkl', 'wb') as f:
+                pkl.dump(volume_half_mass_CHIMERE, f)
+
+
+    #max_volume_plot('GTRC1', 'q', CHIMERE_days)
+
+    #plots(metrics_CHIMERE, metrics_DYNAMICO)
+
+main()