Skip to content

Instantly share code, notes, and snippets.

@bennyistanto

bennyistanto/spei-terraclimate.md

Last active October 10, 2022 14:23
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save bennyistanto/efd8050f02d5e396900ad7f91b329f82 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save bennyistanto/efd8050f02d5e396900ad7f91b329f82 to your computer and use it in GitHub Desktop.
Download ZIP
Procedure to calculate TerraClimate-based Standardized Precipitation-Evapotranspiration Index with area of interest Indonesia
Raw
spei-terraclimate.md

TerraClimate-based Standardized Precipitation-Evapotranspiration Index in Indonesia

This section will explain on how to download TerraClimate's precipitation (ppt) and potential evapotranspiration (pet) monthly data in netCDF format and prepare it as input for Standardized Precipitation-Evapotranspiration Index (SPEI) calculation.

This step-by-step guide was tested using Mac mini Server - Late 2012, 2.3 GHz Quad-Core Intel Core i7, 16 GB 1600 MHz DDR3, running on macOS Catalina 10.15.7 and Windows 10 with Windows Subsystem for Linux enabled running on Thinkpad T480 2019, i7-8650U 1.9GHz, 64 GB 2400 MHz DDR4.

Index

0. Working Directory

For this tutorial, I am working on these folder /Users/bennyistanto/Temp/terraclimate/spei/ (applied to Mac/Linux machine) or Z:/Temp/terraclimate/spei/ (applied to Windows machine) directory. I have some folder inside this directory:

  1. downloads

    1. ppt
      1. nc_original
      2. nc_merge
      3. nc_tiles
      4. nc_subset
      5. nc_llt Place to put netCDF's ppt data that will use as an input
    2. pet
      1. nc_original
      2. nc_merge
      3. nc_tiles
      4. nc_subset
      5. nc_llt Place to put netCDF's pet data that will use as an input

    Place to put downloaded TerraCLimate data, and pre-process temporary files.

  2. outputs

    1. nc_original Output folder for SPEI calculation
      1. gamma
      2. pearson
    2. nc_tll Temporary for nc files from NCO arrange dimension process
    3. nc_merge Merging nc from separate tiles into single global layer (if you are following the global analysis procedures)
    4. geotiff
      1. gamma
        1. idn\spei12\ Final output of SPEI, generate by CDO and GDAL
      2. pearson

Feel free to use your own preferences for this setting/folder arrangements.

1. Software Requirement

If you encounter a problem, please look for a online solution. The installation and configuration described below is mostly performed using a bash shell on macOS. Windows users will need to install and configure a bash shell in order to follow the usage shown below. Try to use Windows Subsystem for Linux for this purpose.

1.1. macOS/Linux

Installing software for macOS/Linux

If you are new to using Bash refer to the following lessons with Software Carpentry: http://swcarpentry.github.io/shell-novice/

  • If you don't have Homebrew, you can install it by pasting below code in your macOS/Linux terminal.

     /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

  • Install wget (for downloading data). Use Hombrew to install it by pasting below code in your macOS terminal.

     brew install wget

  • Download and install Panoply Data Viewer from NASA GISS on your machine for macOS or Linux.

  • Download and install Anaconda Python on your machine for macOS or Linux.

1.2. Windows

Enable the Windows Subsytem for Linux

If you are using Windows machine, it's recomended to follow below step. You will experience an error during SPEI calculation cause by NCO if you use standard Windows 10 and not using Windows Subsytem for Linux.

Guideline below are specific for Windows 10. If you are using Windows Server 2019, please follow Windows Server Installation Guide

Reference: https://docs.microsoft.com/en-us/windows/wsl/install-win10

You must first enable the "Windows Subsystem for Linux - WSL" optional feature before installing any Linux distributions on Windows.

Installing software for Windows

If you have a Bash shell already installed on your Windows OS (e.g. Ubuntu Bash) you can use that for the exercise, but it must be a Bash shell

If you are new to using Bash refer to the following lessons with Software Carpentry: http://swcarpentry.github.io/shell-novice/

  • If you don't have Homebrew, you can install it by pasting below code in your WSL Ubuntu terminal.

     bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install.sh)"

  • Install wget (for downloading data). Use Hombrew to install it by pasting below code in your WSL Ubuntu terminal.

     brew install wget

  • Download and install Panoply Data Viewer from NASA GISS on your machine: Windows.

  • Download and install Anaconda Python on your WSL Ubuntu Linux: Ubuntu Linux on WSL.

    climate-indices python package used for SPEI calculation is rely on netCDF Operator (NCO) and pyNCO wrapper sometimes produce an error in Windows. That's the reason why we will use Anaconda for Linux if you are using Windows machine.

  • Go to https://repo.anaconda.com/archive/ to find the list of Anaconda releases

  • Select the release you want. I have a 64-bit computer, so I chose the latest release ending in x86_64.sh. If I had a 32-bit computer, I'd select the x86.sh version. If you accidentally try to install the wrong one, you'll get a warning in the terminal. I chose Anaconda3-2020.11-Linux-x86_64.sh.

  • From the terminal run wget https://repo.anaconda.com/archive/[YOUR VERSION]. Example:

     wget https://repo.anaconda.com/archive/Anaconda3-2020.11-Linux-x86_64.sh

  • After download process completed, Run the installation script: bash Anaconda[YOUR VERSION].sh

     bash Anaconda3-2020.11-Linux-x86_64.sh

  • Read the license agreement and follow the prompts to press Return/Enter to accept. Later will follow with question on accept the license terms, type yes and Enter. When asks you if you'd like the installer to prepend it to the path, press Return/Enter to confirm the location. Last question will be about initialize Anaconda3, type yes then Enter.

  • Close the terminal and reopen it to reload .bash configs. It will automatically activate base environment.

  • Deactivate base environment then set to false the confirguration of auto activate the base environment by typing

     conda deactivate && conda config --set auto_activate_base false

  • To test that it worked, which python in your Terminal. It should print a path that has anaconda in it. Mine is /home/bennyistanto/anaconda3/bin/python. If it doesn't have anaconda in the path, do the next step.

  • Manually add the Anaconda bin folder to your PATH. To do this, I added "export PATH=/home/bennyistanto/anaconda3/bin:$PATH" to the bottom of my ~/.bashrc file.

  • Optionally install Visual Studio Code when prompted

2. Configure the python environment

The code for calculating SPEI is written in Python 3. It is recommended to use either the Miniconda3 (minimal Anaconda) or Anaconda3 distribution. The below instructions will be Anaconda specific (although relevant to any Python virtual environment), and assume the use of a bash shell.

A new Anaconda environment can be created using the conda environment management system that comes packaged with Anaconda. In the following examples, I’ll use an environment named climate_indices (any environment name can be used instead of climate_indices) which will be created and populated with all required dependencies through the use of the provided setup.py file.

This step must only be done the first time. Once the environment has been created there is no need to do it again.

  • First, open your Terminal (in your macOS/Linux and Ubuntu Linux on WSL), create the Python environment with python3.7 as default:

     conda create -n climate_indices python=3.7

    Proceed with y

  • The environment created can now be ‘activated’:

     conda activate climate_indices

  • Install netCDF Operator (NCO), Climate Data Operator (CDO) from Max-Planck-Institut für Meteorologie and Geospatial Data Abstraction Library (GDAL/OGR) using conda and proceed with y.

     conda install -c conda-forge nco cdo gdal netcdf4

  • Install climate-indices package. Once the environment has been activated then subsequent Python commands will run in this environment where the package dependencies for this project are present. Now the package can be added to the environment along with all required modules (dependencies) via pip:

     pip install climate-indices

3. Preparing input

SPEI requires monthly precipitation and potential evapotranspiration data, for better result, SPEI required minimum 30-years of data.

If you are prefer to use your own dataset also fine, you can still follow this guideline and adjust some steps and code related to filename, unit, format and structure.

3.1. Input requirement

The climate-indices python package enables the user to calculate SPEI using any gridded netCDF dataset. However, there are certain requirements for input files that vary based on input type.

  • Precipitation and potential evapotranspiration unit must be written as millimeters, milimeter, mm, inches, inch or in.

  • Data dimension and order must be written as lat, lon, time (Windows machine required this order) or time, lat, lon (Works tested on Mac/Linux and Linux running on WSL).

  • If your study area are big, it's better to prepare all the data following this dimension order: lat, lon, time as all the data will force following this order during SPEI calculation by NCO module. Let say you only prepare the data as is (leaving the order to lat, lon, time), it also acceptable but it will required lot of memory to use re-ordering the dimension, and usually NCO couldn't handle all the process and failed.

3.2. Download TerraClimate data

  • There are 2 files contains link for downloading data_ppt.sh and data_pet.sh.

  • Download both files and put it in the same location with your working directory.

  • Navigate to downloads/terraclimate/ppt/nc_original and downloads/terraclimate/pet/nc_originalfolder in the working directory. Download using wget all terraclimate in netCDF format from Jan 1958 to Dec 2021 (this is lot of data, ppt +- 7.7GB and pet +- 6.4GB, please make sure you have bandwidth and unlimited data package). Paste and Enter below script in your Terminal.

    If you are in downloads/terraclimate/ppt/nc_original then execute below command

     sh data_ppt.sh

    If you are in downloads/terraclimate/pet/nc_original then execute below command

     sh data_pet.sh

3.3. Clip data using a bounding box based on area of interest and Merge netCDFs into single netCDF

  • Clip your area of interest using bounding box. We will use Indonesia (IDN) as the example case.

    Example: (IDN) bounding box with format lon1,lon2,lat1,lat2 is 94.75,141.25,-11.25,6.35

    Precipitation: Navigate your location to /downloads/ppt/nc_original

    Software Script
    CDO for fl in *.nc; cdo sellonlatbox,94.75,141.25,-11.25,6.35 $fl ../nc_subset/"idn_cli_"$fl; done
    NCO for fl in *.nc; ncks -d latitude,-11.25,6.35 -d longitude,94.75,141.25 $fl -O ../nc_subset/"idn_cli_"$fl; done

    Potential Evapotranspiration: Navigate your location to /downloads/pet/nc_original

    Software Script
    CDO for fl in *.nc; cdo sellonlatbox,94.75,141.25,-11.25,6.35 $fl ../nc_subset/"idn_cli_"$fl; done
    NCO for fl in *.nc; ncks -d latitude,-11.25,6.35 -d longitude,94.75,141.25 $fl -O ../nc_subset/"idn_cli_"$fl; done

  • Merge all annual netcdf in a folder into single netcdf.

    Precipitation: make sure you are in /downloads/ppt/nc_subset.

    Software Script
    CDO cdo mergetime idn_*.nc ../nc_merge/idn_cli_terraclimate_ppt_1958_2021.nc
    NCO Before merging the data, NCO required all the data has time dimension, we will use ncks command to make time the record dimension/variable used for concatenating files. Then use ncrcat to concatenate it.
    • for fl in *.nc; do ncks -O --mk_rec_dmn time $fl $fl; done
    • ncrcat -O -h idn_*.nc ../nc_merge/idn_cli_terraclimate_ppt_1958_2021.nc

    Potential Evapotranspiration: make sure you are in /downloads/pet/nc_subset.

    Software Script
    CDO cdo mergetime idn_*.nc ../nc_merge/idn_cli_terraclimate_pet_1958_2021.nc
    NCO Before merging the data, NCO required all the data has time dimension, we will use ncks command to make time the record dimension/variable used for concatenating files. Then use ncrcat to concatenate it.
    • for fl in *.nc; do ncks -O --mk_rec_dmn time $fl $fl; done
    • ncrcat -O -h idn_*.nc ../nc_merge/idn_cli_terraclimate_pet_1958_2021.nc

3.4. Check variable and attribute

As explain in Step 3.1. Input requirement above, we need to check the variable and attribute on above result to make sure all meet the requirements.

  • Navigate to /downloads/ppt/nc_merge folder in the working directory. Then execute below command.

     ncdump -h idn_cli_terraclimate_ppt_1958_2021.nc

  • Navigate to /downloads/pet/nc_merge folder in the working directory. Then execute below command.

     ncdump -h idn_cli_terraclimate_pet_1958_2021.nc

  • As you can see from above picture, all the requirement is completed: unit is in mm, order dimension for each variables is lat, lon, time, and time dimension is in UNLIMITED. Once this has completed, the dataset can be used as input to climate-indices package for computing SPEI.

4. Calculate the SPEI

Please make sure below points:

  • You are still inside climate_indices environment to start working on SPEI calculation.
  • Variable name on precipitation --var_name_precip, usually terraclimate data use ppt as name while other precipitation data like CHIRPS using precip and IMERG using precipitation as a variable name. To make sure, check using command ncdump -h file.nc then adjust it in SPEI script if needed.
  • Variable name on potential evapotranspiration --var_name_pet, usually terraclimate data use pet as name.
  • Precipitation and potential evapotranspiration unit must be written as millimeters, milimeter, mm, inches, inch or in.
  • Data dimension and order must be written as lat, lon, time (Windows machine required this order) or time, lat, lon (Works tested on Mac/Linux and Linux running on WSL).

