Skip to content

Instantly share code, notes, and snippets.

@LarsSchy
Created June 27, 2019 14:25
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save LarsSchy/9ecb31eb964dd83820c139b2f2769a7c to your computer and use it in GitHub Desktop.
Save LarsSchy/9ecb31eb964dd83820c139b2f2769a7c to your computer and use it in GitHub Desktop.
Example script to cut away nodata and create new tile index
#!/bin/bash
#
# import_raster_retile_v3.sh
#
# Convert and rewrite tif files with specified compression
# remove nodata parts.
# Create tile index for the new data set
#
# Lars Schylberg, Saab AB
# Date: 2019-01-30
# MIT License
#
# Required modules: gdal/ogr
# mapserver-bin (shptree)
# bc (binary calculator)
# parallel
#
CLEAR='\033[0m'
RED='\033[0;31m'
function usage() {
if [ -n "$1" ]; then
echo -e "${RED}---> $1${CLEAR}\n";
fi
echo "Usage: $0 [-sd source_dir] [-td target_dir] [-jq jpeg_quality] "
echo " -sd, --source_dir Source data directory for raster data"
echo " -td, --target_dir Target direcotory for trimmed data"
echo " -bd, --base_dir Base data directory - will be start for shapepath in Mapserver"
echo " -tip, --tile_index_path Path to build tile index data path that is the start at basedir"
echo " -tmp, --tmp_dir Temporary direcotory for processing"
echo " -ca, --compression Compression algorithm "
echo " -ws, --width_size_x Width size x"
echo " -hs, --height_size_y Heigth size y"
echo " -nd, --nodata_value No data value"
echo " -cl, --clean Clean up temp data"
echo " -q, --quiet Be as quiet as possible"
echo ""
echo "Example: $0 \\"
echo "-sd /data/geodata \\"
echo "-td /data/target/geodata \\"
echo "-bd /data \\"
echo "-tip 50k_raster/Geodata \\"
echo "-tmp /data/tmp \\"
echo "-ca DEFLATE \\"
echo "-ws 5000 \\"
echo "-hs 5000 \\"
echo "-nd 0 \\"
echo "-cl \\"
echo "-q"
exit 1
}
# parse params
if [[ "$#" = 0 ]]; then usage "No parameters are supplied" ; fi
while [[ "$#" > 0 ]]; do case $1 in
-sd|--source_dir) SOURCE_DIR="$2";shift;shift;;
-td|--target_dir) TARGET_DIR="$2";shift;shift;;
-bd|--base_dir) BASE_DIR="$2"; shift;shift;;
-tip|--tile_index_path) TIP="$2";shift;shift;;
-tmp|--tmp_dir) TEMP_DIR="$2";shift;shift;;
-ca|--compression_algorithm) COMPRESSION="$2";shift;shift;;
-ws|--width_size) WIDTH_SIZE="$2";shift;shift;;
-hs|--heght_size) HEIGHT_SIZE="$2";shift;shift;;
-nd|--nodata_value) NODATA_VAL="$2";shift;shift;;
-cl|--clean ) CLEAN=1;shift;;
-q|--quiet) QUIET=1;shift;;
*) usage "Unknown parameter passed: $1"; shift; shift;;
esac; done
# verify params
if [ -z "$SOURCE_DIR" ]; then usage "Source directory is not set."; fi;
if [ -z "$TARGET_DIR" ]; then usage "Target directory is not set."; fi;
if [ -z "$BASE_DIR" ]; then usage "Tile index base directory is not set."; fi;
if [ -z "$TIP" ]; then usage "Tile index data path is not set."; fi;
if [ -z "$TEMP_DIR" ]; then usage "Temp dir is not set."; fi;
if [ -z "$COMPRESSION" ]; then usage "Compression algorithm is not set."; fi;
if [ -z "$WIDTH_SIZE" ]; then usage "Width size is not set."; fi;
if [ -z "$HEIGHT_SIZE" ]; then usage "Height size is not set."; fi;
if [ -z "$NODATA_VAL" ]; then usage "No data value is not specified."; fi;
#
# check for needed commands
#
SCRIPTNAME="${0##*/}"
warn() {
printf >&2 "$SCRIPTNAME: $*\n"
}
iscmd() {
command -v >&- "$@"
}
checkdeps() {
local -i not_found
for cmd; do
iscmd "$cmd" || {
warn $"$cmd is not found"
let not_found++
}
done
(( not_found == 0 )) || {
warn $"Install dependencies listed above to use $SCRIPTNAME"
exit 1
}
}
checkdeps shptree bc parallel
#
# Check gdal version
#
REQ_GDAL_VERSION=2.3.0
function version {
echo "$@" | \
awk -F. '{ printf("%d%03d%03d%03d\n", $1,$2,$3,$4); }'
}
if [ $(version $(gdal-config --version)) -ge $(version "${REQ_GDAL_VERSION}") ]; then
if [ "$QUIET" != "1" ]; then echo "GDAL version is up to date" ; fi;
else
echo "GDAL version is $(gdal-config --version)"
echo "GDAL should be updated to at least ${REQ_GDAL_VERSION}"
exit 1
fi
#
# Start timer
#
begin=$(date +"%s")
#
# Check input directory
#
if [ ! -d ${SOURCE_DIR} ] ; then
echo "Error - Input directory is not found: ${SOURCE_DIR}"
exit 1
fi
### Remove previous results and make temp dirs
# clean up previous runs
if [ -d ${TEMP_DIR} ] ; then
/bin/rm -r ${TEMP_DIR}
fi
TEMP_DIR_1=${TEMP_DIR}/tmp_1
TEMP_DIR_2=${TEMP_DIR}/tmp_2
if [ -d ${TEMP_DIR_1} ] ; then
/bin/rm -r ${TEMP_DIR_1}
fi
mkdir -p $TEMP_DIR_1
if [ -d ${TEMP_DIR_2} ] ; then
/bin/rm -r ${TEMP_DIR_2}
fi
mkdir -p $TEMP_DIR_2
### Create dirs for result, remove previous runs
if [ -d ${TARGET_DIR} ] ; then
/bin/rm -r ${TARGET_DIR}
fi
if [ ! -d ${TARGET_DIR} ] ; then
mkdir -p ${TARGET_DIR}
fi
# place tile index file in a directory Tile_index one level up from the tiles
TILE_INDEX_DIR=$(dirname ${TARGET_DIR})/Tile_index
if [ ! -d ${TILE_INDEX_DIR} ] ; then
mkdir -p ${TILE_INDEX_DIR}
fi
## remove previous index files if they exists
if [ -f ${TILE_INDEX_DIR}/tile_index.shp ] ; then
/bin/rm ${TILE_INDEX_DIR}/*
fi
### Declare some useful functions
function gdal_extent() {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " gdal_extent <input_raster>"
return
fi
EXTENT=$(gdalinfo $1 |\
grep "Upper Left\|Lower Right" |\
sed "s/Upper Left //g;s/Lower Right //g;s/).*//g" |\
tr "\n" " " |\
sed 's/ *$//g' |\
tr -d "[(,]")
echo -n "$EXTENT"
# printf "%q" "$EXTENT"
}
export -f gdal_extent
function ogr_extent() {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " ogr_extent <input_vector>"
return
fi
EXTENT=$(ogrinfo -al -so $1 |\
grep Extent |\
sed 's/Extent: //g' |\
sed 's/(//g' |\
sed 's/)//g' |\
sed 's/ - /, /g')
EXTENT_RESULT=`echo $EXTENT | awk -F ',' '{print $1 " " $4 " " $3 " " $2}'`
echo -n "$EXTENT_RESULT"
# printf "%q" "$EXTENT_RESULT"
}
export -f ogr_extent
function gdal_image_size() {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " gdal_image_size <input_raster>"
return
fi
SIZE=$(gdalinfo $1 |\
grep 'Size is' |\
sed "s/Size is //g" |\
tr -d "[(,)]")
echo -n "$SIZE"
}
export -f gdal_image_size
function gdal_compression() {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " gdal_compression <input_raster>"
return
fi
C_COMPRESSION=$(gdalinfo $1 | \
grep COMPRESSION | \
sed 's/COMPRESSION=//g' | \
sed 's/^[ \t]*//;s/[ \t]*$//')
echo -n "${C_COMPRESSION}"
}
export -f gdal_compression
function tile_info {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " tile_info <inout_raster> <tile_width> <tile_height>"
return
fi
input=$1
tileWidth=$2
tileHeight=$3
read -r imageWidth imageHeight <<<$(gdal_image_size $input)
countTilesX=$((imageWidth/tileWidth))
countTilesY=$((imageHeight/tileHeight))
lastTileWidth=$((imageWidth-(countTilesX*tileWidth)))
lastTileHeight=$((imageHeight-(countTilesY*tileHeight)))
if [[ lastTileWidth -gt 0 ]]; then
countTilesX=$((countTilesX+1))
else
lastTileWidth=$tileWidth
fi
if [[ lastTileHeight -gt 0 ]]; then
countTilesY=$((countTilesY+1))
else
lastTileHeight=$tileHeight
fi
echo "$countTilesX $countTilesY $lastTileWidth $lastTileHeight"
}
export -f tile_info
function retile_with_mask() {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " retile_with_mask <input_raster> <output_dir> <width_size_pixels>< <height_size_pixels> <compression>"
return
fi
input=$1
output_dir=$2
output_name=$(basename -s .tif ${input})
WIDTH_SIZE=$3
HEIGHT_SIZE=$4
COMP=$5
read -r imageWidth imageHeight <<<$(gdal_image_size $input)
if (( $(bc <<<"$imageWidth <= $WIDTH_SIZE && $imageHeight <= $HEIGHT_SIZE") )); then
# echo "just rewriting image with specified compression"
gdal_translate \
-q \
-co TILED=YES \
-co COMPRESS=${COMP} \
$input \
${output_dir}/${output_name}.tif
else
read countTilesX countTilesY lastTileWidth lastTileHeight <<<$(tile_info $input $WIDTH_SIZE $HEIGHT_SIZE)
TX=0; TY=0
for COUNTER_Y in $(seq -w 01 $countTilesY); do
for COUNTER_X in $(seq -w 01 $countTilesX); do
if [ "$COUNTER_X" != "$(printf "%2.0d\n" $countTilesX |sed "s/ /0/")" ]; then
WIDTH_PIXELS=$WIDTH_SIZE
else
WIDTH_PIXELS=$lastTileWidth
fi
if [ "$COUNTER_Y" != "$(printf "%2.0d\n" $countTilesY |sed "s/ /0/")" ]; then
HEIGHT_PIXELS=$HEIGHT_SIZE
else
HEIGHT_PIXELS=$lastTileHeight
fi
# Used when doing jpeg with 3 bands:
# -co INTERLEAVE=PIXEL \
# -b 1 -b 2 -b 3 -mask 4 \
# --config GDAL_TIFF_INTERNAL_MASK YES
gdal_translate \
-q \
-co TILED=YES \
-co COMPRESS=${COMP} \
-co BLOCKXSIZE=512 \
-co BLOCKYSIZE=512 \
-srcwin $TY $TX $WIDTH_PIXELS $HEIGHT_PIXELS \
--config GDAL_CACHEMAX 512 \
$input \
${output_dir}/${output_name}_${COUNTER_Y}_${COUNTER_X}.tif
TY=$((TY+HEIGHT_SIZE))
done
TX=$((TX+WIDTH_SIZE))
TY=0
done
fi
}
export -f retile_with_mask
function copy_tif_files() {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " copy_tif_files <input_raster_dir> <output_dir>"
return
fi
INPUT_DIR=$1
OUTPUT_DIR=$2
for file in ${INPUT_DIR}/*.tif
do
file_base=$(basename -s .tif ${file})
cp $file ${OUTPUT_DIR}/${file_base}.tif
done
}
export -f copy_tif_files
function remove_empty_tile() {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " remove_empty_tile <input_raster_dir> <nodata_value>"
return
fi
INPUT_FILE=$1
ND=$2
# if there is a nodata value in the dataset we use that one,
# otherwise we use the nodata value from the command line
NODATA_VALUE=$(gdalinfo ${INPUT_FILE} | \
grep 'NoData Value=' | \
sed 's/NoData Value=//g' | \
awk '{$1=$1;print}')
if [ -n "$NODATA_VALUE" ]; then
gdal_edit.py -unsetnodata ${INPUT_FILE}
STDDEV=$(gdalinfo ${INPUT_FILE} -stats | \
grep STATISTICS_STDDEV | \
sed 's/STATISTICS_STDDEV=//g' | \
awk '{$1=$1;printf "%d", $1}')
else
STDDEV=$(gdalinfo ${INPUT_FILE} -stats | \
grep STATISTICS_STDDEV | \
sed 's/STATISTICS_STDDEV=//g' | \
awk '{$1=$1;printf "%d", $1}')
NODATA_VALUE=$ND
fi
if [ "$STDDEV" = "0" ]; then
/bin/rm ${INPUT_FILE}
else
gdal_edit.py -a_nodata ${NODATA_VALUE} ${INPUT_FILE}
fi
}
export -f remove_empty_tile
function trim_tile() {
if [ -z "$1" ]; then
echo "Missing arguments. Syntax:"
echo " remove_empty_tile <input_file> <output_dir> <compression>"
return
fi
INPUT_FILE=$1
OUTPUT_DIR=$2
COMP=$3
tile_base=$(basename -s .tif ${INPUT_FILE})
tile_dir=$(dirname ${INPUT_FILE})
gdal_calc.py \
-A ${INPUT_FILE} \
--outfile=${tile_dir}/${tile_base}_calc.tif \
--calc="1*(A>=0)" \
--NoDataValue=0 \
--quiet
gdal_polygonize.py \
-q \
${tile_dir}/${tile_base}_calc.tif \
-f "ESRI Shapefile" \
${tile_dir}/${tile_base}_calc.shp
## Just trim if polygon extent is smaller, otherwise copy file to target dir
read -r v1 v2 v3 v4 <<<$(ogr_extent ${tile_dir}/${tile_base}_calc.shp)
read -r r1 r2 r3 r4 <<<$(gdal_extent ${INPUT_FILE})
if (( $(bc <<<"$v1 == $r1 && $v2 == $r2 && $v3 == $r3 && $v4 == $r4") )); then
cp ${INPUT_FILE} ${OUTPUT_DIR}/${tile_base}.tif
else
gdal_translate \
-q \
-projwin $v1 $v2 $v3 $v4 \
-co TILED=YES \
-co BLOCKXSIZE=512 \
-co BLOCKYSIZE=512 \
-co COMPRESS=${COMP} \
--config GDAL_CACHEMAX 512 \
-q \
${INPUT_FILE} \
${OUTPUT_DIR}/${tile_base}.tif
gdaladdo \
-q \
-r nearest \
--config COMPRESS_OVERVIEW ${COMP} \
--config GDAL_TIFF_INTERNAL_MASK YES \
--config GDAL_TIFF_OVR_BLOCKSIZE 256 \
--config TILED YES \
${OUTPUT_DIR}/${tile_base}.tif
fi
}
export -f trim_tile
##### Process raster data ##############################################
# Copy files to first temp dir
if [ "$QUIET" != "1" ]; then echo "Copy files to first temp dir" ; fi;
copy_tif_files ${SOURCE_DIR} ${TEMP_DIR_1}
## check if image contains data, otherwise remove it in first temp dir
if [ "$QUIET" != "1" ]; then echo "Remove empty images" ; fi;
## SHELL=$(type -p bash) parallel -j +1 --will-cite
SHELL=$(type -p bash) parallel -j +1 \
remove_empty_tile \
::: ${TEMP_DIR_1}/*.tif \
::: ${NODATA_VAL}
if [ "$QUIET" != "1" ]; then echo "Start retiling" ; fi;
## SHELL=$(type -p bash) parallel -j +1 --will-cite
SHELL=$(type -p bash) parallel -j +1 \
retile_with_mask \
::: ${TEMP_DIR_1}/*.tif \
::: ${TEMP_DIR_2} \
::: ${WIDTH_SIZE} \
::: ${HEIGHT_SIZE} \
::: ${COMPRESSION}
if [ "$QUIET" != "1" ]; then echo "Remove empty tiles" ; fi;
## SHELL=$(type -p bash) parallel -j +1 --will-cite
SHELL=$(type -p bash) parallel -j +1 \
remove_empty_tile \
::: ${TEMP_DIR_2}/*.tif \
::: ${NODATA_VAL}
if [ "$QUIET" != "1" ]; then echo "Trim tiles" ; fi;
## SHELL=$(type -p bash) parallel -j +1 --will-cite
SHELL=$(type -p bash) parallel -j +1 \
trim_tile \
::: ${TEMP_DIR_2}/*.tif \
::: ${TARGET_DIR} \
::: ${COMPRESSION}
#######################################################################
# Create Mapserver raster image index files
# Move to BASE_DIR to get the tile index paths right for mapserver
BACK_DIR=$(pwd)
cd $BASE_DIR
if [ "$QUIET" != "1" ]; then echo "Craate mapserver image index files"; fi;
if [ "$QUIET" != "1" ]; then
gdaltindex ${TILE_INDEX_DIR}/tile_index.shp ${TIP}/*.tif
else
gdaltindex ${TILE_INDEX_DIR}/tile_index.shp ${TIP}/*.tif 1>/dev/null
fi
if [ "$QUIET" != "1" ]; then
shptree ${TILE_INDEX_DIR}/tile_index.shp
else
shptree ${TILE_INDEX_DIR}/tile_index.shp 1>/dev/null
fi
cd $BACK_DIR
## Remove intermidiate results
if [ "$CLEAN" == "1" ]; then /bin/rm -r ${TEMP_DIR}; fi;
#######################################################################
#
# Print script execution time:
#
termin=$(date +"%s")
difftimelps=$(($termin-$begin))
echo "Script Execution time: $(($difftimelps / 60)) minutes \
and $(($difftimelps % 60)) seconds."
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment