diff --git a/blueprints/test_chemistry/VERSION b/blueprints/test_chemistry/VERSION index 2edeafb09db0093bae6ff060e2dcd2166f5c9387..8fdd954df9831dfd29ceec0d74829b02f3f5d8c3 100644 --- a/blueprints/test_chemistry/VERSION +++ b/blueprints/test_chemistry/VERSION @@ -1 +1 @@ -20 \ No newline at end of file +22 \ No newline at end of file diff --git a/blueprints/test_chemistry/WRF_DOMAIN_PLOT_COAWST.png b/blueprints/test_chemistry/WRF_DOMAIN_PLOT_COAWST.png new file mode 100644 index 0000000000000000000000000000000000000000..52825041e10eb5113f7b33b4bd3e027f73fb3c41 Binary files /dev/null and b/blueprints/test_chemistry/WRF_DOMAIN_PLOT_COAWST.png differ diff --git a/blueprints/test_chemistry/machine_specific/alcc_intel/config.bash b/blueprints/test_chemistry/config.bash similarity index 66% rename from blueprints/test_chemistry/machine_specific/alcc_intel/config.bash rename to blueprints/test_chemistry/config.bash index 0a90406b48b5d55a3b48c27433facc44c1b9a5fd..7dd12edb2fe58da685375b08069b418d9e3d16f4 100644 --- a/blueprints/test_chemistry/machine_specific/alcc_intel/config.bash +++ b/blueprints/test_chemistry/config.bash @@ -19,7 +19,7 @@ WRFDir=${WRF_SRC_PATH} # --- Input data settings --- # path to geogrid input data -geogDir=${WRF_GEOG_PATH} +geogDir=/alcc/gpfs2/home/mbees/data/geog/ # meteo input # Vtable for the chosen meteo input @@ -28,7 +28,7 @@ metVtableFile=${WPS_SRC_PATH}/ungrib/Variable_Tables/Vtable.GFS metInc=1 # full path to a met input file - you can use any "%<>" abbreviations known # to the "date" command -metFilePattern="${WRF_GFS_METEO_PATH}/GF%Y%m%d%H" +metFilePattern="/alcc/gpfs2/home/mbees/data/meteo/GFS/GF%Y%m%d%H" # example: # "/glade/p/rda/data/ds083.2/grib2/%Y/%Y.%m/fnl_%Y%m%d_%H_00.grib2" @@ -57,26 +57,6 @@ restartRootDir=${SCRATCH}/WRF/restart/ # remove run directory after run is finished? removeRunDir=false -# --- MPI settings --- - -mpiCommandPre="mpirun /usr/bin/time -v" -mpiCommandMain="mpirun /usr/bin/time -v" -mpiCommandReal=${mpiCommandPre} - -# --- Batch system --- - -# argument to submit a job in a held state -batchHoldArgument="--hold" -# command to release a held job -batchReleaseCommand="scontrol release" -# command to submit jobs to the queueing system -batchSubmitCommand=sbatch -# dependency argument for chaining runs upon submission -batchDepArgument="--dependency=afterany:__id__" -# sed command ("used as s/__command/\1/") to retrieve job run PID upon -# submission with $batchSubmitCommand -batchPidSedCommand="Submitted batch job \(.*\)" - # --- Chemistry --- withChemistry=true @@ -84,28 +64,28 @@ withChemistry=true # WRF-Chem installation directory WRFChemDir=${WRF_CHEM_SRC_PATH} -# megan_bio_emiss installation directory -WRFMEGANdir=${WRF_CHEM_MEGAN_BIO_EMISS_PATH} -# mozbc installation directory -WRFMOZARTdir=${WRF_CHEM_MOZBC_PATH} -# wesley/exocoldens installation directory -WRFmztoolsdir=${WRF_CHEM_WES_COLDENS_PATH} -# anthro_emiss installation directory -WRFanthrodir=${WRF_CHEM_ANTHRO_EMIS_PATH} -# fire_emis installation directory -WRFfiredir=${WRF_CHEM_FIRE_EMIS_PATH} +# path to utility executables +meganBioEmissBin=megan_bio_emiss +mozbcBin=mozbc +weselyBin=wesely +exo_coldensBin=exo_coldens +anthro_emisBin=anthro_emis +fire_emisBin=fire_emis + +# path to Wesely and Exo_Coldens input data +wesColdensDataPath=${WESELY_EXO_COLDENS_DATA} # path to MEGAN input data MEGANdir=/alcc/gpfs2/home/mbees/data/emissions/biogenic/MEGAN # use anthro_emiss or predefined files? -emissUseAnthroEmiss=true +emissUseAnthroEmiss=false # raw emission input - the files you read in with anthro_emiss -emissDir=/alcc/gpfs2/home/mbees/data/emissions/anthropogenic/EDGAR-HTAP/MOZART_MOSAIC +emissDir=/alcc/gpfs2/home/mbees/data/emissions/anthropogenic/EDGARv5/MOZART_MOSAIC # emission conversion script for anthro_emis - must match emissions in emissDir -emissInpFile=emis_edgarhtap_mozmos.inp +emissInpFile=emis_edgarv5_mozmos.inp # year the emissions are valid for (for offset calculation) -emissYear=2010 +emissYear=2015 # FINN fires fireFilePattern="/alcc/gpfs2/home/mbees/data/emissions/fires/FINN/GLOB_MOZ4_%Y%j.txt" diff --git a/blueprints/test_chemistry/crontab.operational_chemistry b/blueprints/test_chemistry/crontab.operational_chemistry deleted file mode 100644 index 6ef3dc350a80916af08a299fadb2f6b32860dfde..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/crontab.operational_chemistry +++ /dev/null @@ -1,8 +0,0 @@ -SHELL=/bin/bash -MAILTO=christoph.knote@med.uni-augsburg.de -USER=knotechr -# -# ### daily WRF-Chem forecasts ### -# -00 02 * * * /bin/bash -l /alcc/gpfs2/home/u/knotechr/wrfotron/master.bash /alcc/gpfs2/home/u/knotechr/wrfotron/blueprints/operational_chemistry `/bin/date -u --date="1 days ago" "+\%Y \%m \%d 00"` 96 03 -# diff --git a/blueprints/test_chemistry/emis_edgarhtap_mozmos.inp b/blueprints/test_chemistry/emis_edgarhtap_mozmos.inp index 1d9b89e7c4e63934108664816e0af4b5c7f31828..30b8de89cad5717b43a3988c72e3e2f3612dbd08 100644 --- a/blueprints/test_chemistry/emis_edgarhtap_mozmos.inp +++ b/blueprints/test_chemistry/emis_edgarhtap_mozmos.inp @@ -11,7 +11,7 @@ cat_var_prefix = ' ', serial_output = .true., start_output_time = '__startDate__', - stop_output_time = '__endDate__', + stop_output_time = '__startDate__', output_interval = 3600, data_yrs_offset = __emissYearOffset__, emissions_zdim_stag = 1, diff --git a/blueprints/test_chemistry/emis_edgarv5_mozmos.inp b/blueprints/test_chemistry/emis_edgarv5_mozmos.inp new file mode 100644 index 0000000000000000000000000000000000000000..5bd229eb59b94700be346daa82ac08be472d4ad3 --- /dev/null +++ b/blueprints/test_chemistry/emis_edgarv5_mozmos.inp @@ -0,0 +1,29 @@ +&CONTROL + domains = __domains__ + anthro_dir = '__emissDir__' + src_file_prefix = 'edgar_v5_2015_' + src_file_suffix = '_0.1x0.1.nc' + src_names = 'CO(28)','NOx(30)','SO2(64)','NH3(17)','BC(12)','OC(12)','PM2.5(1)','PM10(1)','BENZENE(78)','BIGALK(72)','BIGENE(56)','C2H2(26)','C2H4(28)', + 'C2H5OH(46)','C2H6(30)','C3H6(42)','C3H8(44)','CH2O(30)', + 'CH3CHO(44)','CH3COCH3(58)','CH3OH(32)','CH3COOH(60)','HCOOH(46)', + 'MEK(72)','TOLUENE(92)','XYLENES(106)' + sub_categories = 'TOTAL', + serial_output = .true., + start_output_time = '__startDate__', + stop_output_time = '__startDate__', + output_interval = 3600, + data_yrs_offset = __emissYearOffset__, + emissions_zdim_stag = 1, + emis_map = 'CO->CO','NO->0.8*NOx','NO2->0.2*NOx','SO2->SO2','NH3->NH3','BIGALK->BIGALK','BIGENE->BIGENE', + 'C2H4->C2H4','C2H5OH->C2H5OH','C2H6->C2H6','CH2O->CH2O','CH3CHO->CH3CHO', + 'CH3COCH3->CH3COCH3','CH3OH->CH3OH','MEK->MEK','TOLUENE->TOLUENE', + 'C3H6->C3H6','C3H8->C3H8','ISOP->0.0*CO','C10H16->0.0*CO', + 'SULF->0.0*SO2','C2H2->0.00561790*CO','BENZENE->BENZENE','XYLENE->XYLENES', + 'GLY->0.0*CO','MACR->0.0*CO','MGLY->0.0*CO','MVK->0.0*CO', + 'HCOOH->0.0*CO','HONO->0.0*CO','VOCA->0.04*CO','VOCBB->0.0*CO', + 'ECI(a)->0.1*BC','ECJ(a)->0.9*BC','ORGI(a)->0.1*OC','ORGJ(a)->0.9*OC','PM25I(a)->0.1*PM2.5', + 'PM25J(a)->0.9*PM2.5','PM_10(a)->PM10 + -1.0*PM2.5','SO4I(a)->0.0*PM10','SO4J(a)->0.0*PM10','NO3I(a)->0.0*PM10', + 'NO3J(a)->0.0*PM10','NH4I(a)->0.0*PM10','NH4J(a)->0.0*PM10','NAI(a)->0.0*PM10','NAJ(a)->0.0*PM10', + 'CLI(a)->0.0*PM10','CLJ(a)->0.0*PM10','CO_A->CO','CO_BB->0.0*CO','ORGI_A(a)->0.0*PM10', + 'ORGI_BB(a)->0.0*PM10','ORGJ_A(a)->0.0*PM10','ORGJ_BB(a)->0.0*PM10' +/ diff --git a/blueprints/test_chemistry/iofields.met b/blueprints/test_chemistry/iofields.met index 54f5fbf0547fcbd4e3c3c16abbf6a11f88b4aa8f..10455b4f9351e774714eb0564e92e27393fc1a06 100644 --- a/blueprints/test_chemistry/iofields.met +++ b/blueprints/test_chemistry/iofields.met @@ -1,3 +1,4 @@ +:h:0:alt +:h:0:gsw +:h:0:tke ++:h:0:SWDDIF \ No newline at end of file diff --git a/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/main b/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/main deleted file mode 100644 index 7cea60269ce19abca900b055770ae9dd65cf478d..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/main +++ /dev/null @@ -1,23 +0,0 @@ -#!/bin/bash -l -#SBATCH --partition=alcc1,epyc -#SBATCH -o __runDir__/__mainJobName__.%j.%N.out -#SBATCH -D __runDir__ -#SBATCH -J __mainJobName__ -#SBATCH --nodes=2 -#SBATCH --ntasks-per-node=28 -#SBATCH --mem-per-cpu=2000 -#SBATCH --mail-type=FAIL -#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=12:00:00 - -# -# variable $MACHINEFILE just holds the filename where acquired -# nodes/cores names are written. e.g. -# -MACHINEFILE=__runDir__/slurm.hosts - -# -# Generate Machinefile for openmpi such that hosts are in the same -# order as if run via srun -# -srun hostname -s | sort -n > $MACHINEFILE diff --git a/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/post b/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/post deleted file mode 100644 index e7ed998f3c7a3bcf8f63cde38d0d4614c2d07c49..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/post +++ /dev/null @@ -1,22 +0,0 @@ -#!/bin/bash -l -#SBATCH --partition=alcc1 -#SBATCH -o __runDir__/__postJobName__.%j.%N.out -#SBATCH -D __runDir__ -#SBATCH -J __postJobName__ -#SBATCH --ntasks=1 -#SBATCH --mem=5G -#SBATCH --mail-type=FAIL -#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=02:00:00 - -# -# variable $MACHINEFILE just holds the filename where acquired -# nodes/cores names are written. e.g. -# -MACHINEFILE=__runDir__/slurm.hosts - -# -# Generate Machinefile for openmpi such that hosts are in the same -# order as if run via srun -# -srun hostname -s | sort -n > $MACHINEFILE diff --git a/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/pre b/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/pre deleted file mode 100644 index a82c3dd0cc529e1bac55cf6a6de9f7de6196a4f8..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/pre +++ /dev/null @@ -1,22 +0,0 @@ -#!/bin/bash -l -#SBATCH --partition=alcc1 -#SBATCH -o __runDir__/__preJobName__.%j.%N.out -#SBATCH -D __runDir__ -#SBATCH -J __preJobName__ -#SBATCH --ntasks=1 -#SBATCH --mem=10G -#SBATCH --mail-type=FAIL -#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=03:00:00 - -# -# variable $MACHINEFILE just holds the filename where acquired -# nodes/cores names are written. e.g. -# -MACHINEFILE=__runDir__/slurm.hosts - -# -# Generate Machinefile for openmpi such that hosts are in the same -# order as if run via srun -# -srun hostname -s | sort -n > $MACHINEFILE diff --git a/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/staging b/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/staging deleted file mode 100644 index 7e33cdb7f332dd774d56b135bdf0489df2607733..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_gnu/batch_preambles/staging +++ /dev/null @@ -1,22 +0,0 @@ -#!/bin/bash -l -#SBATCH --partition=alcc1 -#SBATCH -o __runDir__/__stagingJobName__.%j.%N.out -#SBATCH -D __runDir__ -#SBATCH -J __stagingJobName__ -#SBATCH --ntasks=1 -#SBATCH --mem=5G -#SBATCH --mail-type=FAIL -#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=01:00:00 - -# -# variable $MACHINEFILE just holds the filename where acquired -# nodes/cores names are written. e.g. -# -MACHINEFILE=__runDir__/slurm.hosts - -# -# Generate Machinefile for openmpi such that hosts are in the same -# order as if run via srun -# -srun hostname -s | sort -n > $MACHINEFILE diff --git a/blueprints/test_chemistry/machine_specific/alcc_gnu/config.bash b/blueprints/test_chemistry/machine_specific/alcc_gnu/config.bash deleted file mode 100644 index 2cdec29962bcec421e8990f2bcae98349d19c91b..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_gnu/config.bash +++ /dev/null @@ -1,126 +0,0 @@ -#!/bin/bash -# ------------------------------------------------------------------------------ -# WRFOTRON v 2.0b -# Christoph Knote (LMU Munich, Germany) -# 06/2016 -# christoph.knote@lmu.de -# ------------------------------------------------------------------------------ - -# path to the WRFotron installation -chainDir=${HOME}/wrfotron - -# --- Executable locations --- - -# WPS installation directory -WPSDir=${WPS_SRC_PATH} -# WRF installation directory -WRFDir=${WRF_SRC_PATH} - -# --- Input data settings --- - -# path to geogrid input data -geogDir=${WRF_GEOG_PATH} - -# meteo input -# Vtable for the chosen meteo input -metVtableFile=${WPS_SRC_PATH}/ungrib/Variable_Tables/Vtable.GFS -# time increment in hours -metInc=1 -# full path to a met input file - you can use any "%<>" abbreviations known -# to the "date" command -metFilePattern="${WRF_GFS_METEO_PATH}/GF%Y%m%d%H" -# example: -# "/glade/p/rda/data/ds083.2/grib2/%Y/%Y.%m/fnl_%Y%m%d_%H_00.grib2" - -# --- Pre/Postprocessing settings --- - -# prepararation script -preScriptPath=NONEXISTENT.bash - -# postprocessing scripts (arbitrary) -postScriptPath=NONEXISTENT.bash -# postprocessing scripts (actions for each wrfout file) -postPerFileScriptPath=NONEXISTENT.bash - -# --- Working directories --- - -# where the WRF will be run (some fast, large disk like "scratch" or similar) -workDir=${SCRATCH}/WRF/work/ - -# where the unprocessed WRF output will be stored -stagingRootDir=${SCRATCH}/WRF/staging/ -# where the WRF output will be stored -archiveRootDir=${SCRATCH}/archive/WRF/ -# where the WRF restart files will be stored -restartRootDir=${SCRATCH}/WRF/restart/ - -# remove run directory after run is finished? -removeRunDir=false - -# --- MPI settings --- - -# mental note for GNU: -#mpirun -mca pml ucx -x UCX_TLS=rc,shm -N $SLURM_NTASKS_PER_NODE -hostfile $MACHINEFILE -mpiCommandPre="mpirun -mca pml ucx -x UCX_TLS=rc,shm -hostfile $MACHINEFILE /usr/bin/time -v" -mpiCommandMain="mpirun -mca pml ucx -x UCX_TLS=rc,shm -N $SLURM_NTASKS_PER_NODE -hostfile $MACHINEFILE /usr/bin/time -v" -mpiCommandReal=${mpiCommandPre} - -# mental note for INTEL: -#mpiCommandPre="srun /usr/bin/time -v" -#mpiCommandMain="srun /usr/bin/time -v" -#mpiCommandReal="srun /usr/bin/time -v" - -# --- Batch system --- - -# argument to submit a job in a held state -batchHoldArgument="--hold" -# command to release a held job -batchReleaseCommand="scontrol release" -# command to submit jobs to the queueing system -batchSubmitCommand=sbatch -# dependency argument for chaining runs upon submission -batchDepArgument="--dependency=afterany:__id__" -# sed command ("used as s/__command/\1/") to retrieve job run PID upon -# submission with $batchSubmitCommand -batchPidSedCommand="Submitted batch job \(.*\)" - -# --- Chemistry --- - -withChemistry=true - -# WRF-Chem installation directory -WRFChemDir=${WRF_CHEM_SRC_PATH} - -# megan_bio_emiss installation directory -WRFMEGANdir=${WRF_CHEM_MEGAN_BIO_EMISS_PATH} -# mozbc installation directory -WRFMOZARTdir=${WRF_CHEM_MOZBC_PATH} -# wesley/exocoldens installation directory -WRFmztoolsdir=${WRF_CHEM_WES_COLDENS_PATH} -# anthro_emiss installation directory -WRFanthrodir=${WRF_CHEM_ANTHRO_EMIS_PATH} -# fire_emis installation directory -WRFfiredir=${WRF_CHEM_FIRE_EMIS_PATH} - -# path to MEGAN input data -MEGANdir=/alcc/gpfs2/home/mbees/data/emissions/biogenic/MEGAN - -# use anthro_emiss or predefined files? -emissUseAnthroEmiss=true -# raw emission input - the files you read in with anthro_emiss -emissDir=/alcc/gpfs2/home/mbees/data/emissions/anthropogenic/EDGAR-HTAP/MOZART_MOSAIC -# emission conversion script for anthro_emis - must match emissions in emissDir -emissInpFile=emis_edgarhtap_mozmos.inp -# year the emissions are valid for (for offset calculation) -emissYear=2010 - -# FINN fires -fireFilePattern="/alcc/gpfs2/home/mbees/data/emissions/fires/FINN/GLOB_MOZ4_%Y%j.txt" -fireInpFile=finn_fires.inp - -# boundary condition input -chembcFilePattern="/alcc/gpfs2/home/mbees/data/chembc/WACCM/WACCM%Y%m%d" -chembcInpFile=waccm.inp - -# TUV photolysis option 4 data file -TUVDataPath="/alcc/gpfs2/home/mbees/data/tuv/TUV.phot.bz2" diff --git a/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/main b/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/main deleted file mode 100644 index cea7698faec6a40d87620c9e268710ec16c95b28..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/main +++ /dev/null @@ -1,11 +0,0 @@ -#!/bin/bash -l -#SBATCH --partition=alcc1,epyc -#SBATCH -o __runDir__/__mainJobName__.%j.%N.out -#SBATCH -D __runDir__ -#SBATCH -J __mainJobName__ -#SBATCH --nodes=2 -#SBATCH --ntasks-per-node=28 -#SBATCH --mem-per-cpu=2000 -#SBATCH --mail-type=FAIL -#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=12:00:00 diff --git a/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/post b/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/post deleted file mode 100644 index 9cecb891d441d841608a66d6cc6bd2c874488598..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/post +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash -l -#SBATCH --partition=alcc1 -#SBATCH -o __runDir__/__postJobName__.%j.%N.out -#SBATCH -D __runDir__ -#SBATCH -J __postJobName__ -#SBATCH --ntasks=1 -#SBATCH --mem=5G -#SBATCH --mail-type=FAIL -#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=02:00:00 diff --git a/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/pre b/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/pre deleted file mode 100644 index 2bb27db5fd1bfd8c733b1f17bd49b1e301a0c4b2..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/pre +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash -l -#SBATCH --partition=alcc1 -#SBATCH -o __runDir__/__preJobName__.%j.%N.out -#SBATCH -D __runDir__ -#SBATCH -J __preJobName__ -#SBATCH --ntasks=1 -#SBATCH --mem=10G -#SBATCH --mail-type=FAIL -#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=03:00:00 diff --git a/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/staging b/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/staging deleted file mode 100644 index 66b09ee6b48048058b626d7a6642151323cdc194..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/machine_specific/alcc_intel/batch_preambles/staging +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash -l -#SBATCH --partition=alcc1 -#SBATCH -o __runDir__/__stagingJobName__.%j.%N.out -#SBATCH -D __runDir__ -#SBATCH -J __stagingJobName__ -#SBATCH --ntasks=1 -#SBATCH --mem=5G -#SBATCH --mail-type=FAIL -#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=01:00:00 diff --git a/blueprints/test_chemistry/namelist.chem b/blueprints/test_chemistry/namelist.chem index 45b27140b3f4189db45688959dd8475b117046fd..402c2f41c76a1846646f59871f7de0f3bc4519b0 100644 --- a/blueprints/test_chemistry/namelist.chem +++ b/blueprints/test_chemistry/namelist.chem @@ -22,8 +22,8 @@ cldchem_onoff = 0, 0, vertmix_onoff = 1, 1, chem_conv_tr = 1, 1, - conv_tr_wetscav = 0, 0, - conv_tr_aqchem = 0, 0, + conv_tr_wetscav = 1, 1, + conv_tr_aqchem = 1, 1, seas_opt = 2, dust_opt = 3, dmsemis_opt = 1, @@ -36,7 +36,7 @@ biomass_burn_opt = 5, 5, plumerisefire_frq = 30, 30, scale_fire_emiss = .true., .true., - n2o5_hetchem = 0, + n2o5_hetchem = 1, lnox_opt = 1, 1, diagnostic_chem = 0, 0, diagnostic_dep = 0, 0, diff --git a/blueprints/test_chemistry/namelist.wps b/blueprints/test_chemistry/namelist.wps index 18e73b25d885f22fa3068a7954db15e3fd7dbdbb..bdef2dcd93c4191811119e10db5d072b53b2f724 100644 --- a/blueprints/test_chemistry/namelist.wps +++ b/blueprints/test_chemistry/namelist.wps @@ -1,6 +1,6 @@ &share wrf_core = 'ARW', - max_dom = 2, + max_dom = 1, start_date = '__startDate__','__startDate__','__startDate__' end_date = '__endDate__','__startDate__','__startDate__', interval_seconds = __metIncSec__, @@ -10,14 +10,14 @@ &geogrid parent_id = 1, 1, 2, parent_grid_ratio = 1, 9, 5, - i_parent_start = 1, 54, 40, - j_parent_start = 1, 34, 40, - e_we = 91, 82, 196, - e_sn = 93, 73, 196, + i_parent_start = 1, 96, 40, + j_parent_start = 1, 77, 40, + e_we = 46, 163, 196, + e_sn = 47, 136, 196, ! geog_data_res = 'modis_fpar+modis_lai+modis_lakes+modis_30s+modis_15s+30s', 'modis_fpar+modis_lai+modis_lakes+modis_30s+modis_15s+30s', ! 'modis_fpar+modis_lai+modis_lakes+modis_30s+modis_15s+30s', - dx = 20000, - dy = 20000, + dx = 80000, + dy = 80000, map_proj = 'lambert', ref_lat = 50.0, ref_lon = 7.5, diff --git a/blueprints/test_chemistry/namelist.wrf b/blueprints/test_chemistry/namelist.wrf index d2074fe352f3c19697dc8d3b77fb081fcf32df08..49bb51700be10922d2037a551e20771fb289d342 100644 --- a/blueprints/test_chemistry/namelist.wrf +++ b/blueprints/test_chemistry/namelist.wrf @@ -40,18 +40,18 @@ time_step = 120, time_step_fract_num = 0, time_step_fract_den = 1, - max_dom = 2, - e_we = 91, 82, 196, - e_sn = 93, 73, 196, + max_dom = 1, + e_we = 46, 163, 196, + e_sn = 47, 136, 196, e_vert = 33, 33, 33, num_metgrid_levels = 34, num_metgrid_soil_levels = 4, - dx = 20000, 2222, 400, - dy = 20000, 2222, 400, + dx = 80000, 2222, 400, + dy = 80000, 2222, 400, grid_id = 1, 2, 3, parent_id = 0, 1, 2, - i_parent_start = 1, 54, 40, - j_parent_start = 1, 34, 40, + i_parent_start = 1, 96, 40, + j_parent_start = 1, 77, 40, parent_grid_ratio = 1, 9, 5, parent_time_step_ratio = 1, 9, 5, feedback = 0, diff --git a/blueprints/test_chemistry/pp/plotting/geotiffs/run.sh b/blueprints/test_chemistry/pp/plotting/geotiffs/run.sh index 0d8604a358d8a78a5c4239d055bb270ce174a547..dbffc99908251525e4069bc210799e6757712032 100755 --- a/blueprints/test_chemistry/pp/plotting/geotiffs/run.sh +++ b/blueprints/test_chemistry/pp/plotting/geotiffs/run.sh @@ -8,7 +8,7 @@ #SBATCH --mail-user=christoph.knote@med.uni-augsburg.de #SBATCH --time=00:30:00 -module load gnu geos proj +module load gnu gdal geos proj scriptPath=/alcc/gpfs2/home/u/knotechr/wrfotron/tools/nc_2_geotiff.py wrfDataPathPattern="/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout___domain___%Y-%m-%d_%H:%M:%S" @@ -23,10 +23,10 @@ do for anHour in $(seq 0 ${intervalHours} ${rangeHours}) do - aWRFDataPath=$(date -u --date="${dateToProcess} 0 UTC +${anHour} hours" "+${wrfDataPathPattern/__domain__/${domain}}") + aWRFDataPath=$(date -u --date="${firstDateToProcess} 0 UTC +${anHour} hours" "+${wrfDataPathPattern/__domain__/${domain}}") - aPlotDir=$(date -u --date="${dateToProcess} 0 UTC +${anHour} hours" "+${plotDirPattern/__domain__/${domain}}") - aPlotPath=$(date -u --date="${dateToProcess} 0 UTC +${anHour} hours" "+${plotPathPattern}") + aPlotDir=$(date -u --date="${firstDateToProcess} 0 UTC +${anHour} hours" "+${plotDirPattern/__domain__/${domain}}") + aPlotPath=$(date -u --date="${firstDateToProcess} 0 UTC +${anHour} hours" "+${plotPathPattern}") if [ -f ${aWRFDataPath} ] then diff --git a/blueprints/test_chemistry/pp/plotting/maps/plot.py b/blueprints/test_chemistry/pp/plotting/maps/plot.py deleted file mode 100644 index e26b3cde17aa2cc343d4ab057fcbc79fb00498dc..0000000000000000000000000000000000000000 --- a/blueprints/test_chemistry/pp/plotting/maps/plot.py +++ /dev/null @@ -1,354 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -from wrf import to_np, getvar, smooth2d, get_basemap, ll_to_xy, latlon_coords, interplevel, get_cartopy, cartopy_xlim, cartopy_ylim -import numpy as np -import netCDF4 -import datetime -import matplotlib -from matplotlib import pyplot as plt -import cartopy.feature as cfeature -import cartopy.crs as ccrs -from matplotlib.offsetbox import AnchoredText -import matplotlib.gridspec as gridspec -from matplotlib.colors import LinearSegmentedColormap -import sys - -MAP_BUFFER=0.1 # 10 % buffer at edge of maps - -t2datetime = lambda x: datetime.datetime.strptime(x, '%Y-%m-%d_%H:%M:%S') - -def cartopy_lim_buffered(field, direction, buffer=MAP_BUFFER): - if direction == "x": - lims = cartopy_xlim(field) - elif direction == "y": - lims = cartopy_ylim(field) - else: - raise Exception("Moek!!") - dlim = lims[1] - lims[0] - lims[0] += 0.5 * buffer * dlim - lims[1] -= 0.5 * buffer * dlim - return lims - -def wrf_time(nc): - first = ''.join([x.decode('utf-8') for x in nc['Times'][0]]) - return t2datetime(first) - -def total_prec(nc, nc_prev): - prec_vars = ['RAINC', 'RAINSH', 'RAINNC', 'SNOWNC', 'GRAUPELNC', 'HAILNC'] - start = t2datetime(nc.START_DATE) - prec = np.sum([nc[var][0, :, :] for var in prec_vars], axis=0) -# - if start == wrf_time(nc): - return prec - else: - prec_p = np.sum([nc_prev[var][0, :, :] for var in prec_vars], axis=0) - return prec-prec_p - -def add_title_and_description(nc, title, ax): - ''' - title: including {init:s} for timestamp of run initialization, and {fcsthour:d} for forecast hour - ''' - fcstHour = int((wrf_time(nc) - t2datetime(nc.START_DATE)).total_seconds()/3600.) -# - ax.set_title(title.format(init=t2datetime(nc.START_DATE).strftime("%d.%m.%Y %H:%M UTC"), - fcsthour=fcstHour), fontdict={ 'fontsize':8, 'horizontalalignment': 'left'}, loc='left') -# - texttime = AnchoredText(wrf_time(nc).strftime("%H UTC %a"), - loc=2, pad=0.2, prop={'size': 10}, frameon=True) - ax.add_artist(texttime) -# - textcopyright = AnchoredText("$\copyright$ MBEES, Faculty of Medicine, University of Augsburg", - loc=4, pad=0.1, prop={'size': 6, 'color':'white'}, frameon=False) - ax.add_artist(textcopyright) - -def make_fig_and_ax(plt, crs, add_horiz_cbar=False): - fig = plt.figure(figsize=(7,7)) -# - gs = gridspec.GridSpec(2, 2, width_ratios=[20,1], height_ratios=[20,1], wspace=0.05, hspace=0.05) -# - ax = fig.add_subplot(gs[0,0], projection=crs) - ax_cbarv = fig.add_subplot(gs[0,1]) - ax_cbarh = None - if add_horiz_cbar: - ax_cbarh = fig.add_subplot(gs[1,0]) - return fig, ax, ax_cbarv, ax_cbarh - -def slp_map(t, nc, nc_prev, output_path, plevels=np.arange(950, 1050, 5)): -# - slp = getvar(nc, "slp") - mid_clouds = getvar(nc, "mid_cloudfrac") - prec = total_prec(nc, nc_prev) -# - # Smooth the sea level pressure since it tends to be noisy near the - # mountains - smooth_slp = smooth2d(slp, 5, cenweight=4) - # Get the latitude and longitude points - lats, lons = latlon_coords(slp) -# - # Get the WRF projection (CRS) object - crs = get_cartopy(slp) -# - # Create a figure - fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs, add_horiz_cbar=True) -# - # Draw the contours and filled contours - cmap = LinearSegmentedColormap.from_list("", ["white", "silver"]) - cmap2 = matplotlib.cm.viridis_r -# - con = ax.contourf(to_np(lons), to_np(lats), to_np(mid_clouds), - levels=[1/10, 1/4, 1/2, 3/4, 1], cmap=cmap, zorder=2, - transform=ccrs.PlateCarree()) - con2 = ax.contourf(to_np(lons), to_np(lats), to_np(prec), - levels=[0.5, 1, 2, 5, 10, 25, 50], cmap=cmap2, zorder=2, extend='max', - transform=ccrs.PlateCarree()) -# - linewidths = np.ones(plevels.shape) - linewidths[np.where(plevels == 1015)[0]] = 3 - cs = ax.contour(to_np(lons), to_np(lats), to_np(smooth_slp), levels=plevels, - linewidths=linewidths, colors="black", zorder=4, - transform=ccrs.PlateCarree()) - plt.clabel(cs, fmt='%1.0f', fontsize=8) -# - ax.coastlines(linewidth=0.3, zorder=2) - countries = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', scale='50m') - ax.add_feature(countries, facecolor='none', edgecolor='black', linewidth=0.3, zorder=3) -# - cbar = plt.colorbar(con, cax=ax_cbarv, ticks=[1/10, 1/4, 1/2, 3/4, 1]) - cbar2 = plt.colorbar(con2, cax=ax_cbarh, orientation="horizontal") -# - try: - cbar.ax.set_yticklabels(["1/10", "1/4", "1/2", "3/4", "1"], fontsize=10) - except: - print("Plotting tick labels failed") -# - ax.set_xlim(cartopy_lim_buffered(slp, "x")) - ax.set_ylim(cartopy_lim_buffered(slp, "y")) - - add_title_and_description(nc, "Sea Level Pressure [hPa], Mid-Level Cloud Cover, Total prec. in the last hour [mm] \n" - "WRF {init:s} +{fcsthour:d}", ax) -# - plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) - plt.close() - - -def plevel_map(t, nc, output_path, plevel, phi_contour_intervals=np.arange(400, 800, 8)): -# - def interp_var(nc, variablename, p, plevel): - var = getvar(nc, variablename) - return interplevel(var, p, plevel) -# - p = getvar(nc, "p", units='hPa') -# - geop = interp_var(nc, "geopt", p, plevel) - geop = smooth2d(geop, 5, cenweight=4) -# - temp = interp_var(nc, "tc", p, plevel) -# - u, v = interp_var(nc, "uvmet", p, plevel) -# - # Get the latitude and longitude points - lats, lons = latlon_coords(temp) -# - # Get the WRF projection (CRS) object - crs = get_cartopy(temp) -# - # Create a figure - fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs) -# - # Temperature - lab = np.array([-46, -41, -36, -31, -26, -21, -16, -11, -6, -1, 4]) - con = ax.contourf(to_np(lons), to_np(lats), temp, levels=lab, - transform=ccrs.PlateCarree(), - cmap='viridis', zorder=0, extend='both') -# - # Geopotential height - linewidths = np.ones(phi_contour_intervals.shape) - linewidths[np.where(phi_contour_intervals == 552)[0]] = 3 - cs = ax.contour(to_np(lons), to_np(lats), geop/98.0, levels=phi_contour_intervals, - linewidths=linewidths, colors="black", - zorder=1, - transform=ccrs.PlateCarree()) - plt.clabel(cs, fmt='%1.0f', fontsize=10) -# -# - ax.coastlines(linewidth=0.3, zorder=2) - countries = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', scale='50m') - ax.add_feature(countries, facecolor='none', edgecolor='black', linewidth=0.3, zorder=3) -# - nb = 20 - brsb = ax.barbs(to_np(lons[::nb,::nb]), to_np(lats[::nb,::nb]), - to_np(u[::nb,::nb]), to_np(v[::nb,::nb]), - transform=ccrs.PlateCarree(), - length=4.5, zorder=4, linewidth=0.3) -# - ax.set_xlim(cartopy_lim_buffered(temp, "x")) - ax.set_ylim(cartopy_lim_buffered(temp, "y")) -# - cbar = plt.colorbar(con, cax=ax_cbarv) - add_title_and_description(nc, "{:d}".format(plevel) + " hPa level Geopotential [gpdm], Temperature [°C] and Wind \n" - "WRF {init:s} +{fcsthour:d}", ax) -# - plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) - plt.close() - -def sfc_map(t, nc, get_var_fun, levels, ticks, label, pois, output_path, cmap="viridis", extend='both'): -# - var = get_var_fun(nc) - u, v = getvar(nc, 'uvmet10') -# - # Get the latitude and longitude points - lats, lons = latlon_coords(u) -# - # Get the WRF projection (CRS) object - crs = get_cartopy(u) -# - # Create a figure - fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs) -# - # Draw the contours and filled contours - lab = np.array(levels) - con = ax.contourf(to_np(lons), to_np(lats), var, levels=lab, - transform=ccrs.PlateCarree(), - cmap=cmap, zorder=0, extend=extend) -# - linewidths = np.ones(lab.shape)/2 - linewidths[np.where(lab == 0)[0]] = 1 - cs = ax.contour(to_np(lons), to_np(lats), var, levels=lab[::2], - linewidths=linewidths, colors="black", alpha=0.4, - zorder=1, - transform=ccrs.PlateCarree()) - plt.clabel(cs, fmt='%1.0f', fontsize=8) -# - ax.coastlines(linewidth=0.3, zorder=2) - countries = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', scale='50m') - ax.add_feature(countries, facecolor='none', edgecolor='black', linewidth=0.3, zorder=3) -# - nb = 10 - brsb = ax.barbs(to_np(lons[::nb,::nb]), to_np(lats[::nb,::nb]), - to_np(u[::nb,::nb]), to_np(v[::nb,::nb]), - transform=ccrs.PlateCarree(), - length=4.5, zorder=4, linewidth=0.3) -# - def plot_poi(name, lon, lat, ms, fontsize): - ax.plot(lon, lat, 'o', zorder=5, color='darkred', ms=ms, transform=ccrs.PlateCarree()) - ax.annotate(name, xy=(lon, lat), xytext=(5,5), textcoords="offset points", - fontsize=fontsize, zorder=6, color='darkred', - xycoords=ccrs.PlateCarree()._as_mpl_transform(ax)) -# - nul = [ plot_poi(**item) for item in pois ] -# - ax.set_xlim(cartopy_lim_buffered(u, "x")) - ax.set_ylim(cartopy_lim_buffered(u, "y")) -# - cbar = plt.colorbar(con, cax=ax_cbarv, extend=extend) - cbar.set_ticks(ticks) -# - add_title_and_description(nc, label + " and Wind (10m) \n" + "WRF {init:s} +{fcsthour:d}", ax) -# - plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) - plt.close() - -# -------------------------------------------------------------------------------------------------------- -if __name__ == '__main__': - - import datetime - - from argparse import ArgumentParser, ArgumentTypeError - - def valid_date(s): - try: - return datetime.datetime.strptime(s, "%Y-%m-%d") - except ValueError: - raise ArgumentTypeError("Not a valid date: '{0}'.".format(s)) - - parser = ArgumentParser() - parser.add_argument('wrf_data_fpath_pattern', type=str, help='WRF input data path pattern (including strptime pattern for date and {domain:s} for domain (01, 02, ...))') - parser.add_argument('date', type=valid_date, help='Date to process (YYYY-mm-dd)') - parser.add_argument('interval_hours', type=int, help='Time interval in hours for plotting') - parser.add_argument('fcst_time_hours', type=int, help='Forecast time in hours') - parser.add_argument('plot_fpath_pattern', type=str, help='Plot output path pattern (including strptime pattern for date and {plot:s} for plot identifier (e.g., 500hPa_europe...)') - - args = parser.parse_args() -# args = parser.parse_args( [ '/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout_d{domain:02d}_%Y-%m-%d_%H:%M:%S', '2021-02-17', '1', '72', '/alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/fcst_movie/{plot:s}_%Y%m%d%H.png' ]) - times = [ args.date + datetime.timedelta(hours=t) for t in range(0, args.fcst_time_hours, args.interval_hours) ] - - for t in times: - tprev = t + datetime.timedelta(hours=-1) - - input_path_d01 = t.strftime(args.wrf_data_fpath_pattern).format(domain=1) - previnp_path_d01 = tprev.strftime(args.wrf_data_fpath_pattern).format(domain=1) - - try: - nc_d01 = netCDF4.Dataset(input_path_d01) - nc_d01_prev = netCDF4.Dataset(previnp_path_d01) - except Exception as e: - print("Not all data available for step {:s}: {:s}".format(t.strftime("%Y-%m-%d %H"), str(e))) - continue - - print(t) - - pois_southern_germany = [ { 'name':'Munich', 'lon':11.582, 'lat':48.135, 'ms':5, 'fontsize':9 }, - { 'name':'Augsburg', 'lon':10.894446, 'lat':48.366512, 'ms':2.5, 'fontsize':8 }, - { 'name':'Regensburg', 'lon':12.101624, 'lat':49.013432, 'ms':2.5, 'fontsize':8 }, - { 'name':'Ingolstadt', 'lon':11.426311, 'lat':48.761423, 'ms':2.5, 'fontsize':8 }, - { 'name':'Innsbruck', 'lon':11.388397, 'lat':47.258320, 'ms':2.5, 'fontsize':8 } - ] - - # 500 hPa Europe domain - output_path = t.strftime(args.plot_fpath_pattern).format(plot='500hPa_Europe') - plevel_map(t, nc_d01, output_path, 500) - - # Surface Europe domain - # temperature - sfc_map(t, nc_d01, lambda nc: getvar(nc, 'T2') - 273.15, - levels=[-25, -20, -15, -12.5, -10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30, 32.5, 35], - ticks=[-25, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35], - label="Temperature (2m) [°C]", pois={}, - output_path=t.strftime(args.plot_fpath_pattern).format(plot='T_Europe')) - # Ozone - sfc_map(t, nc_d01, lambda nc: getvar(nc, 'o3')[0,:,:] * 1e3, - levels=range(0, 160, 10), - ticks=range(0, 150, 20), - label="Ozone (ground level) [ppbv]", pois={}, - output_path=t.strftime(args.plot_fpath_pattern).format(plot='O3_Europe'), - extend='max') - # NO2 - sfc_map(t, nc_d01, lambda nc: getvar(nc, 'no2')[0,:,:] * 1e3, - levels=range(0, 110, 10), - ticks=range(0, 110, 20), - label="NO$_2$ (ground level) [ppbv]", pois={}, - output_path=t.strftime(args.plot_fpath_pattern).format(plot='NO2_Europe'), - extend='max') - # PM - def get_pm(nc, bins=[ 1, 2, 3 ]) : - pmspecs = [ "so4", "no3", "smpa", "smpbb", "glysoa_sfc", "biog1_c", "biog1_o", "cl", "co3", "nh4", "na", "cl", "oin", "oc", "bc", "water" ] # ug kg-1 - alt = getvar(nc, "ALT")[0,:,:] # inverse density m3 kg-1 - all = [ getvar(nc, spec + "_a{:02d}".format(bin))[0,:,:] / alt for spec in pmspecs for bin in bins ] - return np.sum(all, axis=0) - - sfc_map(t, nc_d01, lambda nc: get_pm(nc, bins=[1]), - levels=range(0, 130, 10), - ticks=range(0, 130, 20), - label="UFP (ground level) [$\mu$g m$^{{-3}}$]", pois={}, - output_path=t.strftime(args.plot_fpath_pattern).format(plot='UFP_Europe'), - extend='max') - sfc_map(t, nc_d01, lambda nc: get_pm(nc, bins=[1,2,3]), - levels=range(0, 130, 10), - ticks=range(0, 130, 20), - label="PM$_{{2.5}}$ (ground level) [$\mu$g m$^{{-3}}$]", pois={}, - output_path=t.strftime(args.plot_fpath_pattern).format(plot='PM2_5_Europe'), - extend='max') - sfc_map(t, nc_d01, lambda nc: get_pm(nc, bins=[1,2,3,4]), - levels=range(0, 130, 10), - ticks=range(0, 130, 20), - label="PM$_{{10}}$ (ground level) [$\mu$g m$^{{-3}}$]", pois={}, - output_path=t.strftime(args.plot_fpath_pattern).format(plot='PM10_Europe'), - extend='max') - # SLP Europe - output_path = t.strftime(args.plot_fpath_pattern).format(plot='SLP_Europe') - slp_map(t, nc_d01, nc_d01_prev, output_path) - - nc_d01.close() - nc_d01_prev.close() - diff --git a/blueprints/test_chemistry/pp/plotting/maps/plots.py b/blueprints/test_chemistry/pp/plotting/maps/plots.py new file mode 100644 index 0000000000000000000000000000000000000000..b44e3a1b5685f90abbf1644ac228271eeb429cea --- /dev/null +++ b/blueprints/test_chemistry/pp/plotting/maps/plots.py @@ -0,0 +1,197 @@ + +from wrf import to_np, getvar, smooth2d, latlon_coords, interplevel, get_cartopy +import numpy as np +import matplotlib +from matplotlib import pyplot as plt +import cartopy.feature as cfeature +import cartopy.crs as ccrs +from matplotlib.colors import LinearSegmentedColormap, ListedColormap + + +from tools import total_prec, cartopy_lim_buffered, make_fig_and_ax, add_title_and_description + +def slp_map(t, nc, prec_prev, output_path, plevels=np.arange(950, 1050, 5)): +# + slp = getvar(nc, "slp") + mid_clouds = getvar(nc, "mid_cloudfrac") + prec = total_prec(nc, prec_prev) +# + # Smooth the sea level pressure since it tends to be noisy near the + # mountains + smooth_slp = smooth2d(slp, 5, cenweight=4) + # Get the latitude and longitude points + lats, lons = latlon_coords(slp) +# + # Get the WRF projection (CRS) object + crs = get_cartopy(slp) +# + # Create a figure + fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs, add_horiz_cbar=True) +# + # Draw the contours and filled contours + cmap = LinearSegmentedColormap.from_list("", ["white", "silver"]) + cmap2 = matplotlib.cm.viridis_r +# + con = ax.contourf(to_np(lons), to_np(lats), to_np(mid_clouds), + levels=[1/10, 1/4, 1/2, 3/4, 1], cmap=cmap, zorder=2, + transform=ccrs.PlateCarree()) + con2 = ax.contourf(to_np(lons), to_np(lats), to_np(prec), + levels=[0.5, 1, 2, 5, 10, 25, 50], cmap=cmap2, zorder=2, extend='max', + transform=ccrs.PlateCarree()) +# + linewidths = np.ones(plevels.shape) + linewidths[np.where(plevels == 1015)[0]] = 3 + cs = ax.contour(to_np(lons), to_np(lats), to_np(smooth_slp), levels=plevels, + linewidths=linewidths, colors="black", zorder=4, + transform=ccrs.PlateCarree()) + plt.clabel(cs, fmt='%1.0f', fontsize=8) +# + ax.coastlines(linewidth=0.3, zorder=2) + countries = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', scale='50m') + ax.add_feature(countries, facecolor='none', edgecolor='black', linewidth=0.3, zorder=3) +# + cbar = plt.colorbar(con, cax=ax_cbarv, ticks=[1/10, 1/4, 1/2, 3/4, 1]) + cbar2 = plt.colorbar(con2, cax=ax_cbarh, orientation="horizontal") +# + try: + cbar.ax.set_yticklabels(["1/10", "1/4", "1/2", "3/4", "1"], fontsize=10) + except: + print("Plotting tick labels failed") +# + ax.set_xlim(cartopy_lim_buffered(slp, "x")) + ax.set_ylim(cartopy_lim_buffered(slp, "y")) + + add_title_and_description(nc, "Sea Level Pressure [hPa], Mid-Level Cloud Cover, Total prec. in the last hour [mm] \n" + "WRF {init:s} +{fcsthour:d}", ax) +# + plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) + plt.close() + + return prec + + +def plevel_map(t, nc, output_path, plevel, phi_contour_intervals=np.arange(400, 800, 8)): +# + def interp_var(nc, variablename, p, plevel): + var = getvar(nc, variablename) + return interplevel(var, p, plevel) +# + p = getvar(nc, "p", units='hPa') +# + geop = interp_var(nc, "geopt", p, plevel) + geop = smooth2d(geop, 5, cenweight=4) +# + temp = interp_var(nc, "tc", p, plevel) +# + u, v = interp_var(nc, "uvmet", p, plevel) +# + # Get the latitude and longitude points + lats, lons = latlon_coords(temp) +# + # Get the WRF projection (CRS) object + crs = get_cartopy(temp) +# + # Create a figure + fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs) +# + # Temperature + lab = np.array([-46, -41, -36, -31, -26, -21, -16, -11, -6, -1, 4]) + con = ax.contourf(to_np(lons), to_np(lats), temp, levels=lab, + transform=ccrs.PlateCarree(), + cmap='viridis', zorder=0, extend='both') +# + # Geopotential height + linewidths = np.ones(phi_contour_intervals.shape) + linewidths[np.where(phi_contour_intervals == 552)[0]] = 3 + cs = ax.contour(to_np(lons), to_np(lats), geop/98.0, levels=phi_contour_intervals, + linewidths=linewidths, colors="black", + zorder=1, + transform=ccrs.PlateCarree()) + plt.clabel(cs, fmt='%1.0f', fontsize=10) +# +# + ax.coastlines(linewidth=0.3, zorder=2) + countries = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', scale='50m') + ax.add_feature(countries, facecolor='none', edgecolor='black', linewidth=0.3, zorder=3) +# + nb = 20 + brsb = ax.barbs(to_np(lons[::nb,::nb]), to_np(lats[::nb,::nb]), + to_np(u[::nb,::nb]), to_np(v[::nb,::nb]), + transform=ccrs.PlateCarree(), + length=4.5, zorder=4, linewidth=0.3) +# + ax.set_xlim(cartopy_lim_buffered(temp, "x")) + ax.set_ylim(cartopy_lim_buffered(temp, "y")) +# + cbar = plt.colorbar(con, cax=ax_cbarv) + add_title_and_description(nc, "{:d}".format(plevel) + " hPa level Geopotential [gpdm], Temperature [°C] and Wind \n" + "WRF {init:s} +{fcsthour:d}", ax) +# + plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) + plt.close() + +# https://stackoverflow.com/questions/37327308/add-alpha-to-an-existing-matplotlib-colormap +def alphafy_colormap(name, mina=0.0, maxa=0.8): + # Choose colormap + cmap = matplotlib.cm.get_cmap(name) + + # Get the colormap colors + my_cmap = cmap(np.arange(cmap.N)) + + # Set alpha + my_cmap[:,-1] = np.linspace(mina, maxa, cmap.N) + + # Create new colormap + return ListedColormap(my_cmap) + + +def sfc_map(t, nc, get_var_fun, levels, ticks, label, bg_function, output_path, cmap="YlOrRd", extend='both'): +# + var = get_var_fun(nc) + u, v = getvar(nc, 'uvmet10') +# + # Get the latitude and longitude points + lats, lons = latlon_coords(u) +# + # Get the WRF projection (CRS) object + crs = get_cartopy(u) +# + # Create a figure + fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs) +# + # Draw the contours and filled contours + lab = np.array(levels) + con = ax.contourf(to_np(lons), to_np(lats), var, levels=lab, + transform=ccrs.PlateCarree(), + cmap=alphafy_colormap(cmap), zorder=1, extend=extend) +# + linewidths = np.ones(lab.shape)/2 + linewidths[np.where(lab == 0)[0]] = 1 + cs = ax.contour(to_np(lons), to_np(lats), var, levels=lab[::2], + linewidths=linewidths, colors="black", alpha=0.5, + zorder=2, + transform=ccrs.PlateCarree()) + plt.clabel(cs, fmt='%1.0f', fontsize=8) +# + # Add background (needs to happen after map is already created...) + bg_function(ax, crs) +# + nb = 10 + brsb = ax.barbs(to_np(lons[::nb,::nb]), to_np(lats[::nb,::nb]), + to_np(u[::nb,::nb]), to_np(v[::nb,::nb]), + transform=ccrs.PlateCarree(), + length=4.5, zorder=4, linewidth=0.5) + + xlim = cartopy_lim_buffered(u, "x") + ylim = cartopy_lim_buffered(u, "y") +# + ax.set_xlim(xlim) + ax.set_ylim(ylim) +# + cbar = plt.colorbar(con, cax=ax_cbarv) + cbar.set_ticks(ticks) +# + add_title_and_description(nc, label + " and Wind (10m) \n" + "WRF {init:s} +{fcsthour:d}", ax) +# + plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) + plt.close() diff --git a/blueprints/test_chemistry/pp/plotting/maps/pois.py b/blueprints/test_chemistry/pp/plotting/maps/pois.py new file mode 100644 index 0000000000000000000000000000000000000000..68244521e48e31c252db0ee3bdb74106d6c2ff1b --- /dev/null +++ b/blueprints/test_chemistry/pp/plotting/maps/pois.py @@ -0,0 +1,268 @@ +southern_germany = [ { 'name':'Munich', 'lon':11.582, 'lat':48.135}, + { 'name':'Salzburg', 'lon':13.0550, 'lat':47.8095}, + { 'name':'Augsburg', 'lon':10.8978, 'lat':48.3705}, + { 'name':'Passau', 'lon':13.4319, 'lat':48.5667}, + { 'name':'Rosenheim', 'lon':12.1181, 'lat':47.8571}, + { 'name':'Ingolstadt', 'lon':11.4258, 'lat':48.7665}, + { 'name':'Garmisch-Partenkirchen', 'lon':11.0948, 'lat':47.4919}, + { 'name':'Ulm', 'lon':9.9876, 'lat':48.4011}, + { 'name':'Stuttgart', 'lon':9.1829, 'lat':48.7758}, + { 'name':'Konstanz', 'lon':9.1737, 'lat':47.6780}, + { 'name':'Zürich', 'lon':8.5417, 'lat':47.3769}, + { 'name':'Nürnberg', 'lon':11.0767, 'lat':49.4521}, + { 'name':'Regensburg', 'lon':12.101624, 'lat':49.013432}, + { 'name':'Innsbruck', 'lon':11.388397, 'lat':47.258320} ] + + +capitals = [ { 'name':'Sukhumi', 'lat':43.001525, 'lon':41.023415}, + { 'name':'Kabul', 'lat':34.575503, 'lon':69.240073}, + { 'name':'Mariehamn', 'lat':60.1, 'lon':19.933333}, + { 'name':'Tirana', 'lat':41.327546, 'lon':19.818698}, + { 'name':'Algiers', 'lat':36.752887, 'lon':3.042048}, + { 'name':'Pago Pago', 'lat':-14.275632, 'lon':-170.702036}, + { 'name':'Andorra la Vella', 'lat':42.506317, 'lon':1.521835}, + { 'name':'Luanda', 'lat':-8.839988, 'lon':13.289437}, + { 'name':'The Valley', 'lat':18.214813, 'lon':-63.057441}, + { 'name':'South Pole', 'lat':-90, 'lon':0}, + { 'name':'St. John\'s', 'lat':17.12741, 'lon':-61.846772}, + { 'name':'Buenos Aires', 'lat':-34.603684, 'lon':-58.381559}, + { 'name':'Yerevan', 'lat':40.179186, 'lon':44.499103}, + { 'name':'Oranjestad', 'lat':12.509204, 'lon':-70.008631}, + { 'name':'Canberra', 'lat':-35.282, 'lon':149.128684}, + { 'name':'Vienna', 'lat':48.208174, 'lon':16.373819}, + { 'name':'Baku', 'lat':40.409262, 'lon':49.867092}, + { 'name':'Nassau', 'lat':25.047984, 'lon':-77.355413}, + { 'name':'Manama', 'lat':26.228516, 'lon':50.58605}, + { 'name':'Dhaka', 'lat':23.810332, 'lon':90.412518}, + { 'name':'Bridgetown', 'lat':13.113222, 'lon':-59.598809}, + { 'name':'Minsk', 'lat':53.90454, 'lon':27.561524}, + { 'name':'Brussels', 'lat':50.85034, 'lon':4.35171}, + { 'name':'Belmopan', 'lat':17.251011, 'lon':-88.75902}, + { 'name':'Porto-Novo', 'lat':6.496857, 'lon':2.628852}, + { 'name':'Hamilton', 'lat':32.294816, 'lon':-64.781375}, + { 'name':'Thimphu', 'lat':27.472792, 'lon':89.639286}, + { 'name':'La Paz', 'lat':-16.489689, 'lon':-68.119294}, + { 'name':'Sarajevo', 'lat':43.856259, 'lon':18.413076}, + { 'name':'Gaborone', 'lat':-24.628208, 'lon':25.923147}, + { 'name':'Bouvet Island', 'lat':-54.43, 'lon':3.38}, + { 'name':'Brasília', 'lat':-15.794229, 'lon':-47.882166}, + { 'name':'Camp Justice', 'lat':21.3419, 'lon':55.4778}, + { 'name':'Road Town', 'lat':18.428612, 'lon':-64.618466}, + { 'name':'Bandar Seri Begawan', 'lat':4.903052, 'lon':114.939821}, + { 'name':'Sofia', 'lat':42.697708, 'lon':23.321868}, + { 'name':'Ouagadougou', 'lat':12.371428, 'lon':-1.51966}, + { 'name':'Bujumbura', 'lat':-3.361378, 'lon':29.359878}, + { 'name':'Phnom Penh', 'lat':11.544873, 'lon':104.892167}, + { 'name':'Yaoundé', 'lat':3.848033, 'lon':11.502075}, + { 'name':'Ottawa', 'lat':45.42153, 'lon':-75.697193}, + { 'name':'Praia', 'lat':14.93305, 'lon':-23.513327}, + { 'name':'George Town', 'lat':19.286932, 'lon':-81.367439}, + { 'name':'Bangui', 'lat':4.394674, 'lon':18.55819}, + { 'name':'N\'Djamena', 'lat':12.134846, 'lon':15.055742}, + { 'name':'Santiago', 'lat':-33.44889, 'lon':-70.669265}, + { 'name':'Beijing', 'lat':39.904211, 'lon':116.407395}, + { 'name':'Flying Fish Cove', 'lat':-10.420686, 'lon':105.679379}, + { 'name':'West Island', 'lat':-12.188834, 'lon':96.829316}, + { 'name':'Bogotá', 'lat':4.710989, 'lon':-74.072092}, + { 'name':'Moroni', 'lat':-11.717216, 'lon':43.247315}, + { 'name':'Kinshasa', 'lat':-4.441931, 'lon':15.266293}, + { 'name':'Brazzaville', 'lat':-4.26336, 'lon':15.242885}, + { 'name':'Avarua', 'lat':-21.212901, 'lon':-159.782306}, + { 'name':'San José', 'lat':9.928069, 'lon':-84.090725}, + { 'name':'Yamoussoukro', 'lat':6.827623, 'lon':-5.289343}, + { 'name':'Zagreb', 'lat':45.815011, 'lon':15.981919}, + { 'name':'Havana', 'lat':23.05407, 'lon':-82.345189}, + { 'name':'Willemstad', 'lat':12.122422, 'lon':-68.882423}, + { 'name':'Nicosia', 'lat':35.185566, 'lon':33.382276}, + { 'name':'Prague', 'lat':50.075538, 'lon':14.4378}, + { 'name':'Copenhagen', 'lat':55.676097, 'lon':12.568337}, + { 'name':'Djibouti', 'lat':11.572077, 'lon':43.145647}, + { 'name':'Roseau', 'lat':15.309168, 'lon':-61.379355}, + { 'name':'Santo Domingo', 'lat':18.486058, 'lon':-69.931212}, + { 'name':'Quito', 'lat':-0.180653, 'lon':-78.467838}, + { 'name':'Cairo', 'lat':30.04442, 'lon':31.235712}, + { 'name':'San Salvador', 'lat':13.69294, 'lon':-89.218191}, + { 'name':'Malabo', 'lat':3.750412, 'lon':8.737104}, + { 'name':'Asmara', 'lat':15.322877, 'lon':38.925052}, + { 'name':'Tallinn', 'lat':59.436961, 'lon':24.753575}, + { 'name':'Addis Ababa', 'lat':8.980603, 'lon':38.757761}, + { 'name':'Stanley', 'lat':-51.697713, 'lon':-57.851663}, + { 'name':'Tórshavn', 'lat':62.007864, 'lon':-6.790982}, + { 'name':'Suva', 'lat':-18.124809, 'lon':178.450079}, + { 'name':'Helsinki', 'lat':60.173324, 'lon':24.941025}, + { 'name':'Paris', 'lat':48.856614, 'lon':2.352222}, + { 'name':'Cayenne', 'lat':4.92242, 'lon':-52.313453}, + { 'name':'Papeete', 'lat':-17.551625, 'lon':-149.558476}, + { 'name':'Saint-Pierre ', 'lat':-21.3419, 'lon':55.4778}, + { 'name':'Libreville', 'lat':0.416198, 'lon':9.467268}, + { 'name':'Banjul', 'lat':13.454876, 'lon':-16.579032}, + { 'name':'Tbilisi', 'lat':41.715138, 'lon':44.827096}, + { 'name':'Berlin', 'lat':52.520007, 'lon':13.404954}, + { 'name':'Accra', 'lat':5.603717, 'lon':-0.186964}, + { 'name':'Gibraltar', 'lat':36.140773, 'lon':-5.353599}, + { 'name':'Athens', 'lat':37.983917, 'lon':23.72936}, + { 'name':'Nuuk', 'lat':64.18141, 'lon':-51.694138}, + { 'name':'St. George\'s', 'lat':12.056098, 'lon':-61.7488}, + { 'name':'Basse-Terre', 'lat':16.014453, 'lon':-61.706411}, + { 'name':'Hagåtña', 'lat':13.470891, 'lon':144.751278}, + { 'name':'Guatemala City', 'lat':14.634915, 'lon':-90.506882}, + { 'name':'St. Peter Port', 'lat':49.455443, 'lon':-2.536871}, + { 'name':'Conakry', 'lat':9.641185, 'lon':-13.578401}, + { 'name':'Bissau', 'lat':11.881655, 'lon':-15.617794}, + { 'name':'Georgetown', 'lat':6.801279, 'lon':-58.155125}, + { 'name':'Port-au-Prince', 'lat':18.594395, 'lon':-72.307433}, + { 'name':'Tegucigalpa', 'lat':14.072275, 'lon':-87.192136}, + { 'name':'Hong Kong', 'lat':22.396428, 'lon':114.109497}, + { 'name':'Budapest', 'lat':47.497912, 'lon':19.040235}, + { 'name':'Reykjavík', 'lat':64.126521, 'lon':-21.817439}, + { 'name':'New Delhi', 'lat':28.613939, 'lon':77.209021}, + { 'name':'Jakarta', 'lat':-6.208763, 'lon':106.845599}, + { 'name':'Tehran', 'lat':35.689198, 'lon':51.388974}, + { 'name':'Baghdad', 'lat':33.312806, 'lon':44.361488}, + { 'name':'Dublin', 'lat':53.349805, 'lon':-6.26031}, + { 'name':'Douglas', 'lat':54.152337, 'lon':-4.486123}, + { 'name':'Tel Aviv', 'lat':32.0853, 'lon':34.781768}, + { 'name':'Rome', 'lat':41.902784, 'lon':12.496366}, + { 'name':'Kingston', 'lat':18.042327, 'lon':-76.802893}, + { 'name':'Tokyo', 'lat':35.709026, 'lon':139.731992}, + { 'name':'St. Helier', 'lat':49.186823, 'lon':-2.106568}, + { 'name':'Amman', 'lat':31.956578, 'lon':35.945695}, + { 'name':'Astana', 'lat':51.160523, 'lon':71.470356}, + { 'name':'Nairobi', 'lat':-1.292066, 'lon':36.821946}, + { 'name':'Tarawa Atoll', 'lat':1.451817, 'lon':172.971662}, + { 'name':'Pristina', 'lat':42.662914, 'lon':21.165503}, + { 'name':'Kuwait City', 'lat':29.375859, 'lon':47.977405}, + { 'name':'Bishkek', 'lat':42.874621, 'lon':74.569762}, + { 'name':'Vientiane', 'lat':17.975706, 'lon':102.633104}, + { 'name':'Riga', 'lat':56.949649, 'lon':24.105186}, + { 'name':'Beirut', 'lat':33.888629, 'lon':35.495479}, + { 'name':'Maseru', 'lat':-29.363219, 'lon':27.51436}, + { 'name':'Monrovia', 'lat':6.290743, 'lon':-10.760524}, + { 'name':'Tripoli', 'lat':32.887209, 'lon':13.191338}, + { 'name':'Vaduz', 'lat':47.14103, 'lon':9.520928}, + { 'name':'Vilnius', 'lat':54.687156, 'lon':25.279651}, + { 'name':'Luxembourg', 'lat':49.611621, 'lon':6.131935}, + { 'name':'Macau', 'lat':22.166667, 'lon':113.55}, + { 'name':'Skopje', 'lat':41.997346, 'lon':21.427996}, + { 'name':'Antananarivo', 'lat':-18.87919, 'lon':47.507905}, + { 'name':'Lilongwe', 'lat':-13.962612, 'lon':33.774119}, + { 'name':'Kuala Lumpur', 'lat':3.139003, 'lon':101.686855}, + { 'name':'Malé', 'lat':4.175496, 'lon':73.509347}, + { 'name':'Bamako', 'lat':12.639232, 'lon':-8.002889}, + { 'name':'Valletta', 'lat':35.898909, 'lon':14.514553}, + { 'name':'Majuro', 'lat':7.116421, 'lon':171.185774}, + { 'name':'Fort-de-France', 'lat':14.616065, 'lon':-61.05878}, + { 'name':'Nouakchott', 'lat':18.07353, 'lon':-15.958237}, + { 'name':'Port Louis', 'lat':-20.166896, 'lon':57.502332}, + { 'name':'Mamoudzou', 'lat':-12.780949, 'lon':45.227872}, + { 'name':'Mexico City', 'lat':19.432608, 'lon':-99.133208}, + { 'name':'Palikir', 'lat':6.914712, 'lon':158.161027}, + { 'name':'Chisinau', 'lat':47.010453, 'lon':28.86381}, + { 'name':'Monaco', 'lat':43.737411, 'lon':7.420816}, + { 'name':'Ulaanbaatar', 'lat':47.886399, 'lon':106.905744}, + { 'name':'Podgorica', 'lat':42.43042, 'lon':19.259364}, + { 'name':'Plymouth', 'lat':16.706523, 'lon':-62.215738}, + { 'name':'Rabat', 'lat':33.97159, 'lon':-6.849813}, + { 'name':'Maputo', 'lat':-25.891968, 'lon':32.605135}, + { 'name':'Naypyidaw', 'lat':19.763306, 'lon':96.07851}, + { 'name':'Stepanakert', 'lat':39.826385, 'lon':46.763595}, + { 'name':'Windhoek', 'lat':-22.560881, 'lon':17.065755}, + { 'name':'Yaren', 'lat':-0.546686, 'lon':166.921091}, + { 'name':'Kathmandu', 'lat':27.717245, 'lon':85.323961}, + { 'name':'Amsterdam', 'lat':52.370216, 'lon':4.895168}, + { 'name':'Willemstad ', 'lat':12.1091242, 'lon':-68.9316546}, + { 'name':'Nouméa', 'lat':-22.255823, 'lon':166.450524}, + { 'name':'Wellington', 'lat':-41.28646, 'lon':174.776236}, + { 'name':'Managua', 'lat':12.114993, 'lon':-86.236174}, + { 'name':'Niamey', 'lat':13.511596, 'lon':2.125385}, + { 'name':'Abuja', 'lat':9.076479, 'lon':7.398574}, + { 'name':'Alofi', 'lat':-19.055371, 'lon':-169.917871}, + { 'name':'Kingston', 'lat':-29.056394, 'lon':167.959588}, + { 'name':'Pyongyang', 'lat':39.039219, 'lon':125.762524}, + { 'name':'Nicosia', 'lat':35.185566, 'lon':33.382276}, + { 'name':'Saipan', 'lat':15.177801, 'lon':145.750967}, + { 'name':'Oslo', 'lat':59.913869, 'lon':10.752245}, + { 'name':'Muscat', 'lat':23.58589, 'lon':58.405923}, + { 'name':'Islamabad', 'lat':33.729388, 'lon':73.093146}, + { 'name':'Ngerulmud', 'lat':7.500384, 'lon':134.624289}, + { 'name':'Ramallah', 'lat':31.9073509, 'lon':35.5354719}, + { 'name':'Panama City', 'lat':9.101179, 'lon':-79.402864}, + { 'name':'Port Moresby', 'lat':-9.4438, 'lon':147.180267}, + { 'name':'Asuncion', 'lat':-25.26374, 'lon':-57.575926}, + { 'name':'Lima', 'lat':-12.046374, 'lon':-77.042793}, + { 'name':'Manila', 'lat':14.599512, 'lon':120.98422}, + { 'name':'Adamstown', 'lat':-25.06629, 'lon':-130.100464}, + { 'name':'Warsaw', 'lat':52.229676, 'lon':21.012229}, + { 'name':'Lisbon', 'lat':38.722252, 'lon':-9.139337}, + { 'name':'San Juan', 'lat':18.466334, 'lon':-66.105722}, + { 'name':'Doha', 'lat':25.285447, 'lon':51.53104}, + { 'name':'Saint-Denis', 'lat':-20.882057, 'lon':55.450675}, + { 'name':'Bucharest', 'lat':44.426767, 'lon':26.102538}, + { 'name':'Moscow', 'lat':55.755826, 'lon':37.6173}, + { 'name':'Kigali', 'lat':-1.957875, 'lon':30.112735}, + { 'name':'St. Pierre', 'lat':46.775846, 'lon':-56.180636}, + { 'name':'Kingstown', 'lat':13.160025, 'lon':-61.224816}, + { 'name':'Apia', 'lat':-13.850696, 'lon':-171.751355}, + { 'name':'San Marino', 'lat':43.935591, 'lon':12.447281}, + { 'name':'São Tomé', 'lat':0.330192, 'lon':6.733343}, + { 'name':'Riyadh', 'lat':24.749403, 'lon':46.902838}, + { 'name':'Dakar', 'lat':14.764504, 'lon':-17.366029}, + { 'name':'Belgrade', 'lat':44.786568, 'lon':20.448922}, + { 'name':'Victoria', 'lat':-4.619143, 'lon':55.451315}, + { 'name':'Freetown', 'lat':8.465677, 'lon':-13.231722}, + { 'name':'Singapore', 'lat':1.280095, 'lon':103.850949}, + { 'name':'Bratislava', 'lat':48.145892, 'lon':17.107137}, + { 'name':'Ljubljana', 'lat':46.056947, 'lon':14.505751}, + { 'name':'Honiara', 'lat':-9.445638, 'lon':159.9729}, + { 'name':'Mogadishu', 'lat':2.046934, 'lon':45.318162}, + { 'name':'Pretoria', 'lat':-25.747868, 'lon':28.229271}, + { 'name':'King Edward Point', 'lat':-54.28325, 'lon':-36.493735}, + { 'name':'Seoul', 'lat':37.566535, 'lon':126.977969}, + { 'name':'Tskhinvali', 'lat':42.22146, 'lon':43.964405}, + { 'name':'Juba', 'lat':4.859363, 'lon':31.57125}, + { 'name':'Madrid', 'lat':40.416775, 'lon':-3.70379}, + { 'name':'Sri Jayawardenepura Kotte', 'lat':6.89407, 'lon':79.902478}, + { 'name':'Gustavia', 'lat':17.896435, 'lon':-62.852201}, + { 'name':'Basseterre', 'lat':17.302606, 'lon':-62.717692}, + { 'name':'Castries', 'lat':14.010109, 'lon':-60.987469}, + { 'name':'Marigot', 'lat':18.067519, 'lon':-63.082466}, + { 'name':'Khartoum', 'lat':15.500654, 'lon':32.559899}, + { 'name':'Paramaribo', 'lat':5.852036, 'lon':-55.203828}, + { 'name':'Longyearbyen ', 'lat':78.062, 'lon':22.055}, + { 'name':'Mbabane', 'lat':-26.305448, 'lon':31.136672}, + { 'name':'Stockholm', 'lat':59.329323, 'lon':18.068581}, + { 'name':'Bern', 'lat':46.947974, 'lon':7.447447}, + { 'name':'Damascus', 'lat':33.513807, 'lon':36.276528}, + { 'name':'Taipei', 'lat':25.032969, 'lon':121.565418}, + { 'name':'Dushanbe', 'lat':38.559772, 'lon':68.787038}, + { 'name':'Dodoma', 'lat':-6.162959, 'lon':35.751607}, + { 'name':'Bangkok', 'lat':13.756331, 'lon':100.501765}, + { 'name':'Dili', 'lat':-8.556856, 'lon':125.560314}, + { 'name':'Lomé', 'lat':6.172497, 'lon':1.231362}, + { 'name':'Nukunonu', 'lat':-9.2005, 'lon':-171.848}, + { 'name':'Nuku\'alofa', 'lat':-21.139342, 'lon':-175.204947}, + { 'name':'Tiraspol', 'lat':46.848185, 'lon':29.596805}, + { 'name':'Port of Spain', 'lat':10.654901, 'lon':-61.501926}, + { 'name':'Edinburgh of the Seven Seas', 'lat':-37.068042, 'lon':-12.311315}, + { 'name':'Tunis', 'lat':36.806495, 'lon':10.181532}, + { 'name':'Ankara', 'lat':39.933364, 'lon':32.859742}, + { 'name':'Ashgabat', 'lat':37.960077, 'lon':58.326063}, + { 'name':'Cockburn Town', 'lat':21.467458, 'lon':-71.13891}, + { 'name':'Funafuti', 'lat':-8.520066, 'lon':179.198128}, + { 'name':'Charlotte Amalie', 'lat':18.3419, 'lon':-64.930701}, + { 'name':'Kampala', 'lat':0.347596, 'lon':32.58252}, + { 'name':'Kiev', 'lat':50.4501, 'lon':30.5234}, + { 'name':'Abu Dhabi', 'lat':24.299174, 'lon':54.697277}, + { 'name':'London', 'lat':51.507351, 'lon':-0.127758}, + { 'name':'Washington', 'lat':38.907192, 'lon':-77.036871}, + { 'name':'Montevideo', 'lat':-34.901113, 'lon':-56.164531}, + { 'name':'Tashkent', 'lat':41.299496, 'lon':69.240073}, + { 'name':'Port Vila', 'lat':-17.733251, 'lon':168.327325}, + { 'name':'Vatican City', 'lat':41.902179, 'lon':12.453601}, + { 'name':'Caracas', 'lat':10.480594, 'lon':-66.903606}, + { 'name':'Hanoi', 'lat':21.027764, 'lon':105.83416}, + { 'name':'Mata-Utu', 'lat':-13.282509, 'lon':-176.176447}, + { 'name':'El Aaiún', 'lat':27.125287, 'lon':-13.1625}, + { 'name':'Sana\'a', 'lat':15.369445, 'lon':44.191007}, + { 'name':'Lusaka', 'lat':-15.387526, 'lon':28.322817}, + { 'name':'Harare', 'lat':-17.825166, 'lon':31.03351} ] \ No newline at end of file diff --git a/blueprints/test_chemistry/pp/plotting/maps/run.py b/blueprints/test_chemistry/pp/plotting/maps/run.py new file mode 100644 index 0000000000000000000000000000000000000000..c9820b99942b44bdd379d2617eaf091122fbc107 --- /dev/null +++ b/blueprints/test_chemistry/pp/plotting/maps/run.py @@ -0,0 +1,110 @@ +import datetime + +from argparse import ArgumentParser, ArgumentTypeError + +def valid_date(s): + try: + return datetime.datetime.strptime(s, "%Y-%m-%d") + except ValueError: + raise ArgumentTypeError("Not a valid date: '{0}'.".format(s)) + +parser = ArgumentParser() +parser.add_argument('wrf_data_fpath_pattern', type=str, help='WRF input data path pattern (including strptime pattern for date and {domain:02d} for domain (01, 02, ...))') +parser.add_argument('date', type=valid_date, help='Date to process (YYYY-mm-dd)') +parser.add_argument('interval_hours', type=int, help='Time interval in hours for plotting') +parser.add_argument('fcst_time_hours', type=int, help='Forecast time in hours') +parser.add_argument('plot_fpath_pattern', type=str, help='Plot output path pattern (including strptime pattern for date, {plot:s} for plot identifier (e.g., 500hPa...) and {domain:02d} for domain (01, 02)') + +args = parser.parse_args() + +from plots import plevel_map, sfc_map, slp_map +import pois +import netCDF4 +from wrf import getvar +import numpy as np +import cartopy.feature as cfeature +import contextily +import rasterio + +times = [ args.date + datetime.timedelta(hours=t) for t in range(0, args.fcst_time_hours, args.interval_hours) ] + +# for calculating delta precip +prec_prev = None + +def bg_function(ax, wrfcrs): + wrfcrs_proj4 = wrfcrs.proj4_init + wrfcrs_proj4 = wrfcrs_proj4.replace('+nadgrids=@null', '') + rasteriocrs = rasterio.crs.CRS.from_proj4(wrfcrs_proj4) + contextily.add_basemap(ax, source=contextily.providers.Stamen.TonerLite, crs=rasteriocrs, zorder=0) + +for t in times: + input_path_d01 = t.strftime(args.wrf_data_fpath_pattern).format(domain=1) + input_path_d02 = t.strftime(args.wrf_data_fpath_pattern).format(domain=2) + + try: + nc_d01 = netCDF4.Dataset(input_path_d01) + nc_d02 = netCDF4.Dataset(input_path_d02) + except Exception as e: + print("Not all data available for step {:s}: {:s}".format(t.strftime("%Y-%m-%d %H"), str(e))) + continue + + print(t) + + # 500 hPa Europe domain +# output_path = t.strftime(args.plot_fpath_pattern).format(plot='500hPa', domain=1) +# plevel_map(t, nc_d01, output_path, 500) + + # Surface Europe domain + # temperature +# sfc_map(t, nc_d01, lambda nc: getvar(nc, 'T2') - 273.15, +# levels=[-25, -20, -15, -12.5, -10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30, 32.5, 35], +# ticks=[-25, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35], +# label="Temperature (2m) [°C]", pois=pois.capitals, +# output_path=t.strftime(args.plot_fpath_pattern).format(plot='T', domain=1)) + + # Ozone + for domain, nc_d, poilist in zip([1,2],[nc_d01, nc_d02], [pois.capitals, pois.southern_germany]): + sfc_map(t, nc_d, lambda nc: getvar(nc, 'o3')[0,:,:] * 1e3, + levels=range(0, 160, 10), + ticks=range(0, 150, 20), + label="Ozone (ground level) [ppbv]", bg_function=bg_function, + output_path=t.strftime(args.plot_fpath_pattern).format(plot='O3', domain=domain), + extend='max') + # NO2 + sfc_map(t, nc_d, lambda nc: getvar(nc, 'no2')[0,:,:] * 1e3, + levels=range(0, 110, 10), + ticks=range(0, 110, 20), + label="NO$_2$ (ground level) [ppbv]", bg_function=bg_function, + output_path=t.strftime(args.plot_fpath_pattern).format(plot='NO2', domain=domain), + extend='max') + # PM + def get_pm(nc, bins=[ 1, 2, 3 ]) : + pmspecs = [ "so4", "no3", "smpa", "smpbb", "glysoa_sfc", "biog1_c", "biog1_o", "cl", "co3", "nh4", "na", "cl", "oin", "oc", "bc", "water" ] # ug kg-1 + alt = getvar(nc, "ALT")[0,:,:] # inverse density m3 kg-1 + all = [ getvar(nc, spec + "_a{:02d}".format(bin))[0,:,:] / alt for spec in pmspecs for bin in bins ] + return np.sum(all, axis=0) + +# sfc_map(t, nc_d, lambda nc: get_pm(nc, bins=[1]), +# levels=range(0, 130, 10), +# ticks=range(0, 130, 20), +# label="UFP (ground level) [$\mu$g m$^{{-3}}$]", bg_function=bg_function, +# output_path=t.strftime(args.plot_fpath_pattern).format(plot='UFP', domain=domain), +# extend='max') + sfc_map(t, nc_d, lambda nc: get_pm(nc, bins=[1,2,3]), + levels=range(0, 130, 10), + ticks=range(0, 130, 20), + label="PM$_{{2.5}}$ (ground level) [$\mu$g m$^{{-3}}$]", bg_function=bg_function, + output_path=t.strftime(args.plot_fpath_pattern).format(plot='PM2_5', domain=domain), + extend='max') + sfc_map(t, nc_d, lambda nc: get_pm(nc, bins=[1,2,3,4]), + levels=range(0, 130, 10), + ticks=range(0, 130, 20), + label="PM$_{{10}}$ (ground level) [$\mu$g m$^{{-3}}$]", bg_function=bg_function, + output_path=t.strftime(args.plot_fpath_pattern).format(plot='PM10', domain=domain), + extend='max') + + # SLP Europe +# output_path = t.strftime(args.plot_fpath_pattern).format(plot='SLP', domain=1) +# prec_prev = slp_map(t, nc_d01, prec_prev, output_path) + + nc_d01.close() diff --git a/blueprints/test_chemistry/pp/plotting/maps/run.sh b/blueprints/test_chemistry/pp/plotting/maps/run.sh index 65688d8fd9d18fa577bb1984d581bd2d81ad39cb..248dc5b5b142ce64455a10d6f2baef35fed36e73 100644 --- a/blueprints/test_chemistry/pp/plotting/maps/run.sh +++ b/blueprints/test_chemistry/pp/plotting/maps/run.sh @@ -1,23 +1,22 @@ #!/bin/bash -l -#SBATCH --partition=alcc1 -#SBATCH -o /alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/maps.%j.%N.out +#SBATCH --partition=alcc1,alcc3,epyc +#SBATCH -o /alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/maps/%j.%N.out #SBATCH -J maps -##SBATCH --nodelist=alcc131 #SBATCH --ntasks=1 -#SBATCH --mem=5G +#SBATCH --mem=10G #SBATCH --mail-type=FAIL #SBATCH --mail-user=christoph.knote@med.uni-augsburg.de -#SBATCH --time=00:30:00 +#SBATCH --time=01:00:00 module load gnu geos ffmpeg export MPLBACKEND="agg" -scriptPath=/alcc/gpfs2/home/u/knotechr/wrfotron/blueprints/operational_chemistry/pp/plotting/maps/plot.py +scriptPath=/alcc/gpfs2/home/u/knotechr/wrfotron/blueprints/operational_chemistry/pp/plotting/maps/run.py wrfDataPathPattern="/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout_d{domain:02d}_%Y-%m-%d_%H:%M:%S" -tmpWorkDir="/alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/fcst_movie" -plotPathPattern="${tmpWorkDir}/{plot:s}_%Y%m%d%H.png" -webDir="/alcc/gpfs2/scratch/mbees/knotechr/plots/WRF/operational_chemistry/maps" +tmpWorkDir="/alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/maps/movie" +plotPathPattern="${tmpWorkDir}/{plot:s}_d{domain:02d}_%Y%m%d%H.png" +webDir="/alcc/gpfs2/scratch/mbees/knotechr/plots/WRF/operational_chemistry/__domain__/maps/%Y/%m" dateToProcess="today" fcstIntervalHours=1 fcstRangeHours=72 @@ -28,45 +27,24 @@ rm -f ${tmpWorkDir}/* python3 $scriptPath ${wrfDataPathPattern} $(date -u --date="${dateToProcess}" "+%Y-%m-%d") ${fcstIntervalHours} ${fcstRangeHours} ${plotPathPattern} -# once ffmpeg gets updated (to something > 2 Nov 2018), we can add still frames at the end with this -vf filter: -#tpad=stop_mode=clone:stop_duration=2 - -# see here for why we need vf argument https://stackoverflow.com/questions/20847674/ffmpeg-libx264-height-not-divisible-by-2 frameRate=3 -# get all plot type names ("500hpa_Europe, o3_Europe", ...) -typs=$(ls -1 ${tmpWorkDir}/*.png | sed -e "s|${tmpWorkDir}/||" | sed -E -e "s/_[0-9]+.png//g" | uniq) -for typ in ${typs} + +for domain in d01 d02 do - # super ugly hack: copy last image seconds*framerate times to another name that ensures they are at the end (using big number as date here) - # to make video "hold" on last image for 'seconds' - allFiles=( $(ls -r1 ${tmpWorkDir}/${typ}_*.png) ) - lastFile=${allFiles[0]} - for iter in $(seq 0 ${nSecondsWait}) + # get all plot type names ("500hpa_d01, o3_d02", ...) + typs=$(ls -1 ${tmpWorkDir}/*${domain}*.png | sed -e "s|${tmpWorkDir}/||" | sed -E -e "s/_[0-9]+.png//g" | uniq) + + for typ in ${typs} do - for frm in $(seq 0 ${frameRate}) - do - let bignum='300000000000+iter*frm' - cp ${lastFile} ${tmpWorkDir}/${typ}_${bignum}.png - done + echo "" + echo $typ $domain + outFilePath=$(date --date="${dateToProcess} 0 UTC" "+${webDir/__domain__/${domain}}/${typ/_${domain}/}_%Y%m%d.mp4") + mkdir -p $(dirname $outFilePath) + rm -f ${outFilePath} + # see here for why we need vf argument https://stackoverflow.com/questions/20847674/ffmpeg-libx264-height-not-divisible-by-2 + cat ${tmpWorkDir}/${typ}_*.png | ffmpeg -y -r ${frameRate} -f image2pipe -i - -vcodec libx264 -crf 18 -preset slow -pix_fmt yuv420p -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2, tpad=stop_mode=clone:stop_duration=2, fps=fps=${frameRate}" ${outFilePath} done - outFilePath=$(date --date="${dateToProcess} 0 UTC" "+${webDir}/%Y%m%d_${typ}.mp4") - ffmpeg -y -r ${frameRate} -pattern_type glob -i "${tmpWorkDir}/${typ}_*.png" -vcodec libx264 -crf 18 -preset slow -pix_fmt yuv420p -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2" ${outFilePath} > /dev/null done +rm -f ${tmpWorkDir}/* -for anHour in $(seq 0 ${fcstIntervalHours} ${fcstRangeHours}) -do - # get path to plots which we created in tmpdir - plotPath=$(date -u --date="${dateToProcess} 0 UTC +${anHour} hours" "+${plotPathPattern}") - # replace "type" placeholder with wildcard to match all - allPlotPath=${plotPath/\{plot:s\}/\*} - - # where to put them on the web - plotWebDir=$(date -u --date="${dateToProcess} 0 UTC +${anHour} hours" "+${webDir}/%Y/%m/%d") - mkdir -p $plotWebDir - - for aPlotPath in ${allPlotPath} - do - mv ${aPlotPath} ${plotWebDir} - done -done diff --git a/blueprints/test_chemistry/pp/plotting/maps/tools.py b/blueprints/test_chemistry/pp/plotting/maps/tools.py new file mode 100644 index 0000000000000000000000000000000000000000..b2a7fdf244453f5afd16089d3842b9c68641624d --- /dev/null +++ b/blueprints/test_chemistry/pp/plotting/maps/tools.py @@ -0,0 +1,68 @@ +from wrf import cartopy_xlim, cartopy_ylim +import matplotlib.gridspec as gridspec + +from matplotlib.offsetbox import AnchoredText + +import datetime +import numpy as np + +t2datetime = lambda x: datetime.datetime.strptime(x, '%Y-%m-%d_%H:%M:%S') + +def wrf_time(nc): + first = ''.join([x.decode('utf-8') for x in nc['Times'][0]]) + return t2datetime(first) + +def total_prec(nc, prec_p): + prec_vars = ['RAINC', 'RAINSH', 'RAINNC', 'SNOWNC', 'GRAUPELNC', 'HAILNC'] + start = t2datetime(nc.START_DATE) + prec = np.sum([nc[var][0, :, :] for var in prec_vars], axis=0) +# + if type(prec_p) == type(None): + prec_p = prec +# + if start == wrf_time(nc): + return prec + else: + return prec-prec_p + +# buffer 0.1 -> 10% border around the edges +def cartopy_lim_buffered(field, direction, buffer=0.1): + if direction == "x": + lims = cartopy_xlim(field) + elif direction == "y": + lims = cartopy_ylim(field) + else: + raise Exception("Direction not understood!") + dlim = lims[1] - lims[0] + lims[0] += 0.5 * buffer * dlim + lims[1] -= 0.5 * buffer * dlim + return lims + +def make_fig_and_ax(plt, crs, add_horiz_cbar=False): + fig = plt.figure(figsize=(7,7)) +# + gs = gridspec.GridSpec(2, 2, width_ratios=[20,1], height_ratios=[20,1], wspace=0.05, hspace=0.05) +# + ax = fig.add_subplot(gs[0,0], projection=crs) + ax_cbarv = fig.add_subplot(gs[0,1]) + ax_cbarh = None + if add_horiz_cbar: + ax_cbarh = fig.add_subplot(gs[1,0]) + return fig, ax, ax_cbarv, ax_cbarh + +def add_title_and_description(nc, title, ax): + ''' + title: including {init:s} for timestamp of run initialization, and {fcsthour:d} for forecast hour + ''' + fcstHour = int((wrf_time(nc) - t2datetime(nc.START_DATE)).total_seconds()/3600.) +# + ax.set_title(title.format(init=t2datetime(nc.START_DATE).strftime("%d.%m.%Y %H:%M UTC"), + fcsthour=fcstHour), fontdict={ 'fontsize':8, 'horizontalalignment': 'left'}, loc='left') +# + texttime = AnchoredText(wrf_time(nc).strftime("%H UTC %a"), + loc=2, pad=0.2, prop={'size': 10}, frameon=True) + ax.add_artist(texttime) +# + textcopyright = AnchoredText("$\copyright$ MBEES, Faculty of Medicine, University of Augsburg", + loc=4, pad=0.1, prop={'size': 6, 'color':'white'}, frameon=False) + ax.add_artist(textcopyright) diff --git a/blueprints/test_chemistry/pp/plotting/meteograms/plot.py b/blueprints/test_chemistry/pp/plotting/meteograms/plot.py index 87d7088bbb95c754a70b97d6a393a71ea22d3c09..278d373363974298fcbb75be02b199cc2c143f57 100755 --- a/blueprints/test_chemistry/pp/plotting/meteograms/plot.py +++ b/blueprints/test_chemistry/pp/plotting/meteograms/plot.py @@ -3,30 +3,35 @@ from wrf import getvar, ll_to_xy, to_np, extract_times import numpy as np + def fT2m(D, xmin, xmax, ymin, ymax): - return D['T2'][0, xmin:xmax, ymin:ymax] - 273.15 + return D["T2"][0, xmin:xmax, ymin:ymax] - 273.15 + def fTd2m(D, xmin, xmax, ymin, ymax): return getvar(D, "td2")[xmin:xmax, ymin:ymax] + def wind10m(D, lat, lon): """Returns the wind speed and the wind direction at 10 m""" try: x, y = ll_to_xy(D, lat, lon) - w_speed, w_dir = getvar(D, 'uvmet10_wspd_wdir') - return w_speed.values[int(x),int(y)], w_dir.values[int(x),int(y)] + w_speed, w_dir = getvar(D, "uvmet10_wspd_wdir") + return w_speed.values[int(x), int(y)], w_dir.values[int(x), int(y)] except: return np.NaN, np.NaN + def wind_uv(D, lat, lon): """Returns the u and v wind components at 10 m""" try: x, y = ll_to_xy(D, lat, lon) - w_u, w_v = getvar(D, 'uvmet10') + w_u, w_v = getvar(D, "uvmet10") return w_u.values[int(x), int(y)], w_v.values[int(x), int(y)] except: return np.NaN, np.NaN + def Mslp(D, lat, lon): try: """Returns the mean sea level pressure""" @@ -36,6 +41,7 @@ def Mslp(D, lat, lon): except: return np.NaN + def cin_cape(D, lat, lon): try: """Returns the 2D max cin and cape""" @@ -45,59 +51,108 @@ def cin_cape(D, lat, lon): except: return np.NaN, np.NaN + def vert_wind(D, lat, lon): try: from wrf import vertcross, CoordPair + x, y = ll_to_xy(D, lat, lon) xy = CoordPair(x=x, y=y) - u, v = getvar(D, 'uvmet') - p = getvar(D, "p", units='hPa') + u, v = getvar(D, "uvmet") + p = getvar(D, "p", units="hPa") levels = np.arange(400, 1060, 60) - u_vert = vertcross(u, p, levels=levels, start_point=xy, end_point=CoordPair(x=x+1, y=y+1), meta=False, missing=np.NaN) - v_vert = vertcross(v, p, levels=levels, start_point=xy, end_point=CoordPair(x=x+1, y=y+1), meta=False, missing=np.NaN) - return list(u_vert[:,0]), list(v_vert[:,0]) + u_vert = vertcross( + u, + p, + levels=levels, + start_point=xy, + end_point=CoordPair(x=x + 1, y=y + 1), + meta=False, + missing=np.NaN, + ) + v_vert = vertcross( + v, + p, + levels=levels, + start_point=xy, + end_point=CoordPair(x=x + 1, y=y + 1), + meta=False, + missing=np.NaN, + ) + return list(u_vert[:, 0]), list(v_vert[:, 0]) except: levels = np.arange(400, 1060, 60) - return [np.NaN for x in range(len(levels))], [np.NaN for x in range(len(levels))] + return [np.NaN for x in range(len(levels))], [ + np.NaN for x in range(len(levels)) + ] + def vert_T(D, lat, lon): try: from wrf import vertcross, CoordPair + x, y = ll_to_xy(D, lat, lon) xy = CoordPair(x=x, y=y) - T = getvar(D, 'tc') - p = getvar(D, "p", units='hPa') + T = getvar(D, "tc") + p = getvar(D, "p", units="hPa") levels = np.arange(400, 1060, 60) - T_vert = to_np(vertcross(T, p, levels=levels, start_point=xy, end_point=CoordPair(x=x + 1, y=y + 1), meta=False, missing=np.NaN)) - return T_vert[:,0][:] + T_vert = to_np( + vertcross( + T, + p, + levels=levels, + start_point=xy, + end_point=CoordPair(x=x + 1, y=y + 1), + meta=False, + missing=np.NaN, + ) + ) + return T_vert[:, 0][:] except: levels = np.arange(400, 1060, 60) - return [np.NaN for x in range(len(levels))], [np.NaN for x in range(len(levels))] + return [np.NaN for x in range(len(levels))], [ + np.NaN for x in range(len(levels)) + ] def fcfrac(D, lat, lon): from wrf import interplevel, CoordPair + levels = np.arange(400, 1030, 30) try: x, y = ll_to_xy(D, lat, lon) - p = getvar(D, "p", units='hPa') - cf = interplevel(D['CLDFRA'], p, levels) + p = getvar(D, "p", units="hPa") + cf = interplevel(D["CLDFRA"], p, levels) return cf.values[:, int(x), int(y)] - + except: return [np.NaN for x in range(len(levels))] def ftot_accum_rain_2d(D, xmin, xmax, ymin, ymax): try: - return np.sum([D[var][0, xmin:xmax, ymin:ymax] for var in - ['RAINC', 'RAINSH', 'RAINNC', - 'SNOWNC', 'GRAUPELNC', 'HAILNC']], axis=0) + return np.sum( + [ + D[var][0, xmin:xmax, ymin:ymax] + for var in [ + "RAINC", + "RAINSH", + "RAINNC", + "SNOWNC", + "GRAUPELNC", + "HAILNC", + ] + ], + axis=0, + ) except: return np.mgrid[xmin:xmax, ymin:ymax][0] * np.NaN + def fdrain(DS, xmin, xmax, ymin, ymax): - tot_accum_rain = np.array([ftot_accum_rain_2d(D, xmin, xmax, ymin, ymax) for D in DS]) + tot_accum_rain = np.array( + [ftot_accum_rain_2d(D, xmin, xmax, ymin, ymax) for D in DS] + ) # weird construct to remove negative rain rates # because timeseries have restarts in them @@ -107,143 +162,180 @@ def fdrain(DS, xmin, xmax, ymin, ymax): dt = np.array([x.total_seconds() for x in (starts[1:] - starts[:-1])]) drain = np.zeros_like(tot_accum_rain) drain[1:] = tot_accum_rain[1:] - tot_accum_rain[:-1] - drain[1:] = np.where(dt > 0, tot_accum_rain[1:].T, (tot_accum_rain[1:] - tot_accum_rain[:-1]).T).T + drain[1:] = np.where( + dt > 0, tot_accum_rain[1:].T, (tot_accum_rain[1:] - tot_accum_rain[:-1]).T + ).T return drain + def ftot_accum_snow_2d(D, xmin, xmax, ymin, ymax): try: - return np.sum([D[var][0, xmin:xmax, ymin:ymax] for var in - ['SNOWNC', 'GRAUPELNC']], axis=0) + return np.sum( + [D[var][0, xmin:xmax, ymin:ymax] for var in ["SNOWNC", "GRAUPELNC"]], axis=0 + ) except: return np.mgrid[xmin:xmax, ymin:ymax][0] * np.NaN - + + def fdsnow(DS, xmin, xmax, ymin, ymax): - tot_accum_snow = np.array([ftot_accum_snow_2d(D, xmin, xmax, ymin, ymax) for D in DS]) + tot_accum_snow = np.array( + [ftot_accum_snow_2d(D, xmin, xmax, ymin, ymax) for D in DS] + ) starts = np.array([t2datetime(D.START_DATE) for D in DS]) dt = np.array([x.total_seconds() for x in (starts[1:] - starts[:-1])]) drain = np.zeros_like(tot_accum_snow) drain[1:] = tot_accum_snow[1:] - tot_accum_snow[:-1] - drain[1:] = np.where(dt > 0, tot_accum_snow[1:].T, (tot_accum_snow[1:] - tot_accum_snow[:-1]).T).T + drain[1:] = np.where( + dt > 0, tot_accum_snow[1:].T, (tot_accum_snow[1:] - tot_accum_snow[:-1]).T + ).T return drain + def t2datetime(x): import datetime - return datetime.datetime.strptime(x, '%Y-%m-%d_%H:%M:%S') + + return datetime.datetime.strptime(x, "%Y-%m-%d_%H:%M:%S") + def ftimes(D): - """ Convert first WRF timestamp description into datetimes string""" - return t2datetime(''.join([x.decode('utf-8') for x in D['Times'][0]])) + """Convert first WRF timestamp description into datetimes string""" + return t2datetime("".join([x.decode("utf-8") for x in D["Times"][0]])) + def corresponding_pixels(D, location_lon, location_lat, deg_round_loc): import numpy as np - xlon = D['XLONG'][0] - xlat = D['XLAT'][0] - around_loc = ((np.abs(xlat - location_lat) < deg_round_loc) & (np.abs(xlon - location_lon) < deg_round_loc)) - + xlon = D["XLONG"][0] + xlat = D["XLAT"][0] + + around_loc = (np.abs(xlat - location_lat) < deg_round_loc) & ( + np.abs(xlon - location_lon) < deg_round_loc + ) + xmin, xmax = np.nonzero(around_loc)[0][[0, -1]] ymin, ymax = np.nonzero(around_loc)[1][[0, -1]] xlon = xlon[xmin:xmax, ymin:ymax] xlat = xlat[xmin:xmax, ymin:ymax] - #print("Using {} pixels for plot".format(np.size(xlon))) + # print("Using {} pixels for plot".format(np.size(xlon))) return (xmin, xmax), (ymin, ymax), xlon, xlat + def haversine(lat1, lon1, lat2, lon2): """ Compute distance between lat lon points """ import numpy as np + R = 6378137 # Earth radius dLat = np.deg2rad(lat2) - np.deg2rad(lat1) dLon = np.deg2rad(lon2) - np.deg2rad(lon1) - a = np.sin(dLat / 2) * np.sin(dLat / 2) + np.cos(np.deg2rad(lat1)) * np.cos(np.deg2rad(lat2)) * np.sin( - dLon / 2.) * np.sin(dLon / 2.) - c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1. - a)); - return R * c; + a = np.sin(dLat / 2) * np.sin(dLat / 2) + np.cos(np.deg2rad(lat1)) * np.cos( + np.deg2rad(lat2) + ) * np.sin(dLon / 2.0) * np.sin(dLon / 2.0) + c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1.0 - a)) + return R * c + def compute_average_resolution(D, xmin, xmax, ymin, ymax): import numpy as np - lon = D['XLONG'][0] - lat = D['XLAT'][0] + + lon = D["XLONG"][0] + lat = D["XLAT"][0] dy_avg = np.average( - haversine(lat[xmin:xmax, ymin:ymax - 1], lon[xmin:xmax, ymin:ymax - 1], lat[xmin:xmax, ymin + 1:ymax], - lon[xmin:xmax, ymin + 1:ymax])) + haversine( + lat[xmin:xmax, ymin : ymax - 1], + lon[xmin:xmax, ymin : ymax - 1], + lat[xmin:xmax, ymin + 1 : ymax], + lon[xmin:xmax, ymin + 1 : ymax], + ) + ) dx_avg = np.average( - haversine(lat[xmin:xmax - 1, ymin:ymax], lon[xmin:xmax - 1, ymin:ymax], lat[xmin + 1:xmax, ymin:ymax], - lon[xmin + 1:xmax, ymin:ymax])) + haversine( + lat[xmin : xmax - 1, ymin:ymax], + lon[xmin : xmax - 1, ymin:ymax], + lat[xmin + 1 : xmax, ymin:ymax], + lon[xmin + 1 : xmax, ymin:ymax], + ) + ) return (dy_avg + dx_avg) / 2 + def datetime_to_float(d): import datetime + return d.timestamp() # -------------------------------------------------------------------------------- -def plot_var(filenames, outfile, deg_round_loc=None, - location_lon=11.5754900, location_lat=48.1374300): + +def plot_var( + filenames, outfile, deg_round_loc=0.2, location_lon=10.8978, location_lat=48.3705 +): """ Plot meteogram for given location (location_lon,location_lat) - Default is Munich Script uses all pixels around location in a radius of ''deg_round_loc'' to attain model spread. """ import matplotlib as mpl - mpl.use('Agg') + + mpl.use("Agg") import matplotlib.pyplot as plt import numpy as np import netCDF4 - import glob import datetime - from matplotlib.dates import DateFormatter, HourLocator - from matplotlib.dates import date2num - from matplotlib.dates import DateFormatter, HourLocator - from matplotlib.dates import date2num from matplotlib.colors import LinearSegmentedColormap from matplotlib.font_manager import FontProperties - def draw_vlines(ax, labels = False): + def draw_vlines(ax, labels=False): for t in time[1:]: - if t.hour == 0. and t.minute == 0.: - ax.vlines(datetime_to_float(t), *ax.get_ylim(), color='black', linestyle=':', lw=1) - for t in time: - if t.hour == 1. and t.minute == 0.: + if t.hour == 0.0 and t.minute == 0.0: + ax.vlines( + datetime_to_float(t), + *ax.get_ylim(), + color="black", + linestyle=":", + lw=1 + ) + for t in time: + if t.hour == 1.0 and t.minute == 0.0: if labels: ycoo = ax.get_ylim() - ax.text(datetime_to_float(t+datetime.timedelta(hours = 11)), -ycoo[1]/3, t.strftime("%A %-d %b"), ha='center') + ax.text( + datetime_to_float(t + datetime.timedelta(hours=11)), + -ycoo[1] / 3, + t.strftime("%A %-d %b"), + ha="center", + ) # -----------------reading input files:------------------------- def get_only_valid_file(f): try: nc = netCDF4.Dataset(f) - nul = extract_times(nc, None) + ts = extract_times(nc, None) + if len(ts) == 0: + nc = None except: nc = None return nc - - DS = [ get_only_valid_file(f) for f in sorted(filenames) ] - DS = [ x for x in DS if not x is None ] - -# import pdb; pdb.set_trace() - if deg_round_loc is None: - # approximate round loc by taking 3 times rough horizontal resolution in degrees - deg_round_loc = 3.0 * DS[0].DX/40000000*360.0 - print("Assume {:5.2f} degrees patch for stats calculation".format(deg_round_loc)) + DS = [get_only_valid_file(f) for f in sorted(filenames)] + DS = [x for x in DS if not x is None] - (xmin, xmax), (ymin, ymax), xlon, xlat = corresponding_pixels(DS[0], location_lon, location_lat, deg_round_loc) + (xmin, xmax), (ymin, ymax), xlon, xlat = corresponding_pixels( + DS[0], location_lon, location_lat, deg_round_loc + ) time = np.array([ftimes(D) for D in DS]) extent = (xmin, xmax, ymin, ymax) - x, y = ll_to_xy(DS[0], location_lat, location_lon) # -----------------reading data from input files----------------------- T2m = [fT2m(D, *extent) for D in DS] @@ -259,17 +351,14 @@ def plot_var(filenames, outfile, deg_round_loc=None, drain = fdrain(DS, *extent) ravg = np.average(drain, axis=(1, 2)) - + dsnow = fdsnow(DS, *extent) snow = np.average(dsnow, axis=(1, 2)) - #snow = [fSW(D, *extent) for D in DS] - #snow = np.average(snow, axis=(1, 2)) - conv = [cin_cape(D, location_lat, location_lon) for D in DS] cape = [x[1] for x in conv] cin = [x[0] for x in conv] - + w_vert = [vert_wind(D, location_lat, location_lon) for D in DS] u_vert = [x[0] for x in w_vert] v_vert = [x[1] for x in w_vert] @@ -279,57 +368,72 @@ def plot_var(filenames, outfile, deg_round_loc=None, cf_vert = [list(fcfrac(D, location_lat, location_lon)) for D in DS] # figure settings: - fig = plt.figure(figsize=(8.27*1.25, 11.69*1.25), num=1) + fig = plt.figure(figsize=(8.27 * 1.25, 11.69 * 1.25), num=1) plt.clf() - mplotkw = {'ls': '-', 'marker': '.'} + mplotkw = {"ls": "-", "marker": "."} + + ax = plt.subplot2grid((8, 20), (0, 19), rowspan=3) + ax0 = plt.subplot2grid((8, 20), (0, 0), colspan=19, rowspan=3) + ax1 = plt.subplot2grid((8, 20), (3, 0), colspan=19, rowspan=2) + ax2 = plt.subplot2grid((8, 20), (5, 0), colspan=19) + ax3 = plt.subplot2grid((8, 20), (6, 0), colspan=19) + ax4 = plt.subplot2grid((8, 20), (7, 0), colspan=19) - ax = plt.subplot2grid((8, 20), (0,19), rowspan=3) - ax0 = plt.subplot2grid((8, 20), (0,0), colspan=19, rowspan=3) - ax1 = plt.subplot2grid((8, 20), (3,0), colspan=19, rowspan=2) - ax2 = plt.subplot2grid((8, 20), (5,0), colspan=19) - ax3 = plt.subplot2grid((8, 20), (6,0), colspan=19) - ax4 = plt.subplot2grid((8, 20), (7,0), colspan=19) - # for the xticks: - hours = np.array(list(map(lambda x:x.hour, time))) + hours = np.array(list(map(lambda x: x.hour, time))) ind = list(np.where(hours % 6 == 0)[0]) # 0.--------------vertical plot------------------------ - - ax0.grid(axis='y') + + ax0.grid(axis="y") # (vertical) cloud fraction levels = np.arange(400, 1030, 30) cf_vert = np.matrix(cf_vert).transpose() cmap = LinearSegmentedColormap.from_list("", ["steelblue", "white"]) - time_mod = list(map(lambda x:datetime_to_float(x), time)) - cs = ax0.contourf(time_mod, levels, cf_vert, levels=[x/8 for x in range(10)], cmap=cmap) - ax0.contourf(time_mod, levels, cf_vert, levels=[x/8 for x in range(9)], cmap=cmap) - cbar = fig.colorbar(cs, cax=ax, orientation='vertical') + time_mod = list(map(lambda x: datetime_to_float(x), time)) + cs = ax0.contourf( + time_mod, levels, cf_vert, levels=[x / 8 for x in range(10)], cmap=cmap + ) + ax0.contourf(time_mod, levels, cf_vert, levels=[x / 8 for x in range(9)], cmap=cmap) + cbar = fig.colorbar(cs, cax=ax, orientation="vertical") cbar.set_ticks([x / 8 for x in range(0, 9)]) try: - cbar.ax.set_yticklabels(['{}/8'.format(x) for x in range(0,9)]) + cbar.ax.set_yticklabels(["{}/8".format(x) for x in range(0, 9)]) except: print("Plotting tick labels failed") - ax.set_ylabel('cloud cover') + ax.set_ylabel("cloud cover") # (vertical) temperature contours levels = np.arange(400, 1060, 60) T_vert = np.matrix(T_vert).transpose() - cs = ax0.contour(time_mod, levels, T_vert, cmap='plasma') - ax0.clabel(cs, fmt='%1.0f', fontsize=10) + cs = ax0.contour(time_mod, levels, T_vert, cmap="plasma") + ax0.clabel(cs, fmt="%1.0f", fontsize=10) # (vertical) wind barbs: for i, j in enumerate(time_mod[:]): - if i % (len(time)//24) == 0: + if i % (len(time) // 24) == 0: for k, l in enumerate(levels): - ax0.barbs(j, l, u_vert[i][k], v_vert[i][k], length=5, pivot='middle', linewidth=0.5) + ax0.barbs( + j, + l, + u_vert[i][k], + v_vert[i][k], + length=5, + pivot="middle", + linewidth=0.5, + ) ax0.set_yticks(np.arange(400, 1050, 50)) ax0.set_ylim(980, 400) - ax0.set_xlim(datetime_to_float(time[0] - datetime.timedelta(hours=3*len(time)/(24*3))), datetime_to_float(time[-1] + datetime.timedelta(hours=3*len(time)/(24*3)))) - ax0.set_ylabel('pressure [hPa]') + ax0.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + ax0.set_ylabel("pressure [hPa]") try: ax0.set_xticklabels([], minor=False) ax0.set_xticklabels([], minor=True) @@ -340,12 +444,11 @@ def plot_var(filenames, outfile, deg_round_loc=None, # 1.--------------temperature plot------------------------ - - ax1.fill_between(time_mod, Tdmin, Tdmax, alpha=.3, color='lightblue') - ax1.fill_between(time_mod, Tmin, Tmax, alpha=.3, color='grey') - ax1.grid(axis='y') - ax1.plot(time_mod, Tdavg, color='blue', label='Td(2m)', **mplotkw) - ax1.plot(time_mod, Tavg, color='black', label='T(2m)', **mplotkw) + ax1.fill_between(time_mod, Tdmin, Tdmax, alpha=0.3, color="lightblue") + ax1.fill_between(time_mod, Tmin, Tmax, alpha=0.3, color="grey") + ax1.grid(axis="y") + ax1.plot(time_mod, Tdavg, color="blue", label="Td(2m)", **mplotkw) + ax1.plot(time_mod, Tavg, color="black", label="T(2m)", **mplotkw) actual_day = time[0].day actual_month = time[0].month @@ -354,26 +457,46 @@ def plot_var(filenames, outfile, deg_round_loc=None, # min and max temperature for the current day for i, j in enumerate(time): - if time[i-1].month != actual_month and (i != 0): + if time[i - 1].month != actual_month and (i != 0): day += 1 - actual_month = time[i-1].month + actual_month = time[i - 1].month - elif (day != time[i-1].day-actual_day) and (i != 0) and (time[i-1].day-actual_day >=0): + elif ( + (day != time[i - 1].day - actual_day) + and (i != 0) + and (time[i - 1].day - actual_day >= 0) + ): day += 1 try: - if Tavg[i] == min(Tavg[day*24:(day+1)*24]): - ax1.text(datetime_to_float(time[i]), Tavg[i] +2, str(np.round(Tavg[i], 1)), color='blue', ha='center') - elif Tavg[i] == max(Tavg[day*24:(day+1)*24]): - ax1.text(datetime_to_float(time[i]), Tavg[i] + 1.5, str(np.round(Tavg[i], 1)), color='red', ha='center') + if Tavg[i] == min(Tavg[day * 24 : (day + 1) * 24]): + ax1.text( + datetime_to_float(time[i]), + Tavg[i] + 2, + str(np.round(Tavg[i], 1)), + color="blue", + ha="center", + ) + elif Tavg[i] == max(Tavg[day * 24 : (day + 1) * 24]): + ax1.text( + datetime_to_float(time[i]), + Tavg[i] + 1.5, + str(np.round(Tavg[i], 1)), + color="red", + ha="center", + ) except: pass - ax1.set_ylabel(r'temp. [$^\circ$C]') - ax1.legend(loc='best', fontsize='small', ncol=3) - ax1.set_ylim(min(min(Tmin), min(Tdmin))-3, - max(max(Tmax), max(Tdmax))+7.5) - ax1.set_xlim(datetime_to_float(time[0] - datetime.timedelta(hours=3*len(time)/(24*3))), datetime_to_float(time[-1] + datetime.timedelta(hours=3*len(time)/(24*3)))) + ax1.set_ylabel(r"temp. [$^\circ$C]") + ax1.legend(loc="best", fontsize="small", ncol=3) + ax1.set_ylim(min(min(Tmin), min(Tdmin)) - 3, max(max(Tmax), max(Tdmax)) + 7.5) + ax1.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) try: ax1.xaxis.set_ticklabels([], minor=True) ax1.xaxis.set_ticklabels([], minor=False) @@ -381,18 +504,37 @@ def plot_var(filenames, outfile, deg_round_loc=None, print("Plotting tick labels failed") ax1.set_xticks([time_mod[i] for i in ind]) draw_vlines(ax1) - + # 2.--------------10m wind plot------------------------ - - ax2.grid(axis='y') - ax2.plot(time_mod, wind_speed, color='black', **mplotkw) - ax2.set_xlim(datetime_to_float(time[0] - datetime.timedelta(hours=3*len(time)/(24*3))), datetime_to_float(time[-1] + datetime.timedelta(hours=3*len(time)/(24*3)))) + + ax2.grid(axis="y") + ax2.plot(time_mod, wind_speed, color="black", **mplotkw) + ax2.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) for i, j in enumerate(time_mod): - if i % (len(time)//24) == 0: - ax2.text(j, wind_speed[i] - 4.9, str(int(round(wind_speed[i]))), color='black', ha='center') - ax2.barbs(j, max(15, max(wind_speed)+5), wind_u[i], wind_v[i], length=5, pivot='middle', linewidth=0.5) - ax2.set_ylabel(r'10 m ws [ms$^{-1}$]') - ax2.set_ylim(-4.9, max(20, max(wind_speed)+8)) + if i % (len(time) // 24) == 0: + ax2.text( + j, + wind_speed[i] - 4.9, + str(int(round(wind_speed[i]))), + color="black", + ha="center", + ) + ax2.barbs( + j, + max(15, max(wind_speed) + 5), + wind_u[i], + wind_v[i], + length=5, + pivot="middle", + linewidth=0.5, + ) + ax2.set_ylabel(r"10 m ws [ms$^{-1}$]") + ax2.set_ylim(-4.9, max(20, max(wind_speed) + 8)) try: ax2.xaxis.set_ticklabels([], minor=True) ax2.xaxis.set_ticklabels([], minor=False) @@ -402,54 +544,94 @@ def plot_var(filenames, outfile, deg_round_loc=None, draw_vlines(ax2) # 3.--------------mslp plot------------------------ - - ax3.grid(axis='y') - ax3.plot(time_mod, mslp, color='red', **mplotkw) - ax3.set_ylabel(r'MSLP [hPa]') + + ax3.grid(axis="y") + ax3.plot(time_mod, mslp, color="red", **mplotkw) + ax3.set_ylabel(r"MSLP [hPa]") try: ax3.xaxis.set_ticklabels([], minor=True) ax3.xaxis.set_ticklabels([], minor=False) except: print("Plotting tick labels failed") - ax3.set_xlim(datetime_to_float(time[0] - datetime.timedelta(hours=3*len(time)/(24*3))), datetime_to_float(time[-1] + datetime.timedelta(hours=3*len(time)/(24*3)))) - ax3.set_ylim(min(mslp)-2, max(mslp)+2) - ax3.ticklabel_format(useOffset=False, axis='y') + ax3.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + ax3.set_ylim(min(mslp) - 2, max(mslp) + 2) + ax3.ticklabel_format(useOffset=False, axis="y") ax3.set_xticks([time_mod[i] for i in ind]) draw_vlines(ax3) # 4.--------------prec. plot------------------------ # rain: - - ax4.grid(axis='y') - ax4.bar(time_mod, ravg, color='c', label='tot. prec.', width=3600) - ax4.set_ylabel(r'total prec. [mm/h]') - ax4.set_ylim(0, max(max(ravg), 3)+1) - ax4.set_xlim(datetime_to_float(time[0] - datetime.timedelta(hours=3*len(time)/(24*3))), datetime_to_float(time[-1] + datetime.timedelta(hours=3*len(time)/(24*3)))) + + ax4.grid(axis="y") + ax4.bar(time_mod, ravg, color="c", label="tot. prec.", width=3600) + ax4.set_ylabel(r"total prec. [mm/h]") + ax4.set_ylim(0, max(max(ravg), 3) + 1) + ax4.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) ax4.set_xticks([time_mod[i] for i in ind], minor=True) - ax4.set_yticks([x for x in range(0, max(int(max(ravg)), 3)+2, max(int(max(ravg)), 3)//3)]) + ax4.set_yticks( + [x for x in range(0, max(int(max(ravg)), 3) + 2, max(int(max(ravg)), 3) // 3)] + ) try: - ax4.xaxis.set_ticklabels(['00', '06', '12', '18']*(len(time)//24), minor = True) + ax4.xaxis.set_ticklabels( + ["00", "06", "12", "18"] * (len(time) // 24), minor=True + ) except: print("Plotting tick labels failed") draw_vlines(ax4, labels=True) # print a red (magenta) thunderstorm symbol if cape > 500 (> 1500) Jkg-1 and cin < 100 Jkg-1 prop = FontProperties() - prop.set_file('/usr/share/fonts/truetype/freefont/FreeMono.ttf') + prop.set_file("/usr/share/fonts/truetype/freefont/FreeMono.ttf") for i, j in enumerate(time_mod): - if (cin[i]< 100 or np.isnan(cin[i])) and (cape[i]>=1500): - ax4.text(datetime_to_float(time[i]), ax4.get_ylim()[1]*0.62, u"\u2608", fontproperties=prop, color='magenta', ha='center', fontsize='large') - elif (cin[i]< 100 or np.isnan(cin[i])) and (cape[i]>=500): - ax4.text(datetime_to_float(time[i]), ax4.get_ylim()[1]*0.62, u"\u2608", fontproperties=prop, color='red', ha='center', fontsize='large') - + if (cin[i] < 100 or np.isnan(cin[i])) and (cape[i] >= 1500): + ax4.text( + datetime_to_float(time[i]), + ax4.get_ylim()[1] * 0.62, + "\u2608", + fontproperties=prop, + color="magenta", + ha="center", + fontsize="large", + ) + elif (cin[i] < 100 or np.isnan(cin[i])) and (cape[i] >= 500): + ax4.text( + datetime_to_float(time[i]), + ax4.get_ylim()[1] * 0.62, + "\u2608", + fontproperties=prop, + color="red", + ha="center", + fontsize="large", + ) # snow: ax5 = ax4.twinx() - ax5.bar(time_mod, snow, color='magenta', label='snow', width=3600) - ax5.set_ylabel(r'snow [cm/h]') + ax5.bar(time_mod, snow, color="magenta", label="snow", width=3600) + ax5.set_ylabel(r"snow [cm/h]") ax5.set_ylim(0, max(max(ravg), 3) + 1) - ax5.set_xlim(datetime_to_float(time[0] - datetime.timedelta(hours=3*len(time)/(24*3))), datetime_to_float(time[-1] + datetime.timedelta(hours=3*len(time)/(24*3)))) - ax5.set_yticks([x for x in range(max(int(max(ravg)), 3)+2) if x%(max(int(max(ravg)), 3)/3)==0]) + ax5.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + ax5.set_yticks( + [ + x + for x in range(max(int(max(ravg)), 3) + 2) + if x % (max(int(max(ravg)), 3) / 3) == 0 + ] + ) ax5.set_xticks([time_mod[i] for i in ind]) try: ax5.yaxis.set_ticklabels([]) @@ -458,21 +640,53 @@ def plot_var(filenames, outfile, deg_round_loc=None, print("Plotting tick labels failed") lines, labels = ax4.get_legend_handles_labels() lines2, labels2 = ax5.get_legend_handles_labels() - ax4.legend(lines + lines2, labels + labels2, loc='best', fontsize='small', ncol=2) + ax4.legend(lines + lines2, labels + labels2, loc="best", fontsize="small", ncol=2) # ---------------additional----------------- - - date = datetime.datetime.utcnow().strftime("%H:%M, %B %d, %Y, UTC") - ax0.text(1, 1.01, "{0:} \n experimental WRF setup \n $\Delta x = {1:.3f}$ km".format(date, compute_average_resolution(DS[0], *extent) / 1e3), fontsize='x-small', transform=ax0.transAxes, horizontalalignment='right') - plt.figtext(0.483, 0.948, r'$\copyright$ Meteorological Institute Munich, LMU, Fabian Jakub & Matjaz Puh', fontsize='x-small', alpha=.4) - ax0.text(0.05, 1.15, "MIM Forecast", transform=ax0.transAxes, fontsize='xx-large', horizontalalignment='left') - ax0.text(0.05, 1.05, " Meteogram @{0:.5f},{1:.5f} \n using {2:} pts for output statistics".format(location_lon, location_lat, np.size(xlon)), transform=ax0.transAxes, fontsize='small', horizontalalignment='left') - - #plt.ticklabel_format(useoffset=False, axis='y') - plt.savefig(outfile, bbox_inches='tight') - -if __name__ == '__main__': + date = datetime.datetime.utcnow().strftime("%H:%M, %B %d, %Y, UTC") + ax0.text( + 1, + 1.01, + "{0:} \n experimental WRF setup \n $\Delta x = {1:.3f}$ km".format( + date, compute_average_resolution(DS[0], *extent) / 1e3 + ), + fontsize="x-small", + transform=ax0.transAxes, + horizontalalignment="right", + ) + + plt.figtext( + 0.483, + 0.948, + r"Christoph Knote (MBEES, Uni Augsburg), Fabian Jakub, Matjaz Puh (Meteorological Institute, LMU Munich)", + fontsize="x-small", + alpha=0.4, + ) + ax0.text( + 0.05, + 1.15, + "MIM Forecast", + transform=ax0.transAxes, + fontsize="xx-large", + horizontalalignment="left", + ) + ax0.text( + 0.05, + 1.05, + " Meteogram @{0:.5f},{1:.5f} \n using {2:} pts for output statistics".format( + location_lon, location_lat, np.size(xlon) + ), + transform=ax0.transAxes, + fontsize="small", + horizontalalignment="left", + ) + + # plt.ticklabel_format(useoffset=False, axis='y') + plt.savefig(outfile, bbox_inches="tight") + + +if __name__ == "__main__": import datetime from argparse import ArgumentParser, ArgumentTypeError @@ -483,39 +697,61 @@ if __name__ == '__main__': raise ArgumentTypeError("Not a valid date: '{0}'.".format(s)) parser = ArgumentParser() - parser.add_argument('wrf_data_fpath_pattern', type=str, help='WRF input data path pattern (including strptime pattern for date and {domain:s} for domain (01, 02, ...))') - parser.add_argument('start_date', type=valid_date, help='First date to process (YYYY-mm-dd)') - parser.add_argument('end_date', type=valid_date, help='Last date to process (YYYY-mm-dd)') - parser.add_argument('interval_hours', type=int, help='Time interval in hours for plotting') - parser.add_argument('domain', type=int, help='Which domain to process') - parser.add_argument('plot_fpath_pattern', type=str, help='Plot output path pattern (including strptime pattern for date and {plot:s} for plot identifier (e.g., 500hPa_europe...)') - parser.add_argument('--spatial_patch_deg', type=float, default=None, help="Size of the patch (in degrees) to use for uncertainty calculations") + parser.add_argument( + "wrf_data_fpath_pattern", + type=str, + help="WRF input data path pattern (including strptime pattern for date and {domain:s} for domain (01, 02, ...))", + ) + parser.add_argument( + "start_date", type=valid_date, help="First date to process (YYYY-mm-dd)" + ) + parser.add_argument( + "end_date", type=valid_date, help="Last date to process (YYYY-mm-dd)" + ) + parser.add_argument( + "interval_hours", type=int, help="Time interval in hours for plotting" + ) + parser.add_argument("domain", type=int, help="Which domain to process") + parser.add_argument( + "plot_fpath_pattern", + type=str, + help="Plot output path pattern (including strptime pattern for date and {plot:s} for plot identifier (e.g., 500hPa_europe...)", + ) + parser.add_argument( + "--spatial_patch_deg", + type=float, + default=0.1, + help="Size of the patch (in degrees) to use for uncertainty calculations", + ) args = parser.parse_args() - #args = parser.parse_args( [ "/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout_d{domain:02d}_%Y-%m-%d_%H:%M:%S", "2021-03-11", "2021-03-14", "3", "1", "{plot:s}_%Y%m%d%H.png" ]) + # args = parser.parse_args( [ "/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout_d{domain:02d}_%Y-%m-%d_%H:%M:%S", "2021-03-11", "2021-03-14", "3", "1", "{plot:s}_%Y%m%d%H.png" ]) import numpy as np - - filenames = [ d.strftime(args.wrf_data_fpath_pattern).format(domain=args.domain) for d in np.arange(args.start_date, args.end_date, datetime.timedelta(hours=args.interval_hours), dtype=datetime.datetime) ] - + + filenames = [ + d.strftime(args.wrf_data_fpath_pattern).format(domain=args.domain) + for d in np.arange( + args.start_date, + args.end_date, + datetime.timedelta(hours=args.interval_hours), + dtype=datetime.datetime, + ) + ] + try: - plot_var(filenames, args.start_date.strftime(args.plot_fpath_pattern.format(plot="met")), deg_round_loc=args.spatial_patch_deg) + plot_var( + filenames, + args.start_date.strftime(args.plot_fpath_pattern.format(plot="met")), + deg_round_loc=args.spatial_patch_deg, + ) except: import matplotlib.font_manager as fm - + print("Rebuilding font cache") fm._rebuild() - plot_var(filenames, args.start_date.strftime(args.plot_fpath_pattern.format(plot="met")), deg_round_loc=args.spatial_patch_deg) - - - - - - - - - - - - - + plot_var( + filenames, + args.start_date.strftime(args.plot_fpath_pattern.format(plot="met")), + deg_round_loc=args.spatial_patch_deg, + ) diff --git a/blueprints/test_chemistry/pp/plotting/meteograms/run.sh b/blueprints/test_chemistry/pp/plotting/meteograms/run.sh index a8e5666667b53d635d4d66622a11627eeb4ca27a..bac917862bea4ad109826ea65f60682b3db868b6 100755 --- a/blueprints/test_chemistry/pp/plotting/meteograms/run.sh +++ b/blueprints/test_chemistry/pp/plotting/meteograms/run.sh @@ -17,7 +17,7 @@ wrfDataPathPattern="/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_c firstDateToProcess=`date --date="2 days ago" "+%Y-%m-%d"` lastDateToProcess=`date --date="3 days" "+%Y-%m-%d"` intervalHours=1 -domain=1 +domain=2 plotDir="/alcc/gpfs2/scratch/mbees/knotechr/plots/WRF/operational_chemistry/meteograms" plotPathPattern="${plotDir}/{plot:s}_%Y%m%d%H.png" diff --git a/tools/CAMS_split.py b/tools/CAMS_split.py new file mode 100644 index 0000000000000000000000000000000000000000..69536da27b4f2c3027cd1ccb83edbb38daf01cf9 --- /dev/null +++ b/tools/CAMS_split.py @@ -0,0 +1,33 @@ +import netCDF4 + +import yaml + +import numpy as np + +with open("CAMS_split.yaml", "r") as stream: + config = yaml.safe_load(stream) + +nc_blueprint_total = netCDF4.Dataset(config["blueprints_path_pattern"].format(varname=config["blueprint_sum_var"])) + +out = {} + +for species in config["blueprint_frac_vars"]: + nc_blueprint_var = netCDF4.Dataset(config["blueprints_path_pattern"].format(varname=species)) + out[species] = {} + for in_sector, out_sector in config["sector_translator"].items(): + if not out_sector in out[species]: + out[species][out_sector] = np.zeros_like(nc_blueprint_total[in_sector][:]) + if in_sector in nc_blueprint_var.variables: + out[species][out_sector] += nc_blueprint_var[in_sector][:] / nc_blueprint_total[in_sector][:] + print(f"{species} / {in_sector} / {out_sector}: { np.mean(out[species][out_sector]) }") + nc_blueprint_var.close() + +nc_blueprint_total.close() + + + + +total_species = "nmvoc" +total_species_varname = "SumAllSectors" + + diff --git a/tools/CAMS_split.yaml b/tools/CAMS_split.yaml new file mode 100644 index 0000000000000000000000000000000000000000..b595e17c22d3899d7efadde8ab1f4d4d40c3af9e --- /dev/null +++ b/tools/CAMS_split.yaml @@ -0,0 +1,35 @@ +blueprints_path_pattern: "EDGARv5/MOZART_MOSAIC/edgar_v5_2015_{varname:s}_0.1x0.1.nc" +data_path_pattern: "CAMS/REG-AP/CAMS-REG-AP_EUR_0.05x0.1_anthro_{varname:s}_v5.1_yearly_2018.nc" + +blueprint_sum_var: NMVOC +blueprint_frac_vars: [ BENZENE, BIGALK, BIGENE, C2H2, C2H4, C2H5OH, C2H6, C3H6, C3H8, CH2O, CH3CHO, CH3COCH3, CH3COOH, CH3OH, HCOOH, MEK, TOLUENE, XYLENES ] + +sector_translator: { + AWB: L_AgriOther, + AGS: L_AgriOther, + CHE: B_Industry, + ENE: A_PublicPower, + FFF: C_OtherStationaryComb, + FOO-PAP: B_Industry, + IND: B_Industry, + IRO: B_Industry, + MNM: K_AgriLivestock, + NMM: B_Industry, + NFE: B_Industry, + NEU: B_Industry, + PRO: D_Fugitives, + PRU-SOL: E_Solvents, + RCO: A_PublicPower, + REF-TRF: B_Industry, + SWD-INC: J_Waste, + SWD-LDF: J_Waste, + TNR-Aviation-CDS: H_Aviation, + TNR-Aviation-CRS: H_Aviation, + TNR-Aviation-LTO: H_Aviation, + TNR-Other: I_OffRoad, + TNR-Ship: G_Shipping, + TOTAL: SumAllSectors, + TRO-RES: F_RoadTransport, + TRO-noRES: F_RoadTransport, + WWT: J_Waste +} diff --git a/tools/cams_emissions_prep.bash b/tools/cams_emissions_prep.bash new file mode 100644 index 0000000000000000000000000000000000000000..594d7543d20b04f19964c2e916040166dda50a96 --- /dev/null +++ b/tools/cams_emissions_prep.bash @@ -0,0 +1,69 @@ +#!/bin/bash + +ls -1 EDGARv5/MOZART_MOSAIC/ | sed -e "s/edgar_v5_2015_//g" | sed -e "s/_0.1x0.1.nc//g" + +BC +BENZENE +BIGALK +BIGENE +C2H2 +C2H4 +C2H5OH +C2H6 +C3H6 +C3H8 +CH2O +CH3CHO +CH3COCH3 +CH3COOH +CH3OH +CO +HCOOH +MEK +NH3 +NMVOC +NOx +OC +PM10 +PM2.5 +SO2 +TOLUENE +XYLENES + +ls -1 CAMS/REG-AP/*2018* | sed -e "s|CAMS/REG-AP/CAMS-REG-AP_EUR_0.05x0.1_anthro_||g" | sed -e "s/_v5.1_yearly_2018.nc//g" + + - ch4 -> ignore + - co + - nh3 + - nmvoc +BENZENE +BIGALK +BIGENE +C2H2 +C2H4 +C2H5OH +C2H6 +C3H6 +C3H8 +CH2O +CH3CHO +CH3COCH3 +CH3COOH +CH3OH +HCOOH +MEK +TOLUENE +XYLENES + + + + - nox + - pm10 + - pm2_5 + - so2 + + +module load nco + +module load python py-numpy py-netcdf4 py-yaml +