Created
June 27, 2019 14:25
-
-
Save LarsSchy/9ecb31eb964dd83820c139b2f2769a7c to your computer and use it in GitHub Desktop.
Example script to cut away nodata and create new tile index
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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