Landsat 8 TOA using data on S3 http://landsat.usgs.gov/Landsat8_Using_Product.php

l8-toa.sh

#!/bin/bash
set -eou pipefail
 
 
function usage() {
  echo ""
  echo "For a given L8 scene, get the top of atmosphere RGB image."
  echo ""
  echo "Usage: ./l8-toa.sh [entityId] [dstpath]"
  echo ""
  echo "e.g. ./l8-toa.sh LC81200282014306LGN00 LC81200282014306LGN00_TOA.TIF"
  echo ""
}


entityId=${1:-}
dstpath=${2:-}


if [[ -z $dstpath ]]; then
  usage;
  exit 1;
fi


function onerror() {
  
  if [ $? -eq 2 ]
  then
    echo ""
    echo "awscli and rasterio are required to run this script."
    echo ""
    echo "  pip install awscli"
    echo "  pip install rasterio"
    echo "  apt-get install parallel (or brew install parallel)"
    echo ""
  else
    usage
  fi

  exit 1
}
trap onerror EXIT


function check_dependencies() {

  dependencies=( parallel aws rio jq )
  for dependency in "${dependencies[@]}"
  do
    (command -v $dependency > /dev/null) || return 2
  done
}


check_dependencies


function get_file() {
  local entityId=$1; local band=$2;

  if [ ! -f  ${entityId}_B${band}.TIF ]; then

    path=${entityId:3:3}
    row=${entityId:6:3}

    aws s3 cp s3://landsat-pds/L8/${path}/${row}/${entityId}/${entityId}_B${band}.TIF ${entityId}_B${band}.TIF
  fi
}
export -f get_file


function get_landsat_toa() {
  local entityId=$1; local dstpath=$2;

  path=${entityId:3:3}
  row=${entityId:6:3}

  parallel get_file $entityId {} ::: `seq 2 4`
  mtl=$(aws s3 cp s3://landsat-pds/L8/${path}/${row}/${entityId}/${entityId}_MTL.json -| cat)

  scale4=$(echo $mtl | jq -r ".L1_METADATA_FILE.RADIOMETRIC_RESCALING.RADIANCE_MULT_BAND_4")
  scale3=$(echo $mtl | jq -r ".L1_METADATA_FILE.RADIOMETRIC_RESCALING.RADIANCE_MULT_BAND_3")
  scale2=$(echo $mtl | jq -r ".L1_METADATA_FILE.RADIOMETRIC_RESCALING.RADIANCE_MULT_BAND_2")

  offset4=$(echo $mtl | jq -r ".L1_METADATA_FILE.RADIOMETRIC_RESCALING.RADIANCE_ADD_BAND_4")
  offset3=$(echo $mtl | jq -r ".L1_METADATA_FILE.RADIOMETRIC_RESCALING.RADIANCE_ADD_BAND_3")
  offset2=$(echo $mtl | jq -r ".L1_METADATA_FILE.RADIOMETRIC_RESCALING.RADIANCE_ADD_BAND_2")

  rio calc --dtype float32 "(+ $offset4 (* $scale4 (read 1)))" ${entityId}_B4.TIF ${entityId}_TOA_B4.TIF
  rio calc --dtype float32 "(+ $offset3 (* $scale3 (read 1)))" ${entityId}_B3.TIF ${entityId}_TOA_B3.TIF
  rio calc --dtype float32 "(+ $offset2 (* $scale2 (read 1)))" ${entityId}_B2.TIF ${entityId}_TOA_B2.TIF

  rio stack --rgb ${entityId}_TOA_B4.TIF ${entityId}_TOA_B3.TIF ${entityId}_TOA_B2.TIF --output tmp.tif
  
  min4=$(echo $mtl | jq -r ".L1_METADATA_FILE.MIN_MAX_RADIANCE.RADIANCE_MINIMUM_BAND_4")
  min3=$(echo $mtl | jq -r ".L1_METADATA_FILE.MIN_MAX_RADIANCE.RADIANCE_MINIMUM_BAND_3")
  min2=$(echo $mtl | jq -r ".L1_METADATA_FILE.MIN_MAX_RADIANCE.RADIANCE_MINIMUM_BAND_2")
  
  max4=$(echo $mtl | jq -r ".L1_METADATA_FILE.MIN_MAX_RADIANCE.RADIANCE_MAXIMUM_BAND_4")
  max3=$(echo $mtl | jq -r ".L1_METADATA_FILE.MIN_MAX_RADIANCE.RADIANCE_MAXIMUM_BAND_3")
  max2=$(echo $mtl | jq -r ".L1_METADATA_FILE.MIN_MAX_RADIANCE.RADIANCE_MAXIMUM_BAND_2")
  
  gdal_translate -ot Uint16 -co photometric=rgb -scale 0 200 0 65535 tmp.tif $dstpath
  rm tmp.tif
}


get_landsat_toa $entityId $dstpath


trap - EXIT
exit 0