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

Commit c3c4743d authored by Thomas Dubos's avatar Thomas Dubos
Browse files

Cleanup half_plume.py

parent f1c5123a
......@@ -226,6 +226,33 @@ 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 concentration_map(read, level, day, filename):
lon, lat, area, volume, density = read(day)
# density = density[level, :, :]
density = (volume*density).sum(axis=0)
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 = np.linspace(0, m2.urcrnrx, density.shape[1])
y = np.linspace(0, m2.urcrnry, density.shape[0])
xx, yy = np.meshgrid(x, y)
m2.pcolormesh(xx, yy, density, zorder=2, alpha=0.5)
plt.colorbar(label='Sulfur dioxide density '+r'$(kg/m^3)$')
plt.savefig(filename%level, dpi=300)
plt.close()
def plot_maps(read_DYN, read_CHI, level, day):
concentration_map(read_DYN, level, day, 'metric_figures/DYN_plume%d_plots.png')
concentration_map(read_CHI, level, day, 'metric_figures/CHI_plume%d_plots.png')
def plot_amount_emitted(read_DYN, read_CHI, ndays=8, hour_step=3):
n=24*ndays//hour_step
......@@ -233,7 +260,7 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=8, hour_step=3):
DYN, CHI = [],[]
for i in range(n):
lon_DYN, lat_DYN, area_DYN, volume_DYN, density_DYN = read_DYN(3, hours[i]) # we shift the DYNAMICO data by 3 days
lon_DYN, lat_DYN, area_DYN, volume_DYN, density_DYN = read_DYN(0, hours[i])
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))
......@@ -257,7 +284,7 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=8, hour_step=3):
plt.title('Total amount of SO2 (should be kg)')
plt.savefig('amount_emitted.png')
plt.close()
# -------------------------------
# Trajectory of max column figure
......@@ -286,49 +313,6 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=8, hour_step=3):
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
......@@ -342,16 +326,14 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=8, hour_step=3):
plt.close()
def half_volume_plot(half_volume_DYN, half_volume_CHI):
def half_volume_plot(days, 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.plot(days, half_volume_DYN, marker='h', color='r')
plt.plot(days, 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.ylim(top = 0.2)
plt.yscale('log')
plt.legend(('DYNAMICO','CHIMERE'), loc='upper left')
plt.tight_layout()
......@@ -360,7 +342,6 @@ def half_volume_plot(half_volume_DYN, half_volume_CHI):
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')
......@@ -370,46 +351,63 @@ def volume_vs_time_plot(volume_DYN, volume_CHI):
plt.savefig('metric_figures/volume_vs_time.png')
plt.close()
def main_half_mass(read_DYN, read_CHI):
volume_half_mass_DYN = []
volume_half_mass_CHI = []
days = range(1,7)
for day in days:
_, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYN(day)
check_area(area_DYNAMICO, volume_DYNAMICO, radius, 'DYNAMICO')
_, _, area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHI(day)
check_area(area_CHIMERE, volume_CHIMERE, radius, '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'))
volume_tot = volume_DYNAMICO.sum()
volume_CHI = volume_half_mass_CHI/volume_tot
volume_DYN = volume_half_mass_DYN/volume_tot
half_volume_plot(days, volume_DYN, volume_CHI)
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("--level", type=int, default=75, help='Level for concentration maps')
getargs.add_argument("--day", type=int, default=7, help='Day for concentration maps')
getargs.add_argument("--maps", action='store_true', help='Concentration maps')
getargs.add_argument("--trajectory", action='store_true', help='Plume trajectory')
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')
getargs.add_argument("--all", action='store_true', help='Plot all')
args = getargs.parse_args()
if args.amount_emitted:
def read_DYN(day,hour):
return read_DYNAMICO(dynamico_nc, dynamico_hybrids, day, hour)
def read_NAT(day,hour):
return read_DYNAMICO_native(native_nc, native_area, dynamico_hybrids, day, hour)
def read_CHI(day,hour):
return read_CHIMERE(wrf_nc, chimere_nc, day, hour)
if args.native:
plot_amount_emitted(read_NAT, read_CHI)
else:
plot_amount_emitted(read_DYN, read_CHI)
def read_CHI(day,hour=0):
return read_CHIMERE(wrf_nc, chimere_nc, day, hour)
# DYNAMICO starts on June 1st
# CHIMERE starts on June 4th
# => we shift the DYNAMICO data by 3 days
if args.native:
def read_DYN(day,hour=0):
return read_DYNAMICO_native(native_nc, native_area, dynamico_hybrids, day+3, hour)
else:
if args.native:
_, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO_native(native_nc, native_area, dynamico_hybrids)
else:
_, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
def read_DYN(day,hour=0):
return read_DYNAMICO(dynamico_nc, dynamico_hybrids, day+3, hour)
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')
if args.all or args.amount_emitted:
plot_amount_emitted(read_DYN, read_CHI)
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'))
if args.all or args.maps:
plot_maps(read_DYN, read_CHI, args.level, args.day)
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
if args.all or args.trajectory:
plot_trajectory(read_DYN, read_CHI)
#max_volume_plot('GTRC1', 'q', CHIMERE_days)
#plots(metrics_CHIMERE, metrics_DYNAMICO)
if args.all or args.half_mass_volume:
main_half_mass(read_DYN, read_CHI)
if __name__ == "__main__":
main()
name: multiscale_transport
channels:
- defaults
dependencies:
- basemap-data-hires
- netcdf4
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