Let's start the calculation!

  • In your Terminal, run the following code.

     process_climate_indices --index spei --periodicity monthly --netcdf_precip /Users/benny/Temp/TERRACLIMATE/spei/downloads/ppt/nc_merge/idn_cli_terraclimate_ppt_1958_2021.nc --var_name_precip ppt --netcdf_pet /Users/benny/Temp/TERRACLIMATE/spei/downloads/pet/nc_merge/idn_cli_terraclimate_pet_1958_2021.nc --var_name_pet pet --output_file_base /Users/benny/Temp/TERRACLIMATE/spei/outputs/nc_original/idn_cli --scales 1 2 3 6 9 12 18 24 36 48 60 72 --calibration_start_year 1958 --calibration_end_year 2021 --multiprocessing all

  • Above code is example for calculating SPEI 1 to 72-months. It's ok if you think you only need some of them. Example: you are interested to calculate SPEI 1 - 3-months or SPEI 12-months, then adjust above code into --scales 1 2 3 or --scales 12.

  • The above command will compute SPEI (both gamma and Pearson Type III distributions) from monthly precipitation dataset and potential evapotranspiration, and the calibration period used will be Jan-1958 through Dec-2021. The index will be computed at 1,2,3,6,9,12,18,24,36,48,60 and 72-month timescales. The output files will be <out_dir>/idn_cli_spei_gamma_xx.nc, and <out_dir>/idn_cli_spei_pearson_xx.nc.

    The output files will be:

    Gamma

    1. 1-month: /outputs/nc_original/idn_cli_spei_gamma_01.nc
    2. 2-month: /outputs/nc_original/idn_cli_spei_gamma_02.nc
    3. 3-month: /outputs/nc_original/idn_cli_spei_gamma_03.nc
    4. 6-month: /outputs/nc_original/idn_cli_spei_gamma_06.nc
    5. 9-month: /outputs/nc_original/idn_cli_spei_gamma_09.nc
    6. 12-month: /outputs/nc_original/idn_cli_spei_gamma_12.nc
    7. 18-month: /outputs/nc_original/idn_cli_spei_gamma_18.nc
    8. 24-month: /outputs/nc_original/idn_cli_spei_gamma_24.nc
    9. 36-month: /outputs/nc_original/idn_cli_spei_gamma_36.nc
    10. 48-month: /outputs/nc_original/idn_cli_spei_gamma_48.nc
    11. 60-month: /outputs/nc_original/idn_cli_spei_gamma_60.nc
    12. 72-month: /outputs/nc_original/idn_cli_spei_gamma_72.nc

    Pearson

    1. 1-month: /outputs/nc_original/idn_cli_spei_pearson_01.nc
    2. 2-month: /outputs/nc_original/idn_cli_spei_pearson_02.nc
    3. 3-month: /outputs/nc_original/idn_cli_spei_pearson_03.nc
    4. 6-month: /outputs/nc_original/idn_cli_spei_pearson_06.nc
    5. 9-month: /outputs/nc_original/idn_cli_spei_pearson_09.nc
    6. 12-month: /outputs/nc_original/idn_cli_spei_pearson_12.nc
    7. 18-month: /outputs/nc_original/idn_cli_spei_pearson_18.nc
    8. 24-month: /outputs/nc_original/idn_cli_spei_pearson_24.nc
    9. 36-month: /outputs/nc_original/idn_cli_spei_pearson_36.nc
    10. 48-month: /outputs/nc_original/idn_cli_spei_pearson_48.nc
    11. 60-month: /outputs/nc_original/idn_cli_spei_pearson_60.nc
    12. 72-month: /outputs/nc_original/idn_cli_spei_pearson_72.nc

Parallelization will occur utilizing all CPUs.

When the SPEI calculation completed, move arrange all the output by moving SPEI files with gamma to gammafolder and with pearson to pearson folder.

For the translation to GeoTIFF as a final output, we only use SPEI gamma version.

5. Visualize the result using Panoply

