plotting_maps

This repo aims to plot standardised maps of ACS climate hazard data. We will develop and provide examples of mapping climate hazards for Australia so that data can be consistently presented.

Examples will include maps for Australia and selected states or regions.

Intended uses include taking netcdf or xarray dataarrays of hazards and indices such as Rx1day, TXx, FFDI and plotting the data on a map of Australia.

Developed by Gen Tolhurst (gentolhurst@gmail.com), Supervised by Mitchell Black (Mitchell.Black@bom.gov.au).\ Funded by ACS

Python environment

This code is designed to work with hh5 analysis3-24.04 virtual environment.

In your terminal, this may look like:

$ module use /g/data/hh5/public/modules
$ module load conda/analysis3-24.04

When starting a new ARE JupyterLab session (https://are.nci.org.au/pun/sys/dashboard/batch_connect/sys/jupyter/ncigadi/session_contexts/new, requires NCI login), selecting the hh5 analysis3-24.04 virtual environment might look like this:

Access shapefiles

This code references shapefiles stored in /g/data/ia39/. You will need to be a member of this project to access the data. Request membership https://my.nci.org.au/mancini/project/ia39

See https://github.com/aus-ref-clim-data-nci/shapefiles for more information on the shapefiles.

Include the projects you need when you start an ARE session. Eg, storage: "gdata/ia39+gdata/hh5+gdata/mn51"

Cloning this repo

Before you can import acs_plotting_maps to use the plotting function plot_acs_hazard, you will need to clone a copy of this repository to your own working directory.

If you are working in your home directory, navigate there:

$ cd ~/

Else, if you are working elsewhere (eg. scratch or project), specify the path:

$ cd /path/to/dir/

$ cd /scratch/PROJECT/USER/

$ cd /g/data/PROJECT/USER/

Then, you can clone this repository to access the Python code and notebooks. \ If you want the new directory to be called anything other than "plotting_maps" please replace the final argument with your choice of directory name:

$ git clone https://github.com/AusClimateService/plotting_maps.git plotting_maps

You will now be able to access the functions, python scripts, and Jupyter notebooks from your user.

Usage in Jupyter Notebook:

See a small, easy-to-follow example here:

• [https://github.com/AusClimateService/plotting_maps/blob/main/minimal_plotting_example_pr.ipynb]
• [https://github.com/AusClimateService/plotting_maps/blob/main/area_statistics_example.ipynb]

Other examples:

• [https://github.com/AusClimateService/plotting_maps/blob/main/plotting_and_stats_examples.ipynb]
• [https://github.com/AusClimateService/plotting_maps/blob/main/acs_plotting_maps_examples.ipynb]

1. Navigate to the directory you cloned to:

   cd ~/plotting_maps

2. Import the ACS plotting maps function and dictionaries and Xarray.

   from acs_plotting_maps import plot_acs_hazard, regions_dict, cmap_dict, tick_dict, plot_acs_hazard_3pp
   import xarray as xr

3. **Load some data.** For example, this will load extratropical storm rx5day rainfall
```python
ds = xr.open_dataset("/g/data/ia39/australian-climate-service/test-data/CORDEX-CMIP6/bias-adjustment-input/AGCD-05i/BOM/ACCESS-CM2/historical/r4i1p1f1/BARPA-R/v1-r1/day/pr/pr_AGCD-05i_ACCESS-CM2_historical_r4i1p1f1_BOM_BARPA-R_v1-r1_day_19600101-19601231.nc")

This data has three dimensions (time, lon, lat). There is a value for every day from 01-01-1960 to 31-12-1960. We can only plot 2D, so next we will calculate a statistic to summarise the data

4. Summarise data into a 2D xr.DataArray. For example, calculate the annual sum:

   var="pr"
   da = ds.sum(dim="time")[var]

5. Finally, use the plotting function.\ You will need to specify:

   • the data (and select the variable eg "pr");
   • suitable arguments for the colorbar including cmap, ticks, cbar_label, and cbar_extend;
   • annotations including title, dataset_name, date_range; and
   • where you want the image outfile saved.

regions = regions_dict['ncra_regions']
plot_acs_hazard(data = da,
                regions = regions,
                cmap=cmap_dict["pr"],
                ticks=tick_dict['pr_annual'],
                cbar_label="annual rainfall [mm]",
                cbar_extend="max",
                title = "Rainfall",
                watermark="EXPERIMENTAL IMAGE ONLY", # When you are making your final figures, remove watermark using: watermark=""
                dataset_name = ds.source_id,
                date_range=f"1 January 1960 to 31 December 1960",
                outfile = "~/figures/out.png");

Plot a three-panel plot

%%time
from plotting_maps.acs_plotting_maps import plot_acs_hazard_3pp

var = "HWAtx"

ds_gwl12 =xr.open_dataset("/g/data/ia39/ncra/heat/data/HWAtx/bias-corrected/ensemble/GWL-average/HWAtx_AGCD-05i_MME50_ssp370_v1-r1-ACS-QME-AGCD-1960-2022_GWL12.nc")
ds_gwl15 = xr.open_dataset("/g/data/ia39/ncra/heat/data/HWAtx/bias-corrected/ensemble/GWL-average/HWAtx_AGCD-05i_MME50_ssp370_v1-r1-ACS-QME-AGCD-1960-2022_GWL15.nc")
ds_gwl20 = xr.open_dataset("/g/data/ia39/ncra/heat/data/HWAtx/bias-corrected/ensemble/GWL-average/HWAtx_AGCD-05i_MME50_ssp370_v1-r1-ACS-QME-AGCD-1960-2022_GWL20.nc")
ds_gwl30 = xr.open_dataset("/g/data/ia39/ncra/heat/data/HWAtx/bias-corrected/ensemble/GWL-average/HWAtx_AGCD-05i_MME50_ssp370_v1-r1-ACS-QME-AGCD-1960-2022_GWL30.nc")

plot_acs_hazard_3pp(ds_gwl15 = ds_gwl15[var], 
                    ds_gwl20 = ds_gwl20[var],
                    ds_gwl30 = ds_gwl30[var],
                    regions = regions_dict['ncra_regions'],
                    cbar_label=f"Temperature [degC]",
                    title=f"Maximum Temperature of Hottest Heatwave for future warming scenarios", 
                    date_range = "Insert subtitle - should include the date range of the data \nand then the dataset below that", 
                    # baseline = "GWL1.2", 
                    dataset_name= "MME50_ssp370",
                    issued_date=None,
                    watermark="EXPERIMENTAL IMAGE ONLY", # When you are making your final figures, remove watermark using: watermark=""
                    watermark_color="k",
                    cmap = cmap_dict["tasmax"],
                    ticks = np.arange(18,53,2),)

6. Calculate summary statistics for the range of models.

# Import needed packages
from acs_area_statistics import acs_regional_stats, get_regions
regions = get_regions(["ncra_regions", "australia"])

this has changed. Previously

"# import needed packages

from acs_area_statistics import acs_regional_stats, regions"

For Calculating the NCRA region stats, we want to compare the regional averages based on different models, eg what is the regional mean value from the coolest/driest model relisation, what is the mean, what is the regional mean from the hottest/wettest model for this, we want ds to have the 10th, median and 90th percentile values from each model, then we can find the range of the models and the MMM.

# calculate the stats using the acs_region_fractional_stats function
# Find the min, mean, max value for each region

ds = xr.open_dataset(filename)
mask_frac = regions.mask_3D_frac_approx(ds)
dims = ("lat", "lon",)
how = ["min", "mean", "max"]

da_summary = acs_regional_stats(ds=ds, infile = filename, mask=mask_frac, dims = dims, how = how,)
da_summary

The dataframe will be saved to: infile.replace(".nc", f"_summary-{'-'.join(how)}_ncra-regions.csv"

For example only, this would make a dataframe in this format:

Time Series extraction

For time series extraction of point locations see https://github.com/AusClimateService/TimeSeriesExtraction

Suggested regions, colormaps and scales for plotting

Using suggested colormaps and scales will improve the consistency across teams producing similar variables. This will support comparison across different plots.

We have provided dictionaries with suggested region shapefiles, cmap colormaps, and tick intervals. Using the recommended items may help make plotting between users more consistent, but if they are not fit for your purpose, you may specify whatever you like.

regions_dict

Region shape files are stored here: /g/data/ia39/aus-ref-clim-data-nci/shapefiles/data/

More information here: [https://aus-ref-clim-data-nci.github.io/aus-ref-clim-data-nci/datasets/shapefiles.html]

Regions include Australian:

• Local Government Areas (LGAs),
• State and Territories,
• land boundary,
• Natural Resource Management (NRM) regions,
• river regions,
• broadacre regions, and
• National Climate Risk Assessment (NCRA) regions.

dict_keys(['aus_local_gov', 'aus_states_territories', 'australia', 'nrm_regions', 'river_regions', 'broadacre_regions', 'ncra_regions'])

cmap_dict

These are suggested colormaps matched with possible variables to plot. This includes color maps for total amount and for anomalies.

tick_dict

# This dictionary gives some suggestions on the scale of the colour map to use for some variables. Some scales are taken from climate maps on bom.gov.au/climate
tick_dict = {"pr_annual":  [0, 50, 100, 200, 300, 400, 600, 1000, 1500, 2000, 3000, 6000],
             "pr_6mon":    [0, 50, 100, 200, 300, 400, 600,  900, 1200, 1800, 2400, 6000],
             "pr_3mon":  [0, 10,  25,  50, 100, 200, 300,  400,  600,  800, 1200, 2500],
             "pr_mon" :  [0,  1,   5,  10,  25,  50, 100,  200,  300,  400,  600, 1200],
             "pr_hour" :  [0, 1,   2,   5,  10,  15,  20,   30,   50,   75,  100,  200,],
             "pr_days": [0, 2, 3, 5, 10, 20, 30, 40, 50, 75, 100, 125, 150, 175],
             "pr_anom_mon": [-1000, -400, -200, -100, -50, -25, -10, 0, 10, 25, 50, 100, 200, 400, 1000],
             "pr_anom_3mon": [-2000, -600, -400, -200, -100, -50, -25, 0, 25, 50, 100, 200, 400, 600, 2000],
             "pr_anom_6mon": [-3000, -1200, -800, -400, -200, -100, -50, 0, 50, 100, 200, 400, 800, 1200, 3000],
             "pr_anom_ann": [-4000, -1800, -1200, -800, -400, -200, -100, 0, 100, 200, 400, 800, 1200, 1800, 4000],
             "pr_diff_mon": [-1000, -400, -200, -100, -50, -25, -10, 10, 25, 50, 100, 200, 400, 1000],
             "pr_diff_ann": [-3000, -1800, -1200, -800, -400, -200, -100, 100, 200, 400, 800, 1200, 1800, 3000],
             "frost_days": [0, 10, 20, 30, 40, 50, 75, 100, 150, 300],
             "frost_days_mon":[0, 2, 5, 10, 15, 20, 25, 31],
             "tas": np.arange(-9,52, 3),
             "tas_anom_day": np.arange(-14, 14.1, 2),
             "tas_anom_mon": np.arange(-7, 7.1, 1),
             "tas_anom_ann": np.arange(-3.5, 3.6, 0.5),
             "apparent_tas": np.arange(-6, 42, 3),
             "percent": np.arange(0,101,10),
             "xts_freq":[0.00, 0.005, 0.01, 0.02, 0.03, 0.05, 0.07, 0.10, 0.12, 0.15]
            }