Worksheet 3: Thresholds and extremes analysis


The following exercises demonstrate analysis of moderate extremes in climate simulated using PRECIS. As with the other worksheets, these are just examples of some of the analysis that you might perform using CDO.



Contents


3.1 Frequency of wet days
3.2 Calculating percentiles
3.3 Calculating extremes Indices


Note: In all of the following exercises, it is assumed that the spin up period and the 8-point rim have been removed.



3.1 Frequency of wet days


CDO can calculate many different kinds of thresholds and indicies, here we begin by finding the frequency of wet days using daily data.

1 a. ) For both cahpa and cahpb calculate the percentage of days in both the baseline and the future period which are wet days (define a wet day as having precipitation >=1 mm/day.)

gec,c this command calculates where a field is greater than or equal to a constant c. It returns a mask of ones where the criteria is met and zeroes where it is not, similar operators exist for less than, equal to etc.
Note: Many of the threshold and index commands may take a few minutes to complete due to the amount of data being processed.

% DATA=/media/sf_kl2/
% cd $DATA/daily
% ls
% mkdir climatology
% for runid in "cahpa" "cahpb"
% do
% cd $DATA/daily/$runid/05216
% cdo timsum -gec,1 -mulc,86400 ${runid}a.pa.6190.05216.rr8.nc $DATA/daily/climatology/${runid}a.wetday.no.baseline.nc
% cdo timsum -gec,1 -mulc,86400 ${runid}a.pa.2150.05216.rr8.nc $DATA/daily/climatology/${runid}a.wetday.no.future.nc
% done

Now as a percentage, there are 10800 days in each period for cahpa which uses a 360 day calendar and 10957 for cahpb which uses a real (365.25 day) calendar.

% cd $DATA/daily/climatology
% cdo mulc,100 -divc,10800 cahpaa.wetday.no.baseline.nc cahpaa.wetday.per.baseline.nc
% cdo mulc,100 -divc,10800 cahpaa.wetday.no.future.nc cahpaa.wetday.per.future.nc
% cdo mulc,100 -divc,10957 cahpba.wetday.no.baseline.nc cahpba.wetday.per.baseline.nc
% cdo mulc,100 -divc,10957 cahpba.wetday.no.future.nc cahpba.wetday.per.future.nc
% xconv -i cahpaa.wetday.per.future.nc

1 b. ) For 1961-1990 APHRODITE observations calculate the percentage of days in the baseline which are wet days.

% cd $DATA/daily/APHRODITE
% cdo timsum -gec,1 aphro.day.6190.nc $DATA/daily/climatology/aphro.wetday.no.nc
% cd $DATA/daily/climatology
% cdo mulc,100 -divc,10957 aphro.wetday.no.nc aphro.wetday.per.nc
% xconv -i aphro.wetday.per.nc

1 c. ) Calculate the difference in future and baseline wet day frequency and also baseline and observation wet day frequency.

% cdo sub cahpba.wetday.per.future.nc cahpba.wetday.per.baseline.nc cahpba.wetday.per.diff.nc
% cdo sub cahpaa.wetday.per.future.nc cahpaa.wetday.per.baseline.nc cahpaa.wetday.per.diff.nc

Remember, we must re-grid model data onto the observation grid (as done in worksheet 2)

% cdo griddes aphro.wetday.per.nc > mygrid
% cdo remapbil,mygrid cahpba.wetday.per.baseline.nc cahpba.wetday.per.baseline.rg.nc
% cdo sub cahpba.wetday.per.baseline.rg.nc aphro.wetday.per.nc cahpba.wetday.per.obs_diff.nc
% cdo remapbil,mygrid cahpaa.wetday.per.baseline.nc cahpaa.wetday.per.baseline.rg.nc
% cdo sub cahpaa.wetday.per.baseline.rg.nc aphro.wetday.per.nc cahpaa.wetday.per.obs_diff.nc

2 a. ) Plot the difference in the percentage of wet day frequency between the models and observations

The script wetday_diff.ncl uses the fields generated in section 3.1. It reads in the NetCDF fields and outputs a .png plot called wetday_diff.png. Open the script for editing.

% cd $DATA/ncl_scripts
% emacs wetday_diff.ncl & # you can use which ever text editor you prefer

Before running this script two changes need to be made. At the top of the script, change data_dir to be the location of the climatology directory which contains the data to be plotted. Also change plot_location to the location of the ncl_scripts directory.

% ncl wetday_diff.ncl
% xv wetday_diff.png

2 b. ) Improve the plot by adding a title and changing the colour map

To make these changes, the script wetday_diff.ncl must be edited.

% emacs wetday_diff.ncl &

To add a main title go to the create panel section at the bottom of the script and change the following resource from ---TITLE--- to an appropriate title:

- resP@txString main title

The colour map currently in use is not suitable. In the setup plot window and colour map section of the script change cols from 'example' to one of the following colour maps, choose which one you think best suits the data and re-run the script.

- GreenYellow
- rainbow
- BlueRed
- WhiteBlueGreenYellowRed

% ncl wetday_diff.ncl

3.2 Calculating percentiles

1 a. ) Calculate in mm/day the baseline (1961-1990) and future (2021-2050) 95th percentile of precipitation. Do this for cahpa, cahpb and also for APHRODITE baseline.

timpctl,c command calculates the cth percentile of each point in a field over all time steps. To use this command you need to have calculated the minimum and maximum of the field over all time steps e.g. cdo timpctl,c [ifile] [min_file] [max_file] [ofile]
timmax/timin commands find the maximum/minimum of a field over all time steps.
This loop will take a few minutes to complete.

% for runid in "cahpa" "cahpb"
% do
% cd $DATA/daily/$runid/05216
% infile=${runid}a.pa.6190.05216.rr8.mmday.nc
% cdo timpctl,95 $infile -timmin $infile -timmax $infile $DATA/daily/climatology/${runid}a.pc95.05216.baseline.mmday.nc
% infile=${runid}a.pa.2150.05216.rr8.mmday.nc
% cdo timpctl,95 $infile -timmin $infile -timmax $infile $DATA/daily/climatology/${runid}a.pc95.05216.future.mmday.nc
% done

% cd $DATA/daily/APHRODITE
% infile=aphro.day.6190.nc
% cdo timpctl,95 $infile -timmin $infile -timmax $infile $DATA/daily/climatology/aphro.pc95.05216.baseline.nc

1 b. ) What is the difference between the future and baseline, and baseline and APHRODITE observations 95th percentile of precipitation?

% cd $DATA/daily/climatology
% cdo sub cahpaa.pc95.05216.future.mmday.nc cahpaa.pc95.05216.baseline.mmday.nc cahpaa.pc95.05216.diff.mmday.nc
% cdo sub cahpba.pc95.05216.future.mmday.nc cahpba.pc95.05216.baseline.mmday.nc cahpba.pc95.05216.diff.mmday.nc

% cdo griddes aphro.pc95.05216.baseline.nc > mygrid
% cdo remapbil,mygrid cahpaa.pc95.05216.baseline.mmday.nc cahpaa.pc95.05216.baseline.mmday.rg.nc
% cdo sub cahpaa.pc95.05216.baseline.mmday.rg.nc aphro.pc95.05216.baseline.nc cahpaa.95pc.05216.obs_diff.mmday.nc
% cdo remapbil,mygrid cahpba.pc95.05216.baseline.mmday.nc cahpba.pc95.05216.baseline.mmday.rg.nc
% cdo sub cahpba.pc95.05216.baseline.mmday.rg.nc aphro.pc95.05216.baseline.nc cahpba.95pc.05216.obs_diff.mmday.nc

2. ) Calculate the baseline (1961-1990) and future (2021-2050) 90th percentile of maximum temperature and the difference between them. Do this for cahpa and cahpb.

To do this edit and run the script worksheet4_percentiles.sh which is currently set up for the 95th percentile of precipitation.

% cd $DATA
% emacs worksheet4_percentiles.sh

and after editing

% ./worksheet4_percentiles.sh

How does the 90th percentile of maximum temperature change in the future?




3 a. ) Plot the future change in the 95th percentile of precipitation for both models

The script pc95_precip_diff.ncl uses the fields generated in section 3.2 of the difference between the future and baseline 95th percentile. It reads in the NetCDF fields and outputs a .png plot called pc95_precip_diff.png. Open the script for editing.

% cd $DATA/ncl_scripts
% emacs pc95_precip_diff.ncl &

Before running this script two changes need to be made. At the top of the script, change data_dir to be the location of the climatology directory which contains the data to be plotted. Also change plot_location to the location of the ncl_scripts directory.

% ncl pc95_precip_diff.ncl
% xv pc95_precip_diff.png

3 b. ) Improve the plot by changing the contour levels and number of panels to reflect the data.

To make these changes, the script pc95_precip_diff.ncl must be edited.

% emacs pc95_precip_diff.ncl &

