Skip to content

GitLab

Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 16bf8a24 authored by Sakina Takache's avatar Sakina Takache
Browse files

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

parent a3840b49
Hide whitespace changes
Inline Side-by-side
Showing with 155 additions and 54 deletions
metrics/half_plume.py
View file @ 16bf8a24
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits.basemap 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?
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)
......@@ -29,6 +34,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
......@@ -36,6 +42,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:
......@@ -58,23 +65,21 @@ 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)
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()))
phi = ds['PHI'][hour, 1:, :, :] # geopotential
lon_deg, lat_deg = np.meshgrid(lon, lat) # degrees
lon, lat = to_radian(lon_deg), to_radian(lat_deg) # radians
......@@ -82,8 +87,6 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
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
......@@ -95,19 +98,17 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
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
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)
# 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
......@@ -115,7 +116,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)
......@@ -127,10 +127,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()
......@@ -143,13 +144,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
......@@ -158,21 +159,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) ))
......@@ -189,6 +185,7 @@ 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
......@@ -196,7 +193,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)]
......@@ -212,9 +210,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')
......@@ -227,27 +225,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():
# 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'
......@@ -258,19 +349,29 @@ def main():
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:
radius, _, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
plot_amount_emitted(read_DYN, read_CHI)
volume_half_mass_DYN = []
volume_half_mass_CHI = []
volume_DYN = []
volume_CHI = []
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)
_, _, 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)
main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment