sentinel2_pipeline/bin/s2_clip

#!/bin/bash

# all outputs can be found in /OUTDIR/log.txt
# creates a .tif named INPUTDATE_l2a.tif
# also returns the path to that .tif file in console

# [ ] weird error while creating the tif: "Warning 1: Unable to export GeoTIFF file with different datatypes per different bands. All bands should have the same types in TIFF."
#     but works anyway

INDIR=''
OUTDIR=''
AOI=''

while test $# -gt 0; do
  case "$1" in
    -h|--help)
      echo "s2_clip - sentinel 2 change detection stitching and clipping"
      echo " "
      echo "s2_clip [options]"
      echo " "
      echo "options:"
      echo "-h, --help                show brief help."
      echo "-a, --aoi                 area of interest as geojson polygon geometry json file. Required."
      echo "-i, --input               space separated list of paths to input SAFE folders. Required."
      echo "-o, --outputdir           path to output directory. Required."
      exit 0
      ;;
    -i|--input)
      shift
      INDIR=$1
      shift
      while [[ ! ${1:0:1} == "-" ]] && [[ ! -z $1 ]];
      do
        INDIR="$INDIR $1"
        shift
      done
      ;;
    -o|--outputdir)
      shift
      OUTDIR=$1
      shift
      ;;
    -a|--aoi)
      shift
      AOI=$1
      shift
      ;;
    *)
      break
      ;;
  esac
done

# check whether input directories & date exist

if [ ! -d "$OUTDIR" ]
  then echo "output directory doesn't exist!"
  exit
fi

if [ ! -f "$AOI" ]
  then echo "aoi file does not exist"
  exit
fi

# make sure OUTDIR doesn't have an appended '/'
OUTDIR=${OUTDIR%/}

# create log file (currently one per directory, not unique)
log=$OUTDIR/log.txt
touch $log
timestamp=$(date +%Y%m%d_%H%M%S)
echo "###script started with INDIR = $INDIR & OUTDIR = $OUTDIR; timestamp: $timestamp" 1>> $log 2>> $log

# create vrtlist to save vrt we need to clip
vrtlist=$OUTDIR/vrt_$timestamp.txt
touch $vrtlist

# call checkfolder to create list of tiles with matching date +-1 day
TILELIST=$INDIR
# if TILELIST isn't empty
if [ "$TILELIST" != "" ]
then
  l2a_tiles=$TILELIST
else echo "no matching tiles for that date ($INPUTDATE) found" 1>> $log 2>> $log
  exit
fi

# output tiff is named inputfolder_inputdate_l2a.tiff
DATE=$(basename ${l2a_tiles%% *})
TIFFNAME=${DATE:11:8}
stitchedpath="${OUTDIR}/${TIFFNAME}_l2a"

function cleanup {
  # cleanup
  [ -f $vrtlist ] && xargs rm -f < $vrtlist
  rm -f $vrtlist ${stitchedpath}.vrt ${stitchedpath}_tmp.tif
  rm -f $log
}
trap cleanup EXIT

function tile_bands {
  # given a path to a MSIL2A tile, a list of all atomic .jp2 band files is returned.
  # where bands are available in multiple resolutions, the highest is used.
  # band order is: 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 8A (note that band 10 is missing!)
  bandlist=`find $tile -path '*IMG_DATA*.jp2' -type f ! -name *AOT* ! -name *WVP* ! -name *TCP* ! -name *TCI* ! -name *SCL* ! -name *B02_20* ! -name *B03_20* ! -name *B04_20* ! -name *B02_60* ! -name *B03_60* ! -name *B04_60* ! -name *B05_60* ! -name *B06_60* ! -name *B07_60* ! -name *B10_60* ! -name *B11_60* ! -name *B12_60* ! -name *B8A_60*`
  echo `python2 -c "import sys; a = sys.argv[1:]; a.sort(key = lambda x: x.split('/')[-1]); print ' '.join(a)" $bandlist`
}

# check if file already exists. If so skip process
if [ -f "${stitchedpath}.tif" ]; then
  echo "${stitchedpath}.tif already exists. Skipping clipping." 1>&2
  echo "${stitchedpath}.tif"
  exit
fi

# for every tile, do
for tile in $l2a_tiles
do
  tilename=$(basename ${tile})
  outtilename=$OUTDIR/$tilename
  echo "-working on $tilename" >> $log

  # filter for bands that we need
  # bandlist is empty when folder doesn't contain an IMG_DATA directory
  bandlist=`tile_bands $tile`

  if [ -z "$bandlist" ]; then
    >&2 echo "skipping invalid tile `basename $tile`"
    continue
  fi

  # build vrt for every tile, write name to $vrtlist for clipping function
  gdalbuildvrt -separate -resolution highest -overwrite $outtilename.vrt $bandlist 1>> $log 2>> $log
  echo $outtilename.vrt >> $vrtlist
done
echo "done with tiles" >> $log

if [ ! -s $vrtlist ]; then
  exit # if we don't have any valid bands for this orbit
fi

# Stitching
echo "start stitching" >> $log
echo "stitching `cat $vrtlist`..." 1>&2
gdalbuildvrt -resolution highest -input_file_list $vrtlist -overwrite $stitchedpath.vrt 1>> $log 2>> $log
gdal_translate -of GTiff $stitchedpath.vrt ${stitchedpath}_tmp.tif 1>> $log 2>> $log
echo "done\n" >> $log

# CLIP image
echo "start clipping" >> $log
echo "clipping $stitchedpath.vrt..." 1>&2
gdalwarp -t_srs "EPSG:4326" -dstnodata 0 -crop_to_cutline -cutline $AOI ${stitchedpath}_tmp.tif -overwrite $stitchedpath.tif 1>> $log 2>> $log
echo "done" >> $log
echo "" >> $log

# output required for pipeline
echo $stitchedpath.tif