To change the contour levels go to the Set up model plots section in the middle of the script and change the following resources to fit the values of the data better:

- res@cnMinLevelValF set min contour level
- res@cnMaxLevelValF set max contour level
- res@cnLevelSpacingF contour spacing

The plot window is currently set up for four plots, to change this to 2 plots, we must make two changes. Firstly, in the setup plot window and colour map section of the script change plot from new(4,graphic) to new(2,graphic) i.e. the number of plots you need.

Secondly at the end of the script in the section create panel change gsn_panel(wks,plot,(/1,4/),resP) to either (/1,2/) or (/2,1/), where (/no. of rows, no. of colums/). Then re-run the script.

% ncl pc95_precip_diff.ncl

3.3 Calculating Extremes Indices

CDO can be used to calculate all kinds of extremes indices, here we calculate a few moderate extremes.

1 . ) Calculate the frequency of warm days in the future (extreme index TX90P.)

To do this we need to use the 90th percentile of baseline maximum temperature which was calculated in 3.2. In the baseline this threshold should be exceeded 10% of the time.
gt command produces a mask of zeroes and ones where ifile1 is greater than ifile2. Similar operators exist for equal to, less than etc.
% cd $DATA/daily/cahpa/03236.max
% cdo timsum -gt -subc,273.15 cahpaa.pa.2150.03236.max.rr8.nc $DATA/daily/climatology/cahpaa.pc90.03236.max.baseline.degC.nc $DATA/daily/climatology/cahpaa.TX90P.nc

% cd $DATA/daily/cahpb/03236.max
% cdo timsum -gt -subc,273.15 cahpba.pa.2150.03236.max.rr8.nc $DATA/daily/climatology/cahpba.pc90.03236.max.baseline.degC.nc $DATA/daily/climatology/cahpba.TX90P.nc

Now calculate as a percentage.

% cd $DATA/daily/climatology
% cdo mulc,100 -divc,10800 cahpaa.TX90P.nc cahpaa.TX90P.per.nc
% cdo mulc,100 -divc,10957 cahpba.TX90P.nc cahpba.TX90P.per.nc

2 a. ) Calculate the percentage of total precipitation which falls on very wet days in the future (where a very wet day is one on which daily rainfall exceeds the 95th percentile of the baseline) over Malaysia.

Note: This loop will take a few mintues to complete, it extracts the Malaysia area from the full field.

% for runid in "cahpa" "cahpb"
% do
% cd $DATA/daily/$runid/05216
% cdo gt -sellonlatbox,99,122,0,9 ${runid}a.pa.2150.05216.rr8.mmday.nc -sellonlatbox,99,122,0,9 $DATA/daily/climatology/${runid}a.pc95.05216.baseline.mmday.nc where_wet.nc
% cdo timsum -mul where_wet.nc -sellonlatbox,99,122,0,9 ${runid}a.pa.2150.05216.rr8.mmday.nc wettot.nc
% cdo mulc,100 -div wettot.nc -timsum -sellonlatbox,99,122,0,9 ${runid}a.pa.2150.05216.rr8.mmday.nc $DATA/daily/climatology/${runid}a.R95pTOT.future.nc
% rm where_wet.nc
% done

3. ) Plot the frequency of warm days in the future and the future percentage of total precipitation falling on very wet days

Edit the script change_precip_temp.ncl used in worksheet 2b to plot the TX90P.per and R95pTOT.future files for cahpa and cahpb files generated in 3.3.

% cd $DATA/ncl_scripts
% cp change_precip_temp.ncl extremes.ncl
% emacs extremes.ncl &

The file names in the read in data precip and read in data temp sections will need to be changed to the files to be plotted, and the name of the plot (plotname) at the top of the script will need to be changed to extremes.

The data will be plotted upside down, add the following resource to the setup plot window and colour map section of the script for both R95pTOT and TX90P:

- res@trYReverse = True reverse y axis

The contour levels will need to be changed to reflect the data, as will the titles. Change the following resources:

- res@cnMinLevelValF set min contour level
- res@cnMaxLevelValF set max contour level
- res@cnLevelSpacingF contour spacing
- resP@txString main title
- res@gsnRightString right title
- res@gsnLeftString left title

You may need to change the color bar, to do this you can change the named colours in precip_cols and temp_cols in the setup plot window and colour map section of the script.

% ncl extremes.ncl
% display extremes.png

You will probably need to edit and run the script multiple times to improve the resulting plot.

How much of total precipitation is falling on very wet days in the future? Is there an increase or decrease in the number of hot days?