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. 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