A script to quickly create a truecolor image from any Landsat scene ID

landsat_speedrun.sh

# Processing Landsat 8 data into truecolor PXM compatibility
# (pansharpening note: a lot of scenes evaluated for this were straight-up missing band 8 :shrug:)
# Usage: bash landsat_speedrun.sh <Landsat scene ID, e.g. "LC08_L1TP_006038_20230925_20231002_02_T1">

# Get the scene ID
# from https://earthexplorer.usgs.gov/
SCENE=$1

# Extract components from SCENE
L_PATH=$(echo $SCENE | cut -c 11-13)
L_ROW=$(echo $SCENE | cut -c 14-16)
YEAR=$(echo $SCENE | cut -c 18-21)

# Collect data from AWS
BANDS=( 2 3 4 8 )
for BAND in ${BANDS[@]}; do
    aws s3 cp \
        --request-payer requester \
        s3://usgs-landsat/collection02/level-1/standard/oli-tirs/${YEAR}/${L_PATH}/${L_ROW}/${SCENE}/${SCENE}_B${BAND}.TIF \
        ${SCENE}_B${BAND}.TIF
done

# Combine and pansharpen for truecolor
gdal_pansharpen.py \
    ${SCENE}_B8.TIF \
    ${SCENE}_B{4,3,2}.TIF \
    ${SCENE}_sharp.tif \
    -r bilinear \
    -co COMPRESS=DEFLATE \
    -co PHOTOMETRIC=RGB

# Determine color min/max thresholds for stretching
eval $(gdalinfo -mm ${SCENE}_sharp.tif -json | jq -r '
  .bands |
  map(.computedMin) as $mins |
  map(.computedMax) as $maxes |
  "MIN=\($mins | min) MAX=\($maxes | max)"
')

# Smoosh colors into 8-bit range, add "pep"
gdal_translate \
    ${SCENE}_sharp.tif \
    ${SCENE}_scaled_255.tif \
    -scale ${MIN} ${MAX} 0 255 \
    -exponent 0.5 \
    -ot Byte \
    -co COMPRESS=DEFLATE \
    -co PHOTOMETRIC=RGB

# clean up everything except the final image
rm ${SCENE}_B{2,3,4,8}.TIF ${SCENE}_sharp.tif