From fab36238f18a1219ed897a8e6acaa893aae836c8 Mon Sep 17 00:00:00 2001 From: George McCabe <23407799+georgemccabe@users.noreply.github.com> Date: Tue, 15 Oct 2024 09:22:57 -0600 Subject: [PATCH] Feature #50 Hurricane Matthew: Add WRF plotting to METplus run (#66) * Initial WRF plotting script, functions, and README for the Hurricane Matthew case. * add Dockerfile to create iwrf-metplus image and update instructions to use that container * make script executable * create directory containing output file if it does not already exist * change default input and output directories to use paths assumed in Docker container * mount directory containing WRF plotting script to container * run WRF plotting script as part of METplus use case * Per #50, update instructions to use I-WRF METplus Docker container that includes METplotpy, METcalcpy, METdataio, and wrf python packages. Mount the Hurricane Matthew Visualization directory so the plotting scripts can be found. The WRF plotting script is called from the METplus use case run * remove commented code * remove characters that should not be rendered in the docs * cleanup usage statement * Per request on @Trumbore on pull request review, updated WRF log output to include updated content. --------- Co-authored-by: Jared A. Lee Co-authored-by: John Halley Gotway --- docker/metplus/Dockerfile | 9 + docs/Users_Guide/matthewjetstream.rst | 32 +- docs/Users_Guide/matthewredcloud.rst | 32 +- docs/Users_Guide/matthewwindows.rst | 32 +- docs/Users_Guide/running.rst | 24 +- .../METplus/PointStat_matthew.conf | 7 +- .../Hurricane_Matthew/Visualization/README.md | 38 + .../Visualization/map_funcs.py | 526 +++++++++++++ .../Visualization/plot_wrf.py | 703 ++++++++++++++++++ 9 files changed, 1359 insertions(+), 44 deletions(-) create mode 100644 docker/metplus/Dockerfile create mode 100644 use_cases/Hurricane_Matthew/Visualization/README.md create mode 100644 use_cases/Hurricane_Matthew/Visualization/map_funcs.py create mode 100755 use_cases/Hurricane_Matthew/Visualization/plot_wrf.py diff --git a/docker/metplus/Dockerfile b/docker/metplus/Dockerfile new file mode 100644 index 0000000..c649ace --- /dev/null +++ b/docker/metplus/Dockerfile @@ -0,0 +1,9 @@ +ARG METPLUS_TAG=6.0.0-beta5 + +FROM dtcenter/metplus-analysis:${METPLUS_TAG} + +# install wrf python package from branch in PR NCAR/wrf-python#242 +# because wrf will not install from pip +RUN python3 -m pip install --upgrade pip \ + && python3 -m pip install git+https://github.com/DWesl/wrf-python.git@cmake-build \ + && python3 -m pip install cartopy diff --git a/docs/Users_Guide/matthewjetstream.rst b/docs/Users_Guide/matthewjetstream.rst index 58e868b..082c6ff 100644 --- a/docs/Users_Guide/matthewjetstream.rst +++ b/docs/Users_Guide/matthewjetstream.rst @@ -181,12 +181,13 @@ make sure that we refer to the same resource names and file paths wherever they Copy and paste the definitions below into your shell to define the variables before proceeding:: WRF_IMAGE=ncar/iwrf:latest - METPLUS_IMAGE=dtcenter/metplus-dev:develop + METPLUS_IMAGE=ncar/iwrf:metplus-latest WORKING_DIR=/home/exouser WRF_DIR=${WORKING_DIR}/wrf/20161006_00 METPLUS_DIR=${WORKING_DIR}/metplus WRF_CONFIG_DIR=${WORKING_DIR}/i-wrf/use_cases/Hurricane_Matthew/WRF METPLUS_CONFIG_DIR=${WORKING_DIR}/i-wrf/use_cases/Hurricane_Matthew/METplus + PLOT_SCRIPT_DIR=${WORKING_DIR}/i-wrf/use_cases/Hurricane_Matthew/Visualization OBS_DATA_VOL=data-matthew-input-obs Any time you open a new shell on your instance, you will need to perform this action @@ -334,16 +335,16 @@ Once completed, you can view the end of an output file to confirm that it succee The output should look something like this:: - Timing for main: time 2016-10-06_11:42:30 on domain 1: 0.23300 elapsed seconds - Timing for main: time 2016-10-06_11:45:00 on domain 1: 0.23366 elapsed seconds - Timing for main: time 2016-10-06_11:47:30 on domain 1: 2.77688 elapsed seconds - Timing for main: time 2016-10-06_11:50:00 on domain 1: 0.23415 elapsed seconds - Timing for main: time 2016-10-06_11:52:30 on domain 1: 0.23260 elapsed seconds - Timing for main: time 2016-10-06_11:55:00 on domain 1: 0.23354 elapsed seconds - Timing for main: time 2016-10-06_11:57:30 on domain 1: 0.23345 elapsed seconds - Timing for main: time 2016-10-06_12:00:00 on domain 1: 0.23407 elapsed seconds - Timing for Writing wrfout_d01_2016-10-06_12:00:00 for domain 1: 0.32534 elapsed seconds - d01 2016-10-06_12:00:00 wrf: SUCCESS COMPLETE WRF + Timing for main: time 2016-10-07_23:50:00 on domain 1: 0.25548 elapsed seconds + Timing for main: time 2016-10-07_23:52:30 on domain 1: 0.25495 elapsed seconds + Timing for main: time 2016-10-07_23:55:00 on domain 1: 0.25066 elapsed seconds + Timing for main: time 2016-10-07_23:57:30 on domain 1: 0.25231 elapsed seconds + Timing for main: time 2016-10-08_00:00:00 on domain 1: 0.25795 elapsed seconds + Timing for Writing wrfout_d01_2016-10-08_00:00:00 for domain 1: 0.68666 elapsed seconds + Timing for Writing wrfout_zlev_d01_2016-10-08_00:00:00 for domain 1: 0.47411 elapsed seconds + Timing for Writing wrfout_plev_d01_2016-10-08_00:00:00 for domain 1: 0.47619 elapsed seconds + Timing for Writing restart for domain 1: 1.54598 elapsed seconds + d01 2016-10-08_00:00:00 wrf: SUCCESS COMPLETE WRF Run METplus =========== @@ -351,12 +352,17 @@ Run METplus After the WRF simulation has finished, you can run the METplus verification to compare the simulated results to the actual weather observations during the hurricane. The verification takes about five minutes to complete. -We use command line options to tell the METplus container several things, including where the observed data is located, -where the METplus configuration can be found, where the WRF output data is located, and where it should create its output files:: +We use command line options to tell the METplus container several things, +including where the observed data is located, +where the METplus configuration can be found, +where the plotting scripts can be found, +where the WRF output data is located, +and where it should create its output files:: docker run --rm -it \ --volumes-from ${OBS_DATA_VOL} \ -v ${METPLUS_CONFIG_DIR}:/config \ + -v ${PLOT_SCRIPT_DIR}:/plot_scripts \ -v ${WORKING_DIR}/wrf:/data/input/wrf \ -v ${METPLUS_DIR}:/data/output ${METPLUS_IMAGE} \ /metplus/METplus/ush/run_metplus.py /config/PointStat_matthew.conf diff --git a/docs/Users_Guide/matthewredcloud.rst b/docs/Users_Guide/matthewredcloud.rst index 78b850f..46c4ed1 100644 --- a/docs/Users_Guide/matthewredcloud.rst +++ b/docs/Users_Guide/matthewredcloud.rst @@ -187,12 +187,13 @@ make sure that we refer to the same resource names and file paths wherever they Copy and paste the definitions below into your shell to define the variables before proceeding:: WRF_IMAGE=ncar/iwrf:latest - METPLUS_IMAGE=dtcenter/metplus-dev:develop + METPLUS_IMAGE=ncar/iwrf:metplus-latest WORKING_DIR=/home/ubuntu WRF_DIR=${WORKING_DIR}/wrf/20161006_00 METPLUS_DIR=${WORKING_DIR}/metplus WRF_CONFIG_DIR=${WORKING_DIR}/i-wrf/use_cases/Hurricane_Matthew/WRF METPLUS_CONFIG_DIR=${WORKING_DIR}/i-wrf/use_cases/Hurricane_Matthew/METplus + PLOT_SCRIPT_DIR=${WORKING_DIR}/i-wrf/use_cases/Hurricane_Matthew/Visualization OBS_DATA_VOL=data-matthew-input-obs Any time you open a new shell on your instance, you will need to perform this action @@ -340,16 +341,16 @@ Once completed, you can view the end of an output file to confirm that it succee The output should look something like this:: - Timing for main: time 2016-10-06_11:42:30 on domain 1: 0.23300 elapsed seconds - Timing for main: time 2016-10-06_11:45:00 on domain 1: 0.23366 elapsed seconds - Timing for main: time 2016-10-06_11:47:30 on domain 1: 2.77688 elapsed seconds - Timing for main: time 2016-10-06_11:50:00 on domain 1: 0.23415 elapsed seconds - Timing for main: time 2016-10-06_11:52:30 on domain 1: 0.23260 elapsed seconds - Timing for main: time 2016-10-06_11:55:00 on domain 1: 0.23354 elapsed seconds - Timing for main: time 2016-10-06_11:57:30 on domain 1: 0.23345 elapsed seconds - Timing for main: time 2016-10-06_12:00:00 on domain 1: 0.23407 elapsed seconds - Timing for Writing wrfout_d01_2016-10-06_12:00:00 for domain 1: 0.32534 elapsed seconds - d01 2016-10-06_12:00:00 wrf: SUCCESS COMPLETE WRF + Timing for main: time 2016-10-07_23:50:00 on domain 1: 0.25548 elapsed seconds + Timing for main: time 2016-10-07_23:52:30 on domain 1: 0.25495 elapsed seconds + Timing for main: time 2016-10-07_23:55:00 on domain 1: 0.25066 elapsed seconds + Timing for main: time 2016-10-07_23:57:30 on domain 1: 0.25231 elapsed seconds + Timing for main: time 2016-10-08_00:00:00 on domain 1: 0.25795 elapsed seconds + Timing for Writing wrfout_d01_2016-10-08_00:00:00 for domain 1: 0.68666 elapsed seconds + Timing for Writing wrfout_zlev_d01_2016-10-08_00:00:00 for domain 1: 0.47411 elapsed seconds + Timing for Writing wrfout_plev_d01_2016-10-08_00:00:00 for domain 1: 0.47619 elapsed seconds + Timing for Writing restart for domain 1: 1.54598 elapsed seconds + d01 2016-10-08_00:00:00 wrf: SUCCESS COMPLETE WRF Run METplus =========== @@ -357,12 +358,17 @@ Run METplus After the WRF simulation has finished, you can run the METplus verification to compare the simulated results to the actual weather observations during the hurricane. The verification takes about five minutes to complete. -We use command line options to tell the METplus container several things, including where the observed data is located, -where the METplus configuration can be found, where the WRF output data is located, and where it should create its output files:: +We use command line options to tell the METplus container several things, +including where the observed data is located, +where the METplus configuration can be found, +where the plotting scripts can be found, +where the WRF output data is located, +and where it should create its output files:: sudo docker run --rm -it \ --volumes-from ${OBS_DATA_VOL} \ -v ${METPLUS_CONFIG_DIR}:/config \ + -v ${PLOT_SCRIPT_DIR}:/plot_scripts \ -v ${WORKING_DIR}/wrf:/data/input/wrf \ -v ${METPLUS_DIR}:/data/output ${METPLUS_IMAGE} \ /metplus/METplus/ush/run_metplus.py /config/PointStat_matthew.conf diff --git a/docs/Users_Guide/matthewwindows.rst b/docs/Users_Guide/matthewwindows.rst index c0a7266..87605b7 100644 --- a/docs/Users_Guide/matthewwindows.rst +++ b/docs/Users_Guide/matthewwindows.rst @@ -62,11 +62,12 @@ or changes the path entirely to use a different location on your computer:: Now you can copy and paste the definitions below into your shell to define the other variables before proceeding:: set WRF_IMAGE=ncar/iwrf:latest - set METPLUS_IMAGE=dtcenter/metplus-dev:develop + set METPLUS_IMAGE=ncar/iwrf:metplus-latest set WRF_DIR=%WORKING_DIR%\wrf\20161006_00 set METPLUS_DIR=%WORKING_DIR%\metplus set WRF_CONFIG_DIR=%WORKING_DIR%\i-wrf-main\use_cases\Hurricane_Matthew\WRF set METPLUS_CONFIG_DIR=%WORKING_DIR%\i-wrf-main\use_cases\Hurricane_Matthew\METplus + set PLOT_SCRIPT_DIR=%WORKING_DIR%\i-wrf-main\use_cases\Hurricane_Matthew\Visualization set OBS_DATA_VOL=data-matthew-input-obs Any time you open a new shell on your instance, you will need to perform this action @@ -223,16 +224,16 @@ Once completed, you can view the end of an output file to confirm that it succee The output should look something like this:: - Timing for main: time 2016-10-06_11:42:30 on domain 1: 0.23300 elapsed seconds - Timing for main: time 2016-10-06_11:45:00 on domain 1: 0.23366 elapsed seconds - Timing for main: time 2016-10-06_11:47:30 on domain 1: 2.77688 elapsed seconds - Timing for main: time 2016-10-06_11:50:00 on domain 1: 0.23415 elapsed seconds - Timing for main: time 2016-10-06_11:52:30 on domain 1: 0.23260 elapsed seconds - Timing for main: time 2016-10-06_11:55:00 on domain 1: 0.23354 elapsed seconds - Timing for main: time 2016-10-06_11:57:30 on domain 1: 0.23345 elapsed seconds - Timing for main: time 2016-10-06_12:00:00 on domain 1: 0.23407 elapsed seconds - Timing for Writing wrfout_d01_2016-10-06_12:00:00 for domain 1: 0.32534 elapsed seconds - d01 2016-10-06_12:00:00 wrf: SUCCESS COMPLETE WRF + Timing for main: time 2016-10-07_23:50:00 on domain 1: 0.25548 elapsed seconds + Timing for main: time 2016-10-07_23:52:30 on domain 1: 0.25495 elapsed seconds + Timing for main: time 2016-10-07_23:55:00 on domain 1: 0.25066 elapsed seconds + Timing for main: time 2016-10-07_23:57:30 on domain 1: 0.25231 elapsed seconds + Timing for main: time 2016-10-08_00:00:00 on domain 1: 0.25795 elapsed seconds + Timing for Writing wrfout_d01_2016-10-08_00:00:00 for domain 1: 0.68666 elapsed seconds + Timing for Writing wrfout_zlev_d01_2016-10-08_00:00:00 for domain 1: 0.47411 elapsed seconds + Timing for Writing wrfout_plev_d01_2016-10-08_00:00:00 for domain 1: 0.47619 elapsed seconds + Timing for Writing restart for domain 1: 1.54598 elapsed seconds + d01 2016-10-08_00:00:00 wrf: SUCCESS COMPLETE WRF Run METplus =========== @@ -240,12 +241,17 @@ Run METplus After the WRF simulation has finished, you can run the METplus verification to compare the simulated results to the actual weather observations during the hurricane. The verification takes about five minutes to complete. -We use command line options to tell the METplus container several things, including where the observed data is located, -where the METplus configuration can be found, where the WRF output data is located, and where it should create its output files:: +We use command line options to tell the METplus container several things, +including where the observed data is located, +where the METplus configuration can be found, +where the plotting scripts can be found, +where the WRF output data is located, +and where it should create its output files:: docker run --rm -it ^ --volumes-from %OBS_DATA_VOL% ^ -v %METPLUS_CONFIG_DIR%:/config ^ + -v %PLOT_SCRIPT_DIR%:/plot_scripts ^ -v %WORKING_DIR%\wrf:/data/input/wrf ^ -v %METPLUS_DIR%:/data/output %METPLUS_IMAGE% ^ /metplus/METplus/ush/run_metplus.py /config/PointStat_matthew.conf diff --git a/docs/Users_Guide/running.rst b/docs/Users_Guide/running.rst index d0f6f9e..6dda4dd 100644 --- a/docs/Users_Guide/running.rst +++ b/docs/Users_Guide/running.rst @@ -17,12 +17,24 @@ Load the apptainer module:: module load apptainer -Change directory to scratch and pull the containers from DockerHub. +Change directory to scratch and pull the images from DockerHub. This will create a `.sif` file in the current directory:: - apptainer pull ${SCRATCH}/metplus-dev_develop.sif docker://dtcenter/metplus-dev:develop + apptainer pull ${SCRATCH}/iwrf-metplus.sif docker://ncar/iwrf:metplus-latest apptainer pull ${SCRATCH}/data-matthew-input-obs.sif oras://registry-1.docker.io/ncar/iwrf:data-matthew-input-obs +.. note:: + + If an error is displayed when attempting to pull the METplus image, + creating a DockerHub account and authenticating through apptainer may be + necessary. + + :: + + apptainer remote login --username {USERNAME} docker://docker.io + + where **{USERNAME}** is your DockerHub username. + Create a directory to store the output:: mkdir ${SCRATCH}/metplus_out @@ -48,6 +60,9 @@ using the --bind argument) * Config directory containing METplus use case configuration file * Local: ${SCRATCH}/i-wrf/use_cases/Hurricane_Matthew/METplus * Container: /config +* Plot script directory containing WRF plotting scripts + * Local: ${SCRATCH}/i-wrf/use_cases/Hurricane_Matthew/Visualization + * Container: /plot_scripts * Output directory to write output * Local: ${SCRATCH}/metplus_out * Container: /data/output @@ -55,14 +70,15 @@ using the --bind argument) :: LOCAL_METPLUS_CONFIG_DIR=${SCRATCH}/i-wrf/use_cases/Hurricane_Matthew/METplus + LOCAL_PLOT_SCRIPT_DIR=${SCRATCH}/i-wrf/use_cases/Hurricane_Matthew/Visualization LOCAL_FCST_INPUT_DIR=/glade/derecho/scratch/jaredlee/nsf_i-wrf/matthew LOCAL_OUTPUT_DIR=${SCRATCH}/metplus_out - export APPTAINER_BIND="${SCRATCH}/data-matthew-input-obs.sif:/data/input/obs:image-src=/,${LOCAL_METPLUS_CONFIG_DIR}:/config,${LOCAL_FCST_INPUT_DIR}:/data/input/wrf,${LOCAL_OUTPUT_DIR}:/data/output" + export APPTAINER_BIND="${SCRATCH}/data-matthew-input-obs.sif:/data/input/obs:image-src=/,${LOCAL_METPLUS_CONFIG_DIR}:/config,${LOCAL_FCST_INPUT_DIR}:/data/input/wrf,${LOCAL_OUTPUT_DIR}:/data/output,${LOCAL_PLOT_SCRIPT_DIR}:/plot_scripts" Execute the run_metplus.py command inside the container to run the use case:: - apptainer exec ${SCRATCH}/metplus-dev_develop.sif /metplus/METplus/ush/run_metplus.py /config/PointStat_matthew.conf + apptainer exec ${SCRATCH}/iwrf-metplus.sif /metplus/METplus/ush/run_metplus.py /config/PointStat_matthew.conf Check that the output data was created locally:: diff --git a/use_cases/Hurricane_Matthew/METplus/PointStat_matthew.conf b/use_cases/Hurricane_Matthew/METplus/PointStat_matthew.conf index fb96f9a..931b7aa 100644 --- a/use_cases/Hurricane_Matthew/METplus/PointStat_matthew.conf +++ b/use_cases/Hurricane_Matthew/METplus/PointStat_matthew.conf @@ -8,7 +8,7 @@ # https://metplus.readthedocs.io/en/latest/Users_Guide/systemconfiguration.html#process-list ### -PROCESS_LIST = MADIS2NC(metar), MADIS2NC(raob), PointStat(surface), PointStat(upper_air) +PROCESS_LIST = MADIS2NC(metar), MADIS2NC(raob), PointStat(surface), PointStat(upper_air), UserScript(wrf_plot) ### @@ -163,3 +163,8 @@ OBTYPE = POINT_STAT_OUTPUT_PREFIX = {instance} POINT_STAT_MASK_GRID = FULL + +[wrf_plot] + +USER_SCRIPT_RUNTIME_FREQ = RUN_ONCE +USER_SCRIPT_COMMAND=/plot_scripts/plot_wrf.py diff --git a/use_cases/Hurricane_Matthew/Visualization/README.md b/use_cases/Hurricane_Matthew/Visualization/README.md new file mode 100644 index 0000000..a25dfcf --- /dev/null +++ b/use_cases/Hurricane_Matthew/Visualization/README.md @@ -0,0 +1,38 @@ +# plot_wrf +Python scripts and functions to make nice plots of WRF model output using Cartopy. + +plot_wrf.py can be called with a set of arguments (use option -h to get the current usage statement): + +``` +> python plot_wrf.py -h +usage: plot_wrf.py [-h] [-w WRF_DIR_PARENT] [-o OUT_DIR_PARENT] [-f CYCLE_DT_FIRST] [-l CYCLE_DT_LAST] + [-i CYCLE_STRIDE_H] [-b BEG_LEAD_TIME] [-e END_LEAD_TIME] [-s STR_LEAD_TIME] [-d DOMAIN] + +options: + -h, --help show this help message and exit + -w WRF_DIR_PARENT, --wrf_dir_parent WRF_DIR_PARENT + string specifying the directory path to the parent WRF output directories, above any + experiment or cycle datetime subdirectories (default: /data/input/wrf) + -o OUT_DIR_PARENT, --out_dir_parent OUT_DIR_PARENT + string specifying the directory path to the parent plot directories (default: /data/output/wrf) + -f CYCLE_DT_FIRST, --cycle_dt_first CYCLE_DT_FIRST + beginning date/time of first WRF simulation [YYYYMMDD_HH] (default: 20161006_00) + -l CYCLE_DT_LAST, --cycle_dt_last CYCLE_DT_LAST + beginning date/time of last WRF simulation [YYYYMMDD_HH] + -i CYCLE_STRIDE_H, --cycle_stride_h CYCLE_STRIDE_H + stride in hours between cycles (default: 24) + -b BEG_LEAD_TIME, --beg_lead_time BEG_LEAD_TIME + beginning lead time for plotting WRF simulations [HH:MM] (default: 00:00) + -e END_LEAD_TIME, --end_lead_time END_LEAD_TIME + ending lead time for plotting WRF simulations [HH:MM] (default: 48:00) + -s STR_LEAD_TIME, --str_lead_time STR_LEAD_TIME + stride to create plots every N minutes (default: 180) + -d DOMAIN, --domain DOMAIN + WRF domain number to be plotted (default: 1) +``` + +The plot_wrf.parse_args function creates a dictionary of options that is then passed to the main routine. Doing this via a dictionary object should make it simpler to add even more customization/options in the future, requiring changes in fewer places than passing numerous positional arguments around. + +The plot_wrf script opens specified wrfout files in sequence, reads in user-specified variables (currently set with options like plot_TERRAIN = True and plot_T2 = False in the main function), creates a dictionary of plotting options that is then passed to map_funcs.map_plot to create and save each plot to a file. Inside the main function there are also user-settable boolean flags to turn on/off plotting surface wind barb overlays, labeled stations/cities, etc. + +Both these requested variables for plotting and other plot customization options could eventually be changed to be passed in on the command line to not require users to modify the script itself before running it, but that is left for future development. diff --git a/use_cases/Hurricane_Matthew/Visualization/map_funcs.py b/use_cases/Hurricane_Matthew/Visualization/map_funcs.py new file mode 100644 index 0000000..70e9e76 --- /dev/null +++ b/use_cases/Hurricane_Matthew/Visualization/map_funcs.py @@ -0,0 +1,526 @@ +""" +map_funcs.py + +Created by: Jared A. Lee (jaredlee@ucar.edu) +Created on: 28 May 2024 + +This file contains common functions and procedures useful for plotting maps with Cartopy, +typically from WRF model output. +""" + +import sys +import os +import numpy as np +import datetime as dt +import wrf +import matplotlib as mpl +import matplotlib.pyplot as plt +import matplotlib.ticker as mticker +import cartopy +import cartopy.crs as ccrs +from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter +from cartopy.mpl.geoaxes import GeoAxes + +def calc_bearing(lon1, lat1, lon2, lat2): + """ + Function to calculate the bearing angle between two points (lon1, lat1) and (lon2, lat2). + -- Inputs: + - lon1, lat1: floats, longitude and latitude of point 1 + - lon2, lat2: floats, longitude and latitude of point 2 + -- Outputs: + - bearing_geog: float, geographical bearing (0 deg = north, 90 deg = east, etc.) + - bearing_math: float, mathematical bearing (0 deg = east, 90 deg = north, etc.) + -- Reference: + - https://www.igismap.com/formula-to-find-bearing-or-heading-angle-between-two-points-latitude-longitude/ + """ + DEG2RAD = np.pi / 180.0 + RAD2DEG = 180.0 / np.pi + x = np.cos(lat2 * DEG2RAD) * np.sin((lon2 - lon1) * DEG2RAD) + y = ((np.cos(lat1 * DEG2RAD) * np.sin(lat2 * DEG2RAD)) - + (np.sin(lat1 * DEG2RAD) * np.cos(lat2 * DEG2RAD) * np.cos((lon2-lon1)*DEG2RAD))) + bearing_geog = np.arctan2(x, y) * RAD2DEG + bearing_math = 90.0 - bearing_geog + # Perhaps unnecessary, but keeps it in the range (-180, 180] + if bearing_math <= -180.0: + bearing_math += 360.0 + + return bearing_geog, bearing_math + +def get_cartopy_features(): + """ + Function to get some commonly used Cartopy features (borders, states, oceans, lakes, rivers, land) + for use in Matplotlib/Cartopy plots. + -- Inputs: None. + -- Outputs: + - borders, states, oceans, lakes, rivers, land + """ + + borders = cartopy.feature.BORDERS + states = cartopy.feature.NaturalEarthFeature(category='cultural', scale='10m', facecolor='none', + name='admin_1_states_provinces_lakes') + oceans = cartopy.feature.OCEAN + lakes = cartopy.feature.LAKES + rivers = cartopy.feature.RIVERS + land = cartopy.feature.LAND + + return borders, states, oceans, lakes, rivers, land + +def truncate_cmap(cmap, minval=0.0, maxval=1.0, n=100): + """ + Function to truncate a matplotlib colormap. Particularly useful when plotting terrain. + -- Required Positional Inputs: + - cmap: Original Matplotlib colormap + -- Optional Inputs: + - minval: Fraction into the original colormap to start the new colormap (default: 0.0) + - maxval: Fraction into the original colormap to end the new colormap (default: 1.0) + - n: Number of colors to create in the new colormap (default: 100) + -- Output: + - new_cmap: New, truncated colormap + """ + new_cmap = mpl.colors.LinearSegmentedColormap.from_list( + 'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval), cmap(np.linspace(minval, maxval, n))) + return new_cmap + +def map_plot(opts): + """ + Procedure to make a map plot with filled contours using matplotlib and Cartopy. + Optionally overlays the map with wind barbs, markers, text labels, a polygon, or cross-section path. + -- Input: + - opts: Dictionary containing plotting options + Required keys: + - fname: string or pathlib object specifying the output file name + - fill_var: 2D variable array to be plotted with filled contours + - suptitle: string for overall plot title (usually one line) + - cbar_lab: string for colorbar label (e.g., 'Wind Speed [m/s]') + - cart_proj: Cartopy object, map projection + - lons: 1D or 2D array of longitude values + - lats: 1D or 2D array of latitude values + - cmap: Matplotlib colormap + - bounds: Matplotlib colormap bounds + - norm: matplotlib colormap norm + Optional keys: + - cart_xlim: Cartopy object, x-axis limits + - cart_ylim: Cartopy object, y-axis limits + - extend: string for colorbar caps ('max', 'min', 'both' [default]) + - fontsize: integer, base fontsize (default: 14) + - figsize: 2D tuple, defining the figure size (default: (10, 8)) + - cbar_loc: string, identifier for positioning of the colorbar ('bottom' [default], 'top', 'right', 'left') + - borders, states, oceans, lakes, rivers, land: Cartopy feature objects for the map + - border_width: numerical line thickness for national borders & coastlines (default: 1.5) + - water_color: string defining water color for the map (default: 'none' [transparent]) + - lat_labels: array of latitude values to label explicitly on the map + - lon_labels: array of longitude values to label explicitly on the map + - suptitle_y: float, y-axis position of the suptitle (default: 0.95) + - title_l: string, plot subtitle (1 or 2 lines) that gets placed above the top-left corner of the plot axes + - title_r: string, plot subtitle (1 or 2 lines) that gets placed above the top-right corner of the plot axes + - title_c: string, plot subtitle (1 or 2 lines) that gets placed above the center of the plot axes + - map_x_thin: integer, thin wind barb location overlays in x-direction (every Nth grid point) (default: 25) + - map_y_thin: integer, thin wind barb location overlays in y-direction (every Nth grid point) (default: 25) + - barb_width: float, linewidth of wind barbs (default: 0.25) + - u: array-like, define the barb directions + - v: array-like, define the barb directions + - mark1_lon: array of longitude values for set 1 of markers + - mark1_lat: array of latitude values for set 1 of markers + - mark1_size: integer specifying marker size for set 1 of markers (default: 100) + - mark1_style: string specifying marker style for set 1 of markers (default: 'o') + - mark1_color: string specifying the marker color for set 1 of markers + - mark1_edgecolor: string specifying the marker edge color for set 1 of markers (default: 'black') + - mark1_width: float specifying linewidth for set 1 of markers + - mark1_val_fill: boolean flag to fill set 1 of markers from data value, cmap, and norm (default: False) + - mark1_var: variable containing data to fill in set 1 of markers (e.g., plot obs station data) + - mark1_zorder: integer indicator of the Matplotlib draw order for set 1 of markers (default: 10) + - mark2_lon: array of longitude values for set 2 of markers + - mark2_lat: array of latitude values for set 2 of markers + - mark2_size: integer specifying marker size for set 2 of markers (default: 100) + - mark2_style: string specifying marker style for set 2 of markers (default: 'o') + - mark2_color: string specifying the marker color for set 2 of markers + - mark2_edgecolor: string specifying the marker edge color for set 2 of markers (default: 'black') + - mark2_width: float specifying linewidth for set 2 of markers + - mark2_val_fill: boolean flag to fill set 2 of markers from data value, cmap, and norm (default: False) + - mark2_var: variable containing data to fill in set 2 of markers (e.g., plot obs station data) + - mark2_zorder: integer indicator of the Matplotlib draw order for set 2 of markers (default: 11) + - text1_lat: array of floats for latitudes for set 1 of text + - text1_lon: array of floats for longitudes for set 1 of text + - text1_lab: array of strings for labels for set 1 of text + - text1_lab_wt: numeric or string indicating font weight for set 1 of text (default: 'normal') + - text2_lat: array of floats for latitudes for set 2 of text + - text2_lon: array of floats for longitudes for set 2 of text + - text2_lab: array of strings for labels for set 2 of text + - text2_lab_wt: numeric or string indicating font weight for set 2 of text (default: 'normal') + - lg_text: array of legend labels + - lg_loc: string, defining legend placement + - lg_fontsize: integer fontsize for legend labels + -- Output: + - generates a plot saved to fname + """ + # Set default values for optional dictionary entries + opts.setdefault('cart_xlim', None) + opts.setdefault('cart_ylim', None) + opts.setdefault('extend', 'both') + opts.setdefault('fontsize', 14) + opts.setdefault('figsize', (10,8)) + opts.setdefault('cbar_loc', 'bottom') + opts.setdefault('borders', None) + opts.setdefault('states', None) + opts.setdefault('oceans', None) + opts.setdefault('lakes', None) + opts.setdefault('rivers', None) + opts.setdefault('land', None) + opts.setdefault('water_color', 'none') + opts.setdefault('border_width', 1.5) + opts.setdefault('lat_labels', None) + opts.setdefault('lon_labels', None) + opts.setdefault('suptitle_y', 0.95) + opts.setdefault('title_l', None) + opts.setdefault('title_r', None) + opts.setdefault('title_c', None) + opts.setdefault('map_x_thin', 25) + opts.setdefault('map_y_thin', 25) + opts.setdefault('barb_width', 0.25) + opts.setdefault('u', None) + opts.setdefault('v', None) + opts.setdefault('mark1_lat', None) + opts.setdefault('mark1_lon', None) + opts.setdefault('mark1_size', 100) + opts.setdefault('mark1_style', 'o') + opts.setdefault('mark1_color', None) + opts.setdefault('mark1_edgecolor', 'black') + opts.setdefault('mark1_width', 1.5) + opts.setdefault('mark1_val_fill', False) + opts.setdefault('mark1_var', None) + opts.setdefault('mark1_zorder', 10) + opts.setdefault('mark2_lat', None) + opts.setdefault('mark2_lon', None) + opts.setdefault('mark2_size', 100) + opts.setdefault('mark2_style', 'o') + opts.setdefault('mark2_color', None) + opts.setdefault('mark2_edgecolor', 'black') + opts.setdefault('mark2_width', 1.5) + opts.setdefault('mark2_val_fill', False) + opts.setdefault('mark2_var', None) + opts.setdefault('mark2_zorder', 11) + opts.setdefault('text1_lat', None) + opts.setdefault('text1_lon', None) + opts.setdefault('text1_lab', None) + opts.setdefault('text1_lab_wt', 'normal') + opts.setdefault('text2_lat', None) + opts.setdefault('text2_lon', None) + opts.setdefault('text2_lab', None) + opts.setdefault('text2_lab_wt', 'normal') + opts.setdefault('lg_text', None) + opts.setdefault('lg_loc', 'lower left') + opts.setdefault('lg_fontsize', 14) + + # Pull everything out of the opts dict into variables for cleaner code later on + fname = opts['fname'] + fill_var = opts['fill_var'] + suptitle = opts['suptitle'] + cbar_lab = opts['cbar_lab'] + cart_proj = opts['cart_proj'] + cart_xlim = opts['cart_xlim'] + cart_ylim = opts['cart_ylim'] + lons = opts['lons'] + lats = opts['lats'] + cmap = opts['cmap'] + bounds = opts['bounds'] + norm = opts['norm'] + extend = opts['extend'] + fontsize = opts['fontsize'] + figsize = opts['figsize'] + cbar_loc = opts['cbar_loc'] + borders = opts['borders'] + states = opts['states'] + oceans = opts['oceans'] + lakes = opts['lakes'] + rivers = opts['rivers'] + land = opts['land'] + water_color = opts['water_color'] + border_width = opts['border_width'] + lat_labels = opts['lat_labels'] + lon_labels = opts['lon_labels'] + suptitle_y = opts['suptitle_y'] + title_l = opts['title_l'] + title_r = opts['title_r'] + title_c = opts['title_c'] + map_x_thin = opts['map_x_thin'] + map_y_thin = opts['map_y_thin'] + barb_width = opts['barb_width'] + u = opts['u'] + v = opts['v'] + mark1_lat = opts['mark1_lat'] + mark1_lon = opts['mark1_lon'] + mark1_size = opts['mark1_size'] + mark1_style = opts['mark1_style'] + mark1_color = opts['mark1_color'] + mark1_edgecolor = opts['mark1_edgecolor'] + mark1_width = opts['mark1_width'] + mark1_val_fill = opts['mark1_val_fill'] + mark1_var = opts['mark1_var'] + mark1_zorder = opts['mark1_zorder'] + mark2_lat = opts['mark2_lat'] + mark2_lon = opts['mark2_lon'] + mark2_size = opts['mark2_size'] + mark2_style = opts['mark2_style'] + mark2_color = opts['mark2_color'] + mark2_edgecolor = opts['mark2_edgecolor'] + mark2_width = opts['mark2_width'] + mark2_val_fill = opts['mark2_val_fill'] + mark2_var = opts['mark2_var'] + mark2_zorder = opts['mark2_zorder'] + text1_lat = opts['text1_lat'] + text1_lon = opts['text1_lon'] + text1_lab = opts['text1_lab'] + text1_lab_wt = opts['text1_lab_wt'] + text2_lat = opts['text2_lat'] + text2_lon = opts['text2_lon'] + text2_lab = opts['text2_lab'] + text2_lab_wt = opts['text2_lab_wt'] + lg_text = opts['lg_text'] + lg_loc = opts['lg_loc'] + lg_fontsize = opts['lg_fontsize'] + + # Set some Matplotlib resources + mpl.rcParams['figure.figsize'] = figsize + mpl.rcParams['grid.color'] = 'gray' + mpl.rcParams['grid.linestyle'] = ':' + mpl.rcParams['font.size'] = fontsize + 2 + mpl.rcParams['figure.titlesize'] = fontsize + 2 + mpl.rcParams['savefig.bbox'] = 'tight' + ll_size = fontsize - 2 + data_crs = ccrs.PlateCarree() + + print('-- Plotting ' + str(fname)) + + # Define the figure and axes + fig = plt.figure() + ax = plt.subplot(projection=cart_proj) + + # If cart_xlim and cart_ylim tuples are not provided, then set plot limits from lat/lon data directly + if cart_xlim is not None and cart_ylim is not None: + ax.set_xlim(cart_xlim) + ax.set_ylim(cart_ylim) + else: + ax.set_extent([np.min(lons), np.max(lons), np.min(lats), np.max(lats)], crs=cart_proj) + + # If lons & lats are 1D, make them into 2D arrays + if lons.ndim == 1 and lats.ndim == 1: + ll2d = np.meshgrid(lons, lats) + lons = ll2d[0] + lats = ll2d[1] + + # Optional: Add various cartopy features + if borders != None: + ax.add_feature(borders, linewidth=border_width, linestyle='-', zorder=3) + if states != None: + ax.add_feature(states, linewidth=border_width/3.0, edgecolor='black', zorder=4) + if oceans != None: + # Drawing the oceans can be VERY slow with Cartopy 0.20+ for some domains, so may want to skip it + # Can set the facecolor for the axes to water_color instead (usually we want this 'none' except for terrain) + ax.add_feature(oceans, facecolor=water_color, zorder=2) + # ax.add_feature(oceans, facecolor='none', zorder=2) + # ax.set_facecolor(opts['water_color']) + if lakes != None: + # Unless facecolor='none', lakes w/ facecolor will appear above filled contour plot, which is undesirable + ax.add_feature(lakes, facecolor=water_color, linewidth=0.25, edgecolor='black', zorder=5) + ax.coastlines(zorder=6, linewidth=border_width) + + # Sometimes longitude labels show up on y-axis, and latitude labels on x-axis in older versions of Cartopy + # Print lat/lon labels only for a specified set (determined by trial & error) to avoid this problem for now + gl = ax.gridlines(draw_labels=True, x_inline=False, y_inline=False) + gl.rotate_labels = False + # If specific lat/lon labels are not specified, then just label the default gridlines + if lon_labels is not None: + gl.xlocator = mticker.FixedLocator(lon_labels) + if lat_labels is not None: + gl.ylocator = mticker.FixedLocator(lat_labels) + gl.top_labels = True + gl.bottom_labels = True + gl.left_labels = True + gl.right_labels = True + gl.xlabel_style = {'size': ll_size} + gl.ylabel_style = {'size': ll_size} + + # Draw the actual filled contour plot + # NOTE: Sometimes a cartopy contourf plot may fail with a Shapely TopologicalError. + # It appears to happen sometimes when projecting a dataset onto a different projection. + # This error occurs more often if nans are present. + # Replacing nans with scipy.interpolate.griddata solves some of these errors, but not all. + # Interestingly, plt.contour still seems to work in these situations as a (suboptimal) workaround. + # This bug occurs with Shapely 1.8.0, Cartopy 0.20.1. + # This issue was resolved with Cartopy 0.20.2 (https://github.com/SciTools/cartopy/issues/1936). + # Using the transform_first argument results in a noticeable speed-up: + # https://scitools.org.uk/cartopy/docs/latest/gallery/scalar_data/contour_transforms.html + # - This also requires the X and Y variables to be 2-D arrays. + # - This also resolves TopologyException: side location conflict errors that can occur when NaNs are present. + + # If the variable has the same shape as lats, then plot the filled contour field + if fill_var.shape == lats.shape: + contourf = True + plt.contourf(wrf.to_np(lons), wrf.to_np(lats), wrf.to_np(fill_var), bounds, + cmap=cmap, norm=norm, extend=extend, transform=data_crs, transform_first=(ax, True)) + # Otherwise, we presumably need to plot an empty map and then plot markers + else: + contourf = False + if mark1_lon is None or mark1_lat is None: + print('ERROR: map_plot in map_funcs.py:') + print(' var does not match the shape of lons or lats.') + print(' If plotting a contour map, fix the mismatch in shape between fill_var and lons/lats.') + print(' Otherwise, this could imply a need for a blank map with markers.') + print(' However, mark1_lon and/or mark1_lat are not provided either.') + print(' If a filled contour map is desired, ensure fill_var is the same shape as lons and lats.') + print(' If a blank map with markers is desired instead, provide mark1_lon and mark1_lat.') + print(' If the markers should be filled according to a data value, then set mark1_val_fill=True.') + print(' Or if a future use case is identified that requires a blank map & no markers, then modify code.') + print(' Exiting!') + sys.exit() + + # Optional: Add marker set 1 to the plot + if mark1_lon is not None and mark1_lat is not None: + if mark1_lat.shape != mark1_lon.shape: + print('ERROR: map_plot in map_funcs.py:') + print( 'mark1_lat and mark1_lon do not have the same shape.') + print(' Exiting!') + sys.exit() + if lg_text is None: + lg_lab1 = None + else: + lg_lab1 = lg_text[0] + if mark1_val_fill: + # Fill markers according to their data value, cmap, and norm + # NOTE: If plotting markers over a blank map, must use plt.scatter, not ax.scatter, to avoid an error about + # the lack of a mappable when drawing the colorbar below. + if mark1_var.shape == mark1_lat.shape: + if contourf: + ax.scatter(mark1_lon, mark1_lat, c=mark1_var, marker=mark1_style, s=mark1_size, + edgecolors=mark1_edgecolor, label=lg_lab1, linewidths=mark1_width, + cmap=cmap, norm=norm, transform=data_crs, zorder=mark1_zorder) + else: + plt.scatter(mark1_lon, mark1_lat, c=mark1_var, marker=mark1_style, s=mark1_size, + edgecolors=mark1_edgecolor, label=lg_lab1, linewidths=mark1_width, + cmap=cmap, norm=norm, transform=data_crs, zorder=mark1_zorder) + else: + print('ERROR: map_plot in map_funcs.py:') + print(' mark1_var, mark1_lat, and mark1_lon are not all the same shape.') + print(' Exiting!') + sys.exit() + # Marker set 1 not filled by any data values + else: + ax.scatter(mark1_lon, mark1_lat, marker=mark1_style, s=mark1_size, color=mark1_color, + edgecolors=mark1_edgecolor, label=lg_lab1, linewidths=mark1_width, + transform=data_crs, zorder=mark1_zorder) + + # Optional: Add marker set 2 to the plot + if mark2_lon is not None and mark2_lat is not None: + if mark2_lat.shape != mark2_lon.shape: + print('ERROR: map_plot in map_funcs.py:') + print( 'mark2_lat and mark2_lon do not have the same shape.') + print(' Exiting!') + sys.exit() + if lg_text is None: + lg_lab2 = None + else: + lg_lab2 = lg_text[1] + # Marker set 2 filled by data values + if mark2_val_fill: + # Fill markers according to their data value, cmap, and norm + # NOTE: If plotting markers over a blank map, must use plt.scatter, not ax.scatter, to avoid an error about + # the lack of a mappable when drawing the colorbar below. + if mark2_var.shape == mark2_lat.shape: + if contourf: + ax.scatter(mark2_lon, mark2_lat, c=mark2_var, marker=mark2_style, s=mark2_size, + edgecolors=mark2_edgecolor, label=lg_lab2, linewidths=mark2_width, + cmap=cmap, norm=norm, transform=data_crs, zorder=mark2_zorder) + else: + plt.scatter(mark2_lon, mark2_lat, c=mark2_var, marker=mark2_style, s=mark2_size, + edgecolors=mark2_edgecolor, label=lg_lab2, linewidths=mark2_width, + cmap=cmap, norm=norm, transform=data_crs, zorder=mark2_zorder) + else: + print('ERROR: map_plot in map_funcs.py:') + print(' mark2_var, mark2_lat, and mark2_lon are not all the same shape.') + print(' Exiting!') + sys.exit() + # Marker set 2 not filled by data values + else: + ax.scatter(mark2_lon, mark2_lat, marker=mark2_style, s=mark2_size, color=mark2_color, + edgecolors=mark2_edgecolor, label=lg_lab2, linewidths=mark2_width, + transform=data_crs, zorder=mark1_zorder) + + # Draw the colorbar + # Credit: https://stackoverflow.com/questions/30030328/correct-placement-of-colorbar-relative-to-geo-axes-cartopy + # Create colorbar axes (temporarily) anywhere + cax = fig.add_axes([0, 0, 0.1, 0.1]) + # Find the location of the main plot axes + posn = ax.get_position() + # Adjust the positioning and orientation of the colorbar, and then draw it + # The colorbar will inherit the norm/extend attributes from the plt.contourf or plt.scatter call above + if cbar_loc == 'bottom': + cax.set_position([posn.x0, posn.y0-0.09, posn.width, 0.05]) + plt.colorbar(cax=cax, orientation='horizontal', label=cbar_lab) + elif cbar_loc == 'right': + cax.set_position([posn.x0+posn.width+0.05, posn.y0, 0.04, posn.height]) + plt.colorbar(cax=cax, orientation='vertical', label=cbar_lab) + elif cbar_loc == 'top' or cbar_loc == 'left': + print('WARNING: cbar_loc=' + cbar_loc + ' requested. Unsupported option. Colorbar will not be drawn.') + print(' Add directives in map_funcs.map_plot to handle that option and draw the colorbar.') + + # Add the overall plot title + plt.suptitle(suptitle, y=suptitle_y) + + # Optional: Add titles to the subplot + if title_l is not None: + ax.set_title(title_l, fontsize=fontsize, loc='left') + if title_r is not None: + ax.set_title(title_r, fontsize=fontsize, loc='right') + if title_c is not None: + ax.set_title(title_c, fontsize=fontsize, loc='center') + + # Optional: Add set 1 of text labels to the plot + if text1_lab is not None and text1_lat is not None and text1_lon is not None: + if len(text1_lab) != len(text1_lat) or len(text1_lab) != len(text1_lon): + print('ERROR: map_plot in map_funcs.py:') + print(' text1_lab, text1_lat, and text1_lon do not all have the same length.') + print(' Exiting!') + sys.exit() + n_text = len(text1_lab) + for xx in range(n_text): + ax.text(text1_lon[xx], text1_lat[xx], text1_lab[xx], horizontalalignment='center', + transform=data_crs, size=fontsize, zorder=13, weight=text1_lab_wt) + + # Optional: Add set 2 of text labels to the plot + if text2_lab is not None and text2_lat is not None and text2_lon is not None: + if len(text2_lab) != len(text2_lat) or len(text2_lab) != len(text2_lon): + print('ERROR: map_plot in map_funcs.py:') + print(' text2_lab, text2_lat, and text2_lon do not all have the same length.') + print(' Exiting!') + sys.exit() + n_text = len(text2_lab) + for xx in range(n_text): + ax.text(text2_lon[xx], text2_lat[xx], text2_lab[xx], horizontalalignment='center', + transform=data_crs, size=fontsize, zorder=14, weight=text2_lab_wt) + + # Optional: Draw wind barbs + if u is not None and v is not None: + if isinstance(lons, np.ndarray): + x_thin = lons[::map_y_thin, ::map_x_thin] + else: + x_thin = lons[::map_y_thin, ::map_x_thin].values + if isinstance(lats, np.ndarray): + y_thin = lats[::map_y_thin, ::map_x_thin] + else: + y_thin = lats[::map_y_thin, ::map_x_thin].values + u_thin = u[::map_y_thin, ::map_x_thin] + v_thin = v[::map_y_thin, ::map_x_thin] + # Assume winds input to here are in m/s instead of kts, so reduce the barb_increments from 5/10/50 to 2.5/5/25 + ax.barbs(x_thin, y_thin, u_thin, v_thin, length=5, transform=data_crs, linewidth=barb_width, + barb_increments={'half': 2.5, 'full': 5, 'flag': 25}) + + # Optional: Add a legend (most useful if 2+ sets of markers) + if lg_text is not None: + ax.legend(loc=lg_loc, fontsize=lg_fontsize).set_zorder(15) + + # create output directory if it does not already exist + os.makedirs(os.path.dirname(fname), exist_ok=True) + + # Save and close the figure + plt.savefig(fname) + plt.close() diff --git a/use_cases/Hurricane_Matthew/Visualization/plot_wrf.py b/use_cases/Hurricane_Matthew/Visualization/plot_wrf.py new file mode 100755 index 0000000..f4e9b90 --- /dev/null +++ b/use_cases/Hurricane_Matthew/Visualization/plot_wrf.py @@ -0,0 +1,703 @@ +#! /usr/bin/env python3 + +""" +plot_wrf.py + +Written by: Jared A. Lee (jaredlee@ucar.edu) +Written on: 13 May 2024 +""" + +import sys +import argparse +import pathlib +import warnings +import datetime as dt +import numpy as np +import pandas as pd +import xarray as xr +import netCDF4 +import wrf +import matplotlib as mpl + +# Import functions from a local file +import map_funcs + +# def main(init_dt_first, init_dt_last, init_stride_h, plot_beg_lead_time, plot_end_lead_time, plot_stride, domain, exp_name): +def main(script_config_opts): + # ============== + # USER SETTINGS: + # ============== + + # Default plot type selection + plot_type = 'png' + # plot_maps = True # Plot 2D maps + plot_subdomain = False # Plot specified subdomain to zoom in on a defined area of interest + plot_stations = True # Plot station markers & labels on the map (e.g., cities of interest) + + # Which variables should be plotted? + plot_TERRAIN = True # terrain height [m] + plot_T2 = True # 2-m temperature [C] + plot_RH2 = True # 2-m relative humidity [%] + plot_SLP = True # sea level pressure [hPa] (NOTE: this requires several other variables in the wrfout file) + plot_WS10 = True # 10-m wind speed [m s-1] + plot_REFL = True # simulated radar reflectivity [dBZ] + plot_RAIN = True # total accumulated rainfall during the simulation [mm] + plot_WS100 = True # 100-m wind speed [m s-1] + plot_GHT500 = True # 500-hPa geopotential height [m] + + # Plot any overlays, like wind barbs? + plot_wind_barbs_sfc = True # overlay 10-m wind barbs for selected plots + plot_wind_barbs_upr = True # overlay upper-air wind barbs for selected plots + + # Default water color (generally use only in terrain plots) + water_color = 'lightblue' + + # Set some other plot options + suptitle_y = 1.00 + plot_fontsize = 13 + barb_thin = 10 + barb_width = 0.5 + + # Set some text labels for demonstration + if plot_stations: + text1_lab = ['Miami', 'Jacksonville', 'Charleston'] + mark1_lat = np.asarray([25.7617, 30.3322, 32.7833]) + mark1_lon = np.asarray([-80.1918, -81.6557, -79.9320]) + text1_lat = np.asarray(mark1_lat) + np.asarray([-0.20, -0.20, -0.40]) + text1_lon = np.asarray(mark1_lon) + np.asarray([1.50, 3.00, 2.70]) + mark1_size = 36 + mark1_color = 'black' + + # Domain plotting ranges in (i,j) space (whole domain by default) + i_beg, i_end = 0, -1 + j_beg, j_end = 0, -1 + # TODO: Implement subdomain plotting + if plot_subdomain: + # Adjust these if you want a different subdomain + i_beg, i_end = 10, 81 + j_beg, j_end = 10, 90 + + lat_labels = [16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40] + lon_labels = [-62, -64, -66, -68, -70, -72, -74, -76, -78, -80, -82, -84, -86, -88] + + # Set map contour plot limits + min_terrain = 0.0 + max_terrain = 1500.1 + int_terrain = 100.0 + + min_slp = 980.0 + max_slp = 1020.1 + int_slp = 2.0 + + min_t2 = 0.0 + max_t2 = 40.1 + int_t2 = 2.0 + + min_rh2 = 0.0 + max_rh2 = 100.1 + int_rh2 = 5.0 + + min_ws10 = 0.0 + max_ws10 = 35.0 + int_ws10 = 2.5 + + min_rain = 0.0 + max_rain = 100.1 + int_rain = 5.0 + + # ======================================= + # CONSTANTS, FORMAT STATEMENTS, AND MORE: + # ======================================= + + G = 9.81 # graviational acceleration [m s-2] + PI = 3.1415926 + DEG2RAD = PI / 180.0 + RAD2DEG = 180.0 / PI + Rd = 297.048 # specific gas constant for dry air [J kg-1 K-1] + Rv = 461.495 # specific gas constant for water vapor [J kg-1 K-1] + C_to_K = 273.15 # additive conversion between degrees Celsius and Kelvin + + missing_val = -9999.0 + + mpl_Wm2 = 'W $\mathregular{m^{-2}}$' + mpl_ms1 = 'm $\mathregular{s^{-1}}$' + mpl_s1 = '$\mathregular{s^{-1}}$' + mpl_Jkg = 'J $\mathregular{kg^{-1}}$' + mpl_um = u'\u03bcm' + mpl_gkg1 = 'g $\mathregular{kg^{-1}}$' + mpl_kgkg1 = 'kg $\mathregular{kg^{-1}}$' + mpl_gm2s1 = 'g $\mathregular{m^{-2}}$ $\mathregular{s^{-1}}$' + mpl_kgm2s1 = 'kg $\mathregular{m^{-2}} $\mathregular{s^{-1}}$' + mpl_kgm2 = 'kg $\mathregular{m^{-2}}$' + mpl_10m3 = '$\mathregular{10^{-3}}$' + mpl_10m4 = '$\mathregular{10^{-4}}$' + mpl_10m5 = '$\mathregular{10^{-5}}$' + mpl_10m6 = '$\mathregular{10^{-6}}$' + + deg_uni = '\u00B0' + en_dash = u'\u2013' + em_dash = u'\u2014' + + fmt_yyyymmdd = '%Y%m%d' + fmt_yyyymmddhh = '%Y%m%d%H' + fmt_yyyymmdd_hh = '%Y%m%d_%H' + fmt_yyyymmdd_hhmm = '%Y%m%d_%H%M' + fmt_dt = '%Y%m%dT%H%M%S' + fmt_yyyy = '%Y' + fmt_mm = '%m' + fmt_dd = '%d' + fmt_hh = '%H' + fmt_nn = '%M' + + fmt_wrf_dt_no_s = '%Y-%m-%d_%H:%M' + fmt_wrf_date = '%Y-%m-%d' + fmt_wrf_time = '%H:%M:%S' + fmt_wrf_dt = fmt_wrf_date + '_' + fmt_wrf_time + fmt_time_file = fmt_yyyymmdd_hhmm + fmt_time_plot = '%d %b %Y/%H%M UTC' + + # Define a custom colormap for radar reflectivity plots + # Modified to add gray for 0–5 dBZ and lightpurple for 75+ dBZ + cmap_radar = np.array([ + [200, 200, 200], [4, 233, 231], [1, 159, 244], [3, 0, 244], + [2, 253, 2], [1, 197, 1], [0, 142, 0], + [253, 248, 2], [229, 188, 0], [253, 149, 0], + [253, 0, 0], [212, 0, 0], [188, 0, 0], + [248, 0, 253], [152, 84, 198], [228, 199, 243]], np.float32) / 255.0 + bounds_radar = np.arange(0., 75.01, 5.0) + # Color names are approximate and only intended for assistance deciphering the RGB table above + colors_radar = np.array([ + 'gray', 'cyan', 'lightblue', 'darkblue', + 'lightgreen', 'green', 'darkgreen', + 'yellow', 'lightorange', 'orange', + 'red', 'darkred', 'brickred', + 'fuschia', 'violet', 'lavender']) + + read_zlev = False + read_plev = False + if plot_WS100: + read_zlev = True + if plot_GHT500: + read_plev = True + + # ============= + # MAIN PROGRAM: + # ============= + + cycle_dt_str_first = script_config_opts['cycle_dt_first'] + cycle_dt_str_last = script_config_opts['cycle_dt_last'] + cycle_stride_h = script_config_opts['cycle_stride_h'] + cycle_dt_first = pd.to_datetime(cycle_dt_str_first, format=fmt_yyyymmdd_hh) + cycle_dt_last = pd.to_datetime(cycle_dt_str_last, format=fmt_yyyymmdd_hh) + cycle_dt_all = pd.date_range(start=cycle_dt_first, end=cycle_dt_last, freq=str(cycle_stride_h) + 'h') + n_cycles = len(cycle_dt_all) + + dom_num = script_config_opts['domain'] + wrf_dom = 'd0' + dom_num + + # Loop over forecast cycles/initializations + for cc in range(n_cycles): + cycle_dt = cycle_dt_all[cc] + cycle_dt_str = cycle_dt.strftime(fmt_yyyymmdd_hh) + cycle_dt_plot = cycle_dt.strftime(fmt_time_plot) + start_time_plot = 'Start: ' + cycle_dt_plot + wrf_dir = script_config_opts['wrf_dir_parent'].joinpath(cycle_dt_str, 'wrfout') + out_dir = script_config_opts['out_dir_parent'].joinpath(cycle_dt_str, 'plots') + + # Build an array of the valid datetimes that need to be read and plotted + beg_lead_h = int(script_config_opts['beg_lead_time'].split(':')[0]) + beg_lead_m = int(script_config_opts['beg_lead_time'].split(':')[1]) + end_lead_h = int(script_config_opts['end_lead_time'].split(':')[0]) + end_lead_m = int(script_config_opts['end_lead_time'].split(':')[1]) + valid_dt_beg = cycle_dt + dt.timedelta(hours=beg_lead_h, minutes=beg_lead_m) + valid_dt_end = cycle_dt + dt.timedelta(hours=end_lead_h, minutes=end_lead_m) + valid_dt_all = pd.date_range(start=valid_dt_beg, end=valid_dt_end, + freq=str(script_config_opts['str_lead_time']) + 'min') + n_valid = len(valid_dt_all) + + # Loop over valid times (this assumes one output time per file) + for vv in range(n_valid): + valid_dt = valid_dt_all[vv] + valid_dt_wrf = valid_dt.strftime(fmt_wrf_dt) + valid_dt_plot = valid_dt.strftime(fmt_time_plot) + valid_dt_file = valid_dt.strftime(fmt_time_file) + valid_time_plot = 'Valid: '+valid_dt_plot + + # suptitle = 'Hurricane Matthew ' + em_dash + ' Domain ' + dom_num + suptitle = 'Hurricane Matthew Test Case' + map_prefix = 'map_wrf_' + wrf_dom + '_' + map_suffix = '_' + valid_dt_file + '.' + plot_type + title_r = start_time_plot + '\n' + valid_time_plot + title_r_no_valid = start_time_plot + title_r_blank = '' + + wrf_fname = wrf_dir.joinpath('wrfout_' + wrf_dom + '_' + valid_dt_wrf) + wrf_fname_zlev = wrf_dir.joinpath('wrfout_zlev_' + wrf_dom + '_' + valid_dt_wrf) + wrf_fname_plev = wrf_dir.joinpath('wrfout_plev_' + wrf_dom + '_' + valid_dt_wrf) + + try: + wrf_fname.is_file() + except FileNotFoundError: + print('WARNING: File '+str(wrf_fname) + 'does not exist. Continuing to the next valid time.') + continue + print('Reading ' + str(wrf_fname)) + # Use NetCDF4-python to open a Dataset, as wrf-python doesn't yet take an xarray Dataset + # wrf.getvar will return an xarray Dataset by default, though + ds_wrf_nc = netCDF4.Dataset(wrf_fname, mode='r') + + # Static fields only need to be read in once + if cc == 0 and vv == 0: + # Latitude, Longitude + da_lat = wrf.getvar(ds_wrf_nc, 'lat', squeeze=False) + da_lon = wrf.getvar(ds_wrf_nc, 'lon', squeeze=False) + wrf_lat = da_lat.values[0, :, :] + wrf_lon = da_lon.values[0, :, :] + n_wrf_lat = wrf_lat.shape[0] + n_wrf_lon = wrf_lon.shape[1] + n_wrf_lev = int(getattr(ds_wrf_nc, 'BOTTOM-TOP_GRID_DIMENSION')) + wrf_lats, wrf_lons = wrf.latlon_coords(da_lat) + + print('Getting cartopy mapping objects') + cart_proj = wrf.get_cartopy(wrfin=ds_wrf_nc) + cart_bounds = wrf.geo_bounds(var=da_lat[0, j_beg:j_end, i_beg:i_end]) + cart_xlim = wrf.cartopy_xlim(wrfin=ds_wrf_nc, geobounds=cart_bounds) + cart_ylim = wrf.cartopy_ylim(wrfin=ds_wrf_nc, geobounds=cart_bounds) + borders, states, oceans, lakes, rivers, land = map_funcs.get_cartopy_features() + + # Start populating dictionary for map plotting options. Update later with other options. + map_opts = { + 'cart_proj': cart_proj, 'cart_xlim': cart_xlim, 'cart_ylim': cart_ylim, + 'borders': borders, 'states': states, 'oceans': oceans, 'lakes': lakes, + 'lons': wrf_lons, 'lats': wrf_lats, 'suptitle': suptitle, 'suptitle_y': suptitle_y, + 'lat_labels': lat_labels, 'lon_labels': lon_labels, 'fontsize': plot_fontsize, + 'map_x_thin': barb_thin, 'map_y_thin': barb_thin, 'barb_width': barb_width, + } + + if plot_stations: + map_opts['mark1_lat'] = mark1_lat + map_opts['mark1_lon'] = mark1_lon + map_opts['text1_lab'] = text1_lab + map_opts['text1_lat'] = text1_lat + map_opts['text1_lon'] = text1_lon + map_opts['mark1_size'] = mark1_size + map_opts['mark1_color'] = mark1_color + + # Terrain + if plot_TERRAIN: + print(' Reading terrain') + da_terrain = wrf.getvar(ds_wrf_nc, 'ter', squeeze=False) + wrf_terrain = da_terrain.values[0, :, :] + + var_file = 'TERRAIN' + var_name = 'Terrain Height' + var_unit = 'm' + wrf_var = wrf_terrain + min_val = np.nanmin(wrf_var[j_beg:j_end, i_beg:i_end]) + max_val = np.nanmax(wrf_var[j_beg:j_end, i_beg:i_end]) + extend = 'both' + cmap = map_funcs.truncate_cmap(mpl.cm.terrain, minval=0.20, maxval=0.95) + bounds = np.arange(min_terrain, max_terrain, int_terrain) + norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend) + title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit + map_opts['fill_var'] = wrf_var + map_opts['water_color'] = water_color + map_opts['extend'] = extend + map_opts['cmap'] = cmap + map_opts['bounds'] = bounds + map_opts['norm'] = norm + map_opts['cbar_lab'] = 'Model ' + var_name + ' [' + var_unit + ']' + map_opts['fname'] = out_dir.joinpath(map_prefix + var_file + '.' + plot_type) + map_opts['title_l'] = title_l + map_opts['title_r'] = title_r_blank + map_funcs.map_plot(map_opts) + + # Make the water color transparent for all subsequent plots + map_opts['water_color'] = 'none' + map_opts['title_r'] = title_r + + if plot_wind_barbs_sfc or plot_WS10: + print(' Reading 10-m wind components (rotated to earth-relative)') + da_uv10 = wrf.getvar(ds_wrf_nc, 'uvmet10', squeeze=False) + wrf_u10 = da_uv10.values[0, 0, :, :] + wrf_v10 = da_uv10.values[1, 0, :, :] + wrf_ws10 = np.sqrt(wrf_u10**2 + wrf_v10**2) + + if plot_WS10: + var_file = 'WS10' + var_name = '10-m Wind Speed' + var_unit = mpl_ms1 + wrf_var = wrf_ws10 + min_val = np.nanmin(wrf_var[j_beg:j_end, i_beg:i_end]) + max_val = np.nanmax(wrf_var[j_beg:j_end, i_beg:i_end]) + extend = 'max' + cmap = mpl.cm.BuGn + bounds = np.arange(min_ws10, max_ws10, int_ws10) + norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend) + map_opts['cbar_lab'] = var_name + ' [' + var_unit + ']' + if plot_wind_barbs_sfc: + var_file = var_file + '+barbs' + # var_name = var_name + '; Barbs' + map_opts['u'] = wrf_u10 + map_opts['v'] = wrf_v10 + else: + map_opts['u'] = None + map_opts['v'] = None + title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit + map_opts['fill_var'] = wrf_var + map_opts['extend'] = extend + map_opts['cmap'] = cmap + map_opts['bounds'] = bounds + map_opts['norm'] = norm + map_opts['fname'] = out_dir.joinpath(map_prefix + var_file + map_suffix) + map_opts['title_l'] = title_l + map_funcs.map_plot(map_opts) + + # Sea level pressuure + if plot_SLP: + print(' Reading sea level pressure') + da_slp = wrf.getvar(ds_wrf_nc, 'slp', squeeze=False) + wrf_slp = da_slp.values[0, :, :] + + var_file = 'SLP' + var_name = 'Sea-Level Pressure' + var_unit = 'hPa' + wrf_var = wrf_slp + min_val = np.nanmin(wrf_var[j_beg:j_end, i_beg:i_end]) + max_val = np.nanmax(wrf_var[j_beg:j_end, i_beg:i_end]) + extend = 'both' + cmap = mpl.cm.viridis + bounds = np.arange(min_slp, max_slp, int_slp) + norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend) + map_opts['cbar_lab'] = var_name + ' [' + var_unit + ']' + if plot_wind_barbs_sfc: + var_file = var_file + '+barbs' + # var_name = var_name + '; 10-m Barbs' + map_opts['u'] = wrf_u10 + map_opts['v'] = wrf_v10 + else: + map_opts['u'] = None + map_opts['v'] = None + title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit + map_opts['fill_var'] = wrf_var + map_opts['extend'] = extend + map_opts['cmap'] = cmap + map_opts['bounds'] = bounds + map_opts['norm'] = norm + map_opts['fname'] = out_dir.joinpath(map_prefix + var_file + map_suffix) + map_opts['title_l'] = title_l + map_funcs.map_plot(map_opts) + + # 2-m air temperature + if plot_T2: + print(' Reading 2-m air temperature') + da_t2 = wrf.getvar(ds_wrf_nc, 'T2', squeeze=False) + if da_t2.attrs['units'] == 'K': + da_t2 = da_t2 - C_to_K + da_t2.attrs['units'] = 'degC' + wrf_t2 = da_t2.values[0, :, :] + + var_file = 'T2' + var_name = '2-m Air Temperature' + var_unit = deg_uni + 'C' + wrf_var = wrf_t2 + min_val = np.nanmin(wrf_var[j_beg:j_end, i_beg:i_end]) + max_val = np.nanmax(wrf_var[j_beg:j_end, i_beg:i_end]) + extend = 'both' + cmap = mpl.cm.rainbow + bounds = np.arange(min_t2, max_t2, int_t2) + norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend) + map_opts['cbar_lab'] = var_name + ' [' + var_unit + ']' + if plot_wind_barbs_sfc: + var_file = var_file + '+barbs' + # var_name = var_name + '; 10-m Barbs' + map_opts['u'] = wrf_u10 + map_opts['v'] = wrf_v10 + else: + map_opts['u'] = None + map_opts['v'] = None + title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit + map_opts['fill_var'] = wrf_var + map_opts['extend'] = extend + map_opts['cmap'] = cmap + map_opts['bounds'] = bounds + map_opts['norm'] = norm + map_opts['fname'] = out_dir.joinpath(map_prefix + var_file + map_suffix) + map_opts['title_l'] = title_l + map_funcs.map_plot(map_opts) + + # 2-m relative humidity + if plot_RH2: + print(' Reading 2-m relative humidity') + da_rh2 = wrf.getvar(ds_wrf_nc, 'rh2', squeeze=False) + wrf_rh2 = da_rh2.values[0, :, :] + + var_file = 'RH2' + var_name = '2-m Relative Humidity' + var_unit = '%' + wrf_var = wrf_rh2 + min_val = np.nanmin(wrf_var[j_beg:j_end, i_beg:i_end]) + max_val = np.nanmax(wrf_var[j_beg:j_end, i_beg:i_end]) + extend = 'max' + cmap = mpl.cm.YlGnBu + bounds = np.arange(min_rh2, max_rh2, int_rh2) + norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend) + map_opts['cbar_lab'] = var_name + ' [' + var_unit + ']' + if plot_wind_barbs_sfc: + var_file = var_file + '+barbs' + # var_name = var_name + '; 10-m Barbs' + map_opts['u'] = wrf_u10 + map_opts['v'] = wrf_v10 + else: + map_opts['u'] = None + map_opts['v'] = None + title_l = var_name + f'\nMin: {min_val:.1f}' + var_unit + f', Max: {max_val:.1f}' + var_unit + map_opts['fill_var'] = wrf_var + map_opts['extend'] = extend + map_opts['cmap'] = cmap + map_opts['bounds'] = bounds + map_opts['norm'] = norm + map_opts['fname'] = out_dir.joinpath(map_prefix + var_file + map_suffix) + map_opts['title_l'] = title_l + map_funcs.map_plot(map_opts) + + # Accumulated rainfall + if plot_RAIN and vv > 0: + print(' Reading accumulated rainfall') + da_rainc = wrf.getvar(ds_wrf_nc, 'RAINC', squeeze=False) + da_rainnc = wrf.getvar(ds_wrf_nc, 'RAINNC', squeeze=False) + wrf_rain = da_rainc.values[0, :, :] + da_rainnc.values[0, :, :] + # Mask RAIN=0.0 for plotting + wrf_rain_plot = np.ma.masked_equal(np.where(wrf_rain == 0.0, missing_val, wrf_rain), missing_val) + + var_file = 'RAIN' + var_name = 'Accumulated Precipitation' + var_unit = 'mm' + wrf_var1 = wrf_rain + wrf_var2 = wrf_rain_plot + min_val = np.nanmin(wrf_var1[j_beg:j_end, i_beg:i_end]) + max_val = np.nanmax(wrf_var1[j_beg:j_end, i_beg:i_end]) + extend = 'max' + cmap = mpl.cm.GnBu + bounds = np.arange(min_rain, max_rain, int_rain) + norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend) + map_opts['cbar_lab'] = var_name + ' [' + var_unit + ']' + title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit + map_opts['fill_var'] = wrf_var2 + map_opts['extend'] = extend + map_opts['cmap'] = cmap + map_opts['bounds'] = bounds + map_opts['norm'] = norm + map_opts['fname'] = out_dir.joinpath(map_prefix + var_file + map_suffix) + map_opts['title_l'] = title_l + # print(wrf_var1.shape) + # print(wrf_var2.shape) + map_funcs.map_plot(map_opts) + + # Radar reflectivity + if plot_REFL and vv > 0: + print(' Reading radar reflectivity') + da_refl = wrf.getvar(ds_wrf_nc, 'dbz', squeeze=False) + wrf_refl = da_refl.values[0, 0, :, :] + # Mask REFL <= 0.0 for plotting + wrf_refl_plot = np.ma.masked_equal(np.where(wrf_refl <= 0.0, missing_val, wrf_refl), missing_val) + + var_file = 'REFL' + var_name = 'Radar Reflectivity' + var_unit = 'dBZ' + wrf_var1 = wrf_refl + wrf_var2 = wrf_refl_plot + min_val = np.nanmin(wrf_var1[j_beg:j_end, i_beg:i_end]) + max_val = np.nanmax(wrf_var1[j_beg:j_end, i_beg:i_end]) + extend = 'max' + refl_rgb = cmap_radar + bounds = bounds_radar + cmap, norm = mpl.colors.from_levels_and_colors(bounds, refl_rgb, extend=extend) + map_opts['cbar_lab'] = var_name + ' [' + var_unit + ']' + if plot_wind_barbs_sfc: + var_file = var_file + '+barbs' + var_name = var_name + '; 10-m Barbs' + map_opts['u'] = wrf_u10 + map_opts['v'] = wrf_v10 + else: + map_opts['u'] = None + map_opts['v'] = None + title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit + map_opts['fill_var'] = wrf_var2 + map_opts['extend'] = extend + map_opts['cmap'] = cmap + map_opts['bounds'] = bounds + map_opts['norm'] = norm + map_opts['fname'] = out_dir.joinpath(map_prefix + var_file + map_suffix) + map_opts['title_l'] = title_l + # print(wrf_var1.shape) + # print(wrf_var2.shape) + map_funcs.map_plot(map_opts) + + if read_zlev: + try: + wrf_fname_zlev.is_file() + except FileNotFoundError: + print('WARNING: File ' + str(wrf_fname_zlev) + 'does not exist. Continuing to the next valid time.') + continue + print('Reading ' + str(wrf_fname_zlev)) + # Use NetCDF4-python to open a Dataset, as wrf-python doesn't yet take an xarray Dataset + # wrf.getvar will return an xarray Dataset by default, though + ds_wrf_zlev_nc = netCDF4.Dataset(wrf_fname_zlev, mode='r') + wrf_z_zlev = wrf.getvar(ds_wrf_zlev_nc, 'Z_ZL', squeeze=False) + + # 100-m wind speed + if plot_WS100: + ind_z = np.where(wrf_z_zlev == -100)[0][0] + wrf_ws100 = wrf.getvar(ds_wrf_zlev_nc, 'S_ZL', squeeze=False)[0, ind_z, :, :] + + var_file = 'WS100' + var_name = '100-m Wind Speed' + var_unit = mpl_ms1 + wrf_var = wrf_ws100 + min_val = np.nanmin(wrf_var[j_beg:j_end, i_beg:i_end]) + max_val = np.nanmax(wrf_var[j_beg:j_end, i_beg:i_end]) + extend = 'max' + cmap = mpl.cm.BuGn + bounds = np.arange(min_ws10, max_ws10, int_ws10) + norm = mpl.colors.BoundaryNorm(bounds, cmap.N, extend=extend) + map_opts['cbar_lab'] = var_name + ' [' + var_unit + ']' + if plot_wind_barbs_upr: + var_file = var_file + '+barbs' + var_name = var_name + '; Barbs' + wrf_u100 = wrf.getvar(ds_wrf_zlev_nc, 'U_ZL', squeeze=False).values[0, ind_z, :, :] + wrf_v100 = wrf.getvar(ds_wrf_zlev_nc, 'V_ZL', squeeze=False).values[0, ind_z, :, :] + map_opts['u'] = wrf_u100 + map_opts['v'] = wrf_v100 + else: + map_opts['u'] = None + map_opts['v'] = None + title_l = var_name + f'\nMin: {min_val:.1f} ' + var_unit + f', Max: {max_val:.1f} ' + var_unit + map_opts['fill_var'] = wrf_var + map_opts['extend'] = extend + map_opts['cmap'] = cmap + map_opts['bounds'] = bounds + map_opts['norm'] = norm + map_opts['fname'] = out_dir.joinpath(map_prefix + var_file + map_suffix) + map_opts['title_l'] = title_l + map_funcs.map_plot(map_opts) + + + +def parse_args(): + parser = argparse.ArgumentParser() + parser.add_argument('-w', '--wrf_dir_parent', default='/data/input/wrf', + help='string specifying the directory path to the parent WRF output directories, ' + 'above any experiment or cycle datetime subdirectories (default: /data/input/wrf)') + parser.add_argument('-o', '--out_dir_parent', default='/data/output/wrf', + help='string specifying the directory path to the parent plot directories (default: /data/output/wrf)') + parser.add_argument('-f', '--cycle_dt_first', default='20161006_00', + help='beginning date/time of first WRF simulation [YYYYMMDD_HH] (default: 20161006_00)') + parser.add_argument('-l', '--cycle_dt_last', default=None, + help='beginning date/time of last WRF simulation [YYYYMMDD_HH]') + parser.add_argument('-i', '--cycle_stride_h', default=24, type=int, + help='stride in hours between cycles (default: 24)') + parser.add_argument('-b', '--beg_lead_time', default='00:00', + help='beginning lead time for plotting WRF simulations [HH:MM] (default: 00:00)') + parser.add_argument('-e', '--end_lead_time', default='48:00', + help='ending lead time for plotting WRF simulations [HH:MM] (default: 48:00)') + parser.add_argument('-s', '--str_lead_time', default=180, type=int, + help='stride to create plots every N minutes (default: 180)') + parser.add_argument('-d', '--domain', default='1', help='WRF domain number to be plotted (default: 1)') + # parser.add_argument('-x', '--exp_name', default=None, + # help='WRF experiment name(s), if applicable. If requesting plots for multiple experiments, ' + # 'separate them by commas (e.g., exp01,exp02).') + + args = parser.parse_args() + wrf_dir_parent = args.wrf_dir_parent + out_dir_parent = args.out_dir_parent + cycle_dt_first = args.cycle_dt_first + cycle_dt_last = args.cycle_dt_last + cycle_stride_h = args.cycle_stride_h + beg_lead_time = args.beg_lead_time + end_lead_time = args.end_lead_time + str_lead_time = args.str_lead_time + domain = args.domain + # exp_names_inp = args.exp_name + + # if exp_names_inp is None: + # exp_name = None + # else: + # exp_name = exp_names_inp.split(',') + + if out_dir_parent is None: + out_dir_parent = wrf_dir_parent + + # Make both paths into pathlib objects + wrf_dir_parent = pathlib.Path(wrf_dir_parent) + out_dir_parent = pathlib.Path(out_dir_parent) + + if len(cycle_dt_first) != 11: + print('ERROR! Incorrect length for positional argument init_dt_first. Exiting!') + parser.print_help() + sys.exit() + elif cycle_dt_first[8] != '_': + print('ERROR! Incorrect format for positional argument init_dt_first. Exiting!') + parser.print_help() + sys.exit() + + if cycle_dt_last != None: + if len(cycle_dt_last) != 11: + print('ERROR! Incorrect length for positional argument init_dt_last. Exiting!') + parser.print_help() + sys.exit() + elif cycle_dt_last[8] != '_': + print('ERROR! Incorrect format for positional argument init_dt_last. Exiting!') + parser.print_help() + sys.exit() + else: + cycle_dt_last = cycle_dt_first + + if len(beg_lead_time) != 5: + print('ERROR! Incorrect length for optional argument -b (beg_lead_time). Exiting!') + parser.print_help() + sys.exit() + elif beg_lead_time[2] != ':': + print('ERROR! Incorrect format for optional argument -b (beg_lead_time). Exiting!') + parser.print_help() + sys.exit() + + if len(end_lead_time) != 5: + print('ERROR! Incorrect length for optional argument -e (end_lead_time). Exiting!') + parser.print_help() + sys.exit() + elif end_lead_time[2] != ':': + print('ERROR! Incorrect format for optional argument -e (end_lead_time). Exiting!') + parser.print_help() + sys.exit() + + # Put all these configuration options into a dictionary, to make further development or customization easier + script_config_opts = { + 'wrf_dir_parent': wrf_dir_parent, + 'out_dir_parent': out_dir_parent, + 'cycle_dt_first': cycle_dt_first, + 'cycle_dt_last': cycle_dt_last, + 'cycle_stride_h': cycle_stride_h, + 'beg_lead_time': beg_lead_time, + 'end_lead_time': end_lead_time, + 'str_lead_time': str_lead_time, + 'domain': domain, + # 'exp_name': exp_name, + } + + return script_config_opts + # return init_dt_first, init_dt_last, init_stride_h, beg_lead_time, end_lead_time, plot_stride, domain, exp_name + +if __name__ == '__main__': + now_time_beg = dt.datetime.utcnow() + # init_dt_first, init_dt_last, init_stride_h, plot_beg_lead_time, plot_end_lead_time, plot_stride, domain, exp_name = parse_args() + # main(init_dt_first, init_dt_last, init_stride_h, plot_beg_lead_time, plot_end_lead_time, plot_stride, domain, exp_name) + script_config_opts = parse_args() + main(script_config_opts) + now_time_end = dt.datetime.utcnow() + run_time_tot = now_time_end - now_time_beg + now_time_beg_str = now_time_beg.strftime('%Y-%m-%d %H:%M:%S') + now_time_end_str = now_time_end.strftime('%Y-%m-%d %H:%M:%S') + print('\nScript completed successfully.') + print(' Beg time: '+now_time_beg_str) + print(' End time: '+now_time_end_str) + print(' Run time: '+str(run_time_tot)+'\n')