gistfile1.sh

#!/usr/bin/env bash
#
# clip_usgs_to_shape.sh
# Clip individual band files from USGS-distributed remote-senisng images
# to extent of polygon provided in a shapefile, using the GDAL command line
# tools.
#
# Chris Waigl, 2014-04-15
#
# Input:
#  * search pattern for imagery scene directories in current directory
#  * path to shapefile
# 
# Usage:
# clip_usgs_to_shape.sh [search_pattern] [path_to_shapefile.shp]
# 
# Assumptions:
#  * GDAL command line tools are installed (duh).
#  * All GeoTIFF files, assumed to be all with a .TIF or .tif extension,
#    within the target directories will be clipped. This includes those
#    that come from custom processing after unzipping the scene archive
#  * For each GeoTIFF file, a new file called [basename]_clip.tif will
#    be generated.
#  * The shapefile and the GeoTIFF files are both in the same coordinate
#    reference system. Usually this means the shapefile has to be converted
#    to UTM (appropriate zone) first, using GDAL/OGR. For example for UTM6N on 
#    the WGS84 datum, pre-process the shapefile like this:
#        ogr2ogr -f "ESRI Shapefile" polyUTM6.shp poly.shp -t_SRS EPSG:32606
#  * The command is run from the directory immediately above the scene 
#    direcories. The directory strucutre looks somewhat like this:
#    ├── EO1Spath_row_day...
#    │   ├── EO1Spath_row_day..._B01_L1GST.TIF
#    │   ├── EO1Spath_row_day..._B02_L1GST.TIF
#    │   ├── EO1Spath_row_day..._B03_L1GST.TIF
#    │   ├── ...
#    │   ├── EO1Spath_row_day..._processed1.tif
#    │   ├── EO1Spath_row_day..._processed2.tif
#    │   ├── EO1Spath_row_day..._MTL_L1GST.TXT
#    │   ├── EO1Spath_row_day..._PF1_01.fgdc
#    │   └── README.txt
#    ├── LTnpath_row_day...
#    │   ├── LTnpath_row_day..._B1.TIF
#    │   ├── LTnpath_row_day..._B2.TIF
#    │   ├── LTnpath_row_day..._B3.TIF
#    │   ├── ...
#    │   ├── LTnpath_row_day..._GCP.txt
#    │   ├── LTnpath_row_day..._MTL.txt
#    │   ├── LTnpath_row_day..._VER.jpg
#    │   ├── LTnpath_row_day..._VER.txt
#    │   └── README.GTF

# Check correct number of arguments
if [ $# -ne 2 ]; then
    echo "*** This script requires two arguments. Usage:"
    echo "*** $0 [search_pattern] [path_to_shapefile.shp]"
    exit 1
fi

# Check whether shapefile exist
if [ ! -f $2 ]; then
    echo "*** Shapefile not found: $2"
    exit 1
fi

# Assemble directories
PATTERN=$1
SHPFILE=$2
STARTDIR=$(pwd)
DIRLIST=$(find . -name "*$1*" -mindepth 1 -maxdepth 1 -type d)
if [ -z $DIRLIST ]; then
    echo "*** No directories found that match the pattern $1"
    exit 1
fi

# Calculate extent of clipping 
BASE=`basename $SHPFILE .shp`
EXTENT=`ogrinfo -so $SHPFILE $BASE | grep Extent \
| sed 's/Extent: //g' | sed 's/(//g' | sed 's/)//g' \
| sed 's/ - /, /g'`
EXTENT=`echo $EXTENT | awk -F ',' '{print $1 " " $4 " " $3 " " $2}'`
# Testing Wood River 2009 in Hyperion, result should be:
# EXTENT="448286.298777 7161025.255319 465335.082492 7132951.853959"

for DIR in $DIRLIST
do
    echo "*** Working in directory ${DIR}."
    cd ${DIR}
    # support both Landsat and EO-1, upper and lowercase extension
    FILELIST=`ls [LE]*.[Tt][Ii][Ff]`
    for TIFFILE in $FILELIST
    do
        # do not clip already clipped files
        CHKCLIP=$(echo ${TIFFILE} | grep "_clip\.tif$")
        if [[ ! -z "${CHKCLIP}" ]]; then
            echo "*** Skipping ${TIFFILE}"
            continue
        fi
        NEWNAME="$(basename ${TIFFILE%.*})_clip.tif"
        echo "*** Trying to clip $TIFFILE -> $NEWNAME"
        if [ -f $NEWNAME ]; then
            read -p "$NEWNAME already exists. Overwrite? [yn]" yn
            case $yn in
                [Yy]* ) ;;
                [Nn]* ) NEWNAME="$(basename ${NEWNAME} .tif)_clip.tif";;
                * ) echo "Please answer yes or no.";;
            esac
        fi
        gdal_translate -projwin $EXTENT -of GTiff $TIFFILE $NEWNAME
    done
    cd $STARTDIR
done