How to download, mosaic and resample Modis NDVI data

Modis NDVI data are fully described here

Briefly: these data are intended to be the continuation of the well-known NOAA-AVHRR Vegetation Index and are currently used in a wide range of applications.

The following procedure (very basic!!! use at your own risk! no waranty!!) works only in a Linux environment (even if it has been ported on a Windows machine under Cygwin).
Pre-requisities: wget and MRT (Modis Reprojection Tool) that needs to be installed properly (that is MRT's bin directory needs to be added to the environmental variable $PATH)
It consists in three Bash script and one auxiliary file, namely:
  • wrapper_ModisNDVI.sh
  • get.ndvi.sh
  • mrt.ndvi.sh
  • custom.env
Copy all the files in your workin directory and cd to that directory.

Usage: sh wrapper_ModisNDVI.sh YYYY.MM.DD file_env
where YYYY.MM.DD is the year.month.day of the NDVI time step you're interested in as listed in https://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.006/ and the file file_env (custom.env) contains several settings to run correctly the procedure.
The outputs of the procedure are a NDVi tif layer, a VI_Quality tif layer and a pixel_reliability layer.
The file file_env (custom.env) should be self explanatory, if not write me and I'll try to detail the meaning of the variables.

Example custom.env:
 # Prerequisities: wget, MRT
 # Tiles settings  
 PRODUCT_ID='MOD13Q1'  
 PRODUCT_VERSION='006'  
 htiles='17 18 19 20'  # horizontal tiles you want to download and mosaic
 vtiles='02 03 04 05'  # vertical tiles you want to download and mosaic
 # Remote settings  
 http_remote_server='http://e4ftl01.cr.usgs.gov/MOLT/'${PRODUCT_ID}'.'${PRODUCT_VERSION}  
 remote_user='your_user_name'  
 remote_pass='your_password'  
 # Local settings  
 LOGDIR=$HOME'/log'  # this is where you want to store the log files
 RAW_DATA_DIR=$HOME'/data'  # this is where you want to store the data
 MAXTRY=100  # a download setting
 SLEEP_TIME=60  # a download setting
 # MRT settings  
 MRT_BIN_DIR=$HOME'/build/MRT/installdir/bin'  # where you installed the MRT software
 lonmin=-10.0  # need to change for your custom ROI
 lonmax=30.0  # need to change for your custom ROI
 latmin=30.0  # need to change for your custom ROI
 latmax=65.0  # need to change for your custom ROI
 pixel_size='0.0020833'  # approximately 250 m
 #flag_mosaic='TRUE'  
 flag_del='TRUE'  # delete the hdf files
 flag_proj='TRUE'  
 datum='WGS84'  
 bands_subset='1 0 1 0 0 0 0 0 0 0 0 1'  
 proj_type='GEO'  
 res_type='NEAREST_NEIGHBOR'  
 spatial_subset='INPUT_LAT_LONG'

The wrapper Bash file is:

 #!bin/bash  
 #############################################################################  
 #  
 # Project : Wrapper for NDVI MODIS Terra  
 #      Download, mosaic and resample MOD13Q1 data product  
 # Author : valcap74@gmail.com  
 # Date  : 20171024  
 # Version : 1.0  
 #  
 #############################################################################  
 if [ $# -ne 2 ]; then  
  echo "Not enough arguments."  
  echo  
  echo "Usage:"  
  echo "     $0 YYYY.MM.DD file_env"  
  echo "Example: $0 2000.04.06 ./eurodeer.env"  
  exit 1  
 fi  
 source /etc/profile  
 # Enviromnental variables  
 file_env=$2  
 if [ ! -e $file_env ]; then  
  echo "$file_env is missing"  
  exit 1  
 else  
  source $2  
 fi  
 # Log  
 #exec 1>$LOGDIR/${0}_${PRODUCT_ID}_${1}.log  
 #exec 2>&1  
 # Check whether MRT is installed  
 if [ -e $MRT_BIN_DIR/resample ]; then  
  echo resample OK  
 else  
  echo ops resample is missing and is mandatory  
 fi  
 if [ -e $MRT_BIN_DIR/mrtmosaic ]; then  
  echo mrtmosaic OK  
 else  
  echo ops mrtmosaic is missing and is mandatory  
 fi  
 res=`which wget`  
 if [ $? -ne 0 ]; then  
  echo 'ops wget is missing and is mandatory'  
 fi  
 # Loop over procedures  
 for proc in get.ndvi.sh mrt.ndvi.sh  
 do  
  if [ ! -e $proc ]; then  
   echo "ops $proc is missing and mandatory..."  
   exit 1  
  fi  
  echo `date +%H:%M:%S`" - Running $proc $1 $2"  
  sh ./$proc $1 $2  
  if [ $? -eq 0 ]; then  
   echo `date +%H:%M:%S`" - $proc $1 $2 ended successfully"   
  else  
   echo "ops problem in $proc $1 $2"  
   exit 1  
  fi  
 done  
 exit 0;

The downloading Bash file is:

 #!/bin/bash  
 #############################################################################  
 #  
 # Project : Get MOD13Q1 data from http://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.005  
 # Author : valcap74@gmail.com  
 # Date  : 20171002  
 # Version : 1.0  
 #  
 #############################################################################  
 if [ $# -ne 2 ]; then  
  echo "Not enough arguments."  
  echo  
  echo "Usage:"  
  echo "     $0 YYYY.MM.DD file_env"  
  echo "Example: $0 2000.04.06 ./eurodeer.env"  
  exit 1  
 fi  
 source /etc/profile  
 # Enviromnental variables  
 file_env=$2  
 if [ ! -e $file_env ]; then  
  echo "$file_env is missing"  
  exit 1  
 else  
  source $2  
 fi  
 # Parse input date  
 YYYY=`echo $1 | cut -c1-4`  
 MM=`echo $1 | cut -c6-7`  
 DD=`echo $1 | cut -c9-10`  
 # Log  
 rm -f $LOGDIR/${0}_${PRODUCT_ID}_${1}.log  
 exec 1>$LOGDIR/${0}_${PRODUCT_ID}_${1}.log  
 exec 2>&1  
 # check if file exists  
 cd $RAW_DATA_DIR  
 for htile in $htiles; do  
  for vtile in $vtiles; do  
   echo "Checking if $1 is available at $http_remote_server"  
   wget -O - $http_remote_server | grep $1 >&/dev/null  
   if [ $? -eq 0 ]; then  
    echo "OK!"  
    day=`date -d "${YYYY}/${MM}/${DD}" +%j`  
    wget -O - $http_remote_server/${1} | grep A${YYYY}${jday} | grep h${htile} | grep v${vtile} | grep hdf | grep -v xml > contains_file_name  
    filename=`cat contains_file_name | cut -d"=" -f4 | cut -d '>' -f1 | sed -e "s/\"//g"`  
    if [ ! -z "$filename" ]; then  
      rm -f contains_file_name  
      echo `date +%H:%M:%S`" - Downloading $filename from $http_remote_server/${1}"  
      ntry=1  
      for (( itry=$ntry; itry<=$MAXTRY; itry++ ));do  
       echo `date +%H:%M:%S`" - Try num: "$itry  
       wget --user $remote_user --password $remote_pass $http_remote_server/${1}/$filename >&/dev/null  
       if [ $? -eq 0 ];then  
       echo `date +%H:%M:%S`" - $filename OK!"  
        break  
      fi  
      sleep $SLEEP_TIME  
      done  
      echo "NO! Failed to download $filename"  
    else  
     echo "NO! $1 is present at $http_remote_server but it does not contain any h${htile} and v${vtile} tiles or it is not present the hdf file format"  
    fi  
   else  
    echo "NO! $1 is not available at $http_remote_server"  
   fi  
  done  
 done  
 exit 0;

The MRT Bash file is:

 #!bin/bash  
 #############################################################################  
 #  
 # Project : Process MOD13Q1 data using MRT  
 # Author : valcap74@gmail.com  
 # Date  : 20171024  
 # Version : 1.0  
 #  
 #############################################################################  
 if [ $# -ne 2 ]; then  
  echo "Not enough arguments."  
  echo  
  echo "Usage:"  
  echo "     $0 YYYY.MM.DD file_env"  
  echo "Example: $0 2000.04.06 ./eurodeer.env"  
  exit 1  
 fi  
 source /etc/profile  
 # Enviromnental variables  
 file_env=$2  
 if [ ! -e $file_env ]; then  
  echo "$file_env is missing"  
  exit 1  
 else  
  source $2  
 fi  
 # Parse input date  
 YYYY=`echo $1 | cut -c1-4`  
 MM=`echo $1 | cut -c6-7`  
 DD=`echo $1 | cut -c9-10`  
 # Log  
 exec 1>$LOGDIR/${0}_${PRODUCT_ID}_${1}.log  
 exec 2>&1  
 cd $RAW_DATA_DIR  
 export MRTDATADIR=$RAW_DATA_DIR  
 ############################  
 # First mosaic all the tiles  
 echo "Running mrtmosaic"  
 rm -f ./TmpMosaic.prm Mosaic_${YYYY}-${MM}-${DD}.hdf  
 jday=`date -d "${YYYY}/${MM}/${DD}" +%j`  
 for htile in $htiles; do  
  for vtile in $vtiles; do  
  res_file=`ls ${PRODUCT_ID}.A${YYYY}${jday}.h${htile}v${vtile}.${PRODUCT_VERSION}.*.hdf`  
  if [ $? -eq 0 ]; then  
   cat << EOF >> ./TmpMosaic.prm  
 $res_file  
 EOF  
  else  
   echo "ops file $htile $vtile ${YYYY}${jday} is missing and it is mandatory"  
   exit 1;  
  fi  
  done  
 done  
 res=`${MRT_BIN_DIR}/mrtmosaic -i ./TmpMosaic.prm -s "$bands_subset" -o Mosaic_${YYYY}-${MM}-${DD}.hdf`  
 if [ $? -eq 0 ]; then  
  echo OK!  
  rm -f ./TmpMosaic.prm  
 else  
  echo "ops problem in mrtmosaic"  
  echo "output is $res"  
 fi  
 ############################  
 # Second resample the output tile  
 echo "Running resample..."  
 cat << EOF > ./resample.prm  
 INPUT_FILENAME = $RAW_DATA_DIR/Mosaic_${YYYY}-${MM}-${DD}.hdf  
 SPECTRAL_SUBSET = ( 1 1 1 )  
 SPATIAL_SUBSET_TYPE = $spatial_subset  
 SPATIAL_SUBSET_UL_CORNER = ( $latmax $lonmin )  
 SPATIAL_SUBSET_LR_CORNER = ( $latmin $lonmax )  
 OUTPUT_FILENAME = $RAW_DATA_DIR/${PRODUCT_ID}_${YYYY}_${MM}_${DD}.tif  
 RESAMPLING_TYPE = $res_type  
 OUTPUT_PROJECTION_TYPE = $proj_type   
 OUTPUT_PROJECTION_PARAMETERS = (   
 0.0 0.0 0.0  
 0.0 0.0 0.0  
 0.0 0.0 0.0  
 0.0 0.0 0.0 )  
 DATUM = $datum  
 OUTPUT_PIXEL_SIZE = $pixel_size  
 EOF  
 res=`${MRT_BIN_DIR}/resample -p ./resample.prm`  
 if [ $? -eq 0 ]; then  
  echo OK!  
  rm -f ./resample.prm  
 else  
  echo "ops problem in resample"  
  echo "output is $res"  
 fi  
 # del .hdf files  
 if [ "$flag_del" == "TRUE" ]; then  
  rm -fv $RAW_DATA_DIR/Mosaic_${YYYY}-${MM}-${DD}.hdf  
  jday=`date -d "${YYYY}/${MM}/${DD}" +%j`  
  for htile in $htiles; do  
   for vtile in $vtiles; do  
   res_file=`ls $RAW_DATA_DIR/${PRODUCT_ID}.A${YYYY}${jday}.h${htile}v${vtile}.${PRODUCT_VERSION}.*.hdf`  
   if [ $? -eq 0 ]; then  
    rm -fv $RAW_DATA_DIR/$res_file  
   fi  
   done  
  done  
 fi  
 exit 0;

No quality check of the NDVI pixel values is performed (for additional post-processing you have the VI_Quality layer).

A sample image of the final output is below:




Commenti

  1. valcap26 ottobre 2017 alle ore 15:04

    Bug fixes: there's an error in the downloading procedure...
    line "filename=`cat contains_file_name | cut -d"=" -f4 | cut -d '>' -f1 | sed -e "s/\"//g"`" needs editing to parse properly the html string

    RispondiElimina

Posta un commento

Post popolari in questo blog

How to run the WRF model using ERA5 (on model levels) as initial and boundary conditions

This procedure was tested with the WRF model Version 3.8.* and ERA5 data. 1. ERA5 data were downloaded as described in How+to+download+ERA5+data+via+the+ECMWF+Web+API   Here I'm attaching a sample script (launched on ecgate) to download surface data: # General settings Nort='90.0' West='-90.0' Sout='0.0' East='90.0' resX='0.25' resY='0.25' initdate='2012-11-01' paralist='msl/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4' timelist='00:00:00/01:00:00/02:00:00/03:00:00/04:00:00/05:00:00/06:00:00/07:00:00/08:00:00/09:00:00/10:00:00/11:00:00/12:00:00/13:00:00/14:00:00/15:00:00/16:00:00/17:00:00/18:00:00/19:00:00/20:00:00/21:00:00/22 :00:00/23:00:00' gribfile='ERA5_sfc_'${initdate} mars <<EOF retrieve, class=ea, dataset=era5, date=${initdate}, expver=1, levtype=sfc, param=${paralist}, stream=oper, time=${t
Continua a leggere

Yet another tutorial on how to install WRF

Immagine
Several tutorials on how to install WRF are available online. Basically the common root is the Official Tutorial from mmm.ucar . Some other noticeable websites are: WRF on Ubuntu server 16.04 again WRF on Ubuntu server 16.04 WRF on CENTOs WRF on Archlinux (with some nice ncl examples) WRF on a generic Linux machine   WRF Forum (section Installation)   etcetc...   Among the several tutorials you find online I'd like to add another (useless?) guide (or script) I used to automate the compilation of the WRF/WPS systems, with some special notes. It is definitely a non-exhaustive script and it might have some (hopefully non-destructive errors), so....  >>>use at your own risk<<< The prerequisites are: the netcdf libray (version 4.x) successfully compiled and installed (here in /usr/local) gcc, gfortran (version 5.4.0 or higher) #!/bin/bash # It is supposed that you have the following sources files in $SRCDIR # ARWpost_V3.tar.gz mpic
Continua a leggere

Code published via github and zenodo

Today, thanks to my colleague alfcrisci , I found pretty straightforward to connect some of the several (too many??) accounts I have (github, zenodo, orcid) into a comprehensive house to publish and protect the code I developed making it citable. You get this  by following the instructions you find here. valcap/csv2little_r: First release : this is the first trial... valcap/regens: First release of regens : this is the second trial... Hope it will be useful in the future!
Continua a leggere