Let see the result.

  • From the /outputs/nc_original/gamma directory, right-click file idn_cli_spei_gamma_12.nc and Open With Panoply.

  • From the Datasets tab select spei_gamma_12_month and click Create Plot

  • In the Create Plot window select option Georeferenced Longitude-Latitude.

  • When the Plot window opens:

    • Array tab: Change the time into 744 to view data on 1 December 2020
    • Scale tab: Change value on Min -3.09, Max 3.09, Major 6, Color Table CB_RdBu_09.cpt
    • Map tab: Change value on Center on Lon 120.0 Lat -2.0, then Zoom in the map through menu-editor Plot > Zoom - Plot In few times until Sao Tome and Principe appear proportionally. Set grid spacing 1.0 and Labels on every grid lines.
    • Overlays tab: Change Overlay 1 to MWDB_Coasts_Countries_1.cnob

6. Convert the result to GeoTIFF

We need CDO to do a conversion of the result into GeoTIFF format, and CDO required the variable should be in time,lat,lon, while the output from SPEI: idn_cli_spei_gamma_x_month.nc in lat,lon,time, you can check this via ncdump -h idn_cli_spei_gamma_12.nc

  • Navigate your Terminal to folder /outputs/nc_original/gamma/

  • Let's re-order the variables into time,lat,lon using ncpdq command from NCO and save the result to folder /outputs/nc_tll/idn/

     ncpdq -a time,lat,lon idn_cli_spei_gamma_12.nc ../../../outputs/nc_tll/idn/tp_cli_spei_gamma_12.nc

  • Navigate your Terminal to folder /outputs/nc_tll/idn/

  • Check result and metadata to make sure everything is correct.

     ncdump -h idn_cli_spei_gamma_12.nc

  • Then convert all idn_cli_spei_gamma_12.nc value into GeoTIFF with time dimension information as the filename using CDO and GDAL. Usually the geotiff image will not have projection information, so we will add that information via the script: -a_ullr ulx uly lrx lry -a_srs EPSG:4326

  • Execute below script and save the result to folder /outputs/geotiff/gamma/idn/spei12

     for t in `cdo showdate idn_cli_spei_gamma_12.nc`; do
     cdo seldate,$t idn_cli_spei_gamma_12.nc dummy.nc     
     gdal_translate -of GTiff -a_ullr 94.75 6.35 141.25 -11.25 -a_srs EPSG:4326 -co COMPRESS=LZW -co PREDICTOR=1 dummy.nc ../../geotiff/gamma/idn/spei12/idn_cli_terraclimate_spei12_$t.tif
 done

  • Next, you can continue to translate other SPEI files.

Congrats, now you are able to calculate SPEI based on monthly rainfall in netCDF and translate the output into GeoTIFF format.

7. Indonesia SPEI data

Indonesia SPEI-12 data using terraclimate in GeoTIFF format has been extracted, other SPEI for different time period will follow. Now please check or open SPEI-12 data for December 2020 using QGIS software. If you don't have the software, you can download it for free via this link.

8. Comparison with other data provider

SPEI Global Drought Monitor

Santiago Begueria and friend from Spanish National Research Council who invented SPEI, released the SPEI Global Drought Monitor which offers near real-time information about drought conditions at the global scale, with a 1 degree spatial resolution and a monthly time resolution.

Link for SPEI 12-month, December 2020 from Global Drought Monitor - https://spei.csic.es/map/maps.html#months=4#month=11#year=2020

Climate Engine

Climate Engine is a free web application powered by Google Earth Engine that can be used to create on-demand maps and charts from publicly available satellite and climate data using a standard web browser. Climate Engine allows users to analyze and interact with climate and earth observations for decision support related to drought, water use, agricultural, wildfire, and ecology.

One of the product that could generate easily using Climate Engine is SPEI and using terraclimate data. Link https://climengine.page.link/Aym1 - Register to Climate Engine platform may required to access the link

9. Reference

  1. https://climatedataguide.ucar.edu/climate-data/standardized-precipitation-evapotranspiration-index-spei
  2. https://spei.csic.es
  3. http://www.climatologylab.org/terraclimate.html
  4. https://pypi.org/project/climate-indices/
  5. https://climate-indices.readthedocs.io/en/latest/
  6. https://code.mpimet.mpg.de/projects/cdo
  7. http://nco.sourceforge.net
  8. https://appliedsciences.nasa.gov/join-mission/training/english/arset-applications-gpm-imerg-reanalysis-assessing-extreme-dry-and-wet
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment