prep_toposm_data

#!/bin/bash

if [[ -z $TOPOSM_ENV_SET ]]; then
    echo "Error: TopOSM environment not set."
    exit 1
fi

SRS900913="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over"

# NOTE: X values are negative, since we're west of the 0 meridian
LEFT=$(($XMINLL * -1))
RIGHT=$(($XMAXLL * -1))
BOTTOM=$YMINLL
TOP=$YMAXLL

mkdir -pv $HILLSHADE_DIR
mkdir -pv $COLORMAP_DIR
mkdir -pv $HYPSORELIEF_DIR

for x in `seq $RIGHT $LEFT`; do
for y in `seq $BOTTOM $TOP`; do

    BASE=`printf "n%02dw%03d" $y $x`
    FILE=$BASE.tif
    NEDPATH=$NED13_DIR/$FILE

    if [ -f $NEDPATH ] ; then
    echo; echo "***** $NEDPATH"

    if [[ ! -f $HYPSORELIEF_DIR/$BASE.tif ]] ; then

    # NOTE: First hillshading, then reprojecting seems to yield a better
    # result than vice versa.

    # NOTE: JPEG compression in the TIFFs saves a lot of disk space, but
    # ImageMagick doesn't like it.  Thus, we make JPEG-compressed TIFFs
    # to keep around, but use uncompressed files as inputs to IM and
    # delete them when we're done.

    # Render and reproject hillshade
    echo "Rendering hillshade..."
    echo -n "shade: "
    gdaldem hillshade $NEDPATH $HILLSHADE_DIR/$BASE.unproj.tif -z 0.00001
    echo -n "warp: "
    gdalwarp -t_srs "$SRS900913" -r cubicspline -multi -of GTiff \
        $HILLSHADE_DIR/$BASE.unproj.tif $HILLSHADE_DIR/$BASE.uncomp.tif
    rm $HILLSHADE_DIR/$BASE.unproj.tif
    echo -n "compress: "
    gdal_translate -of GTiff -co COMPRESS=JPEG \
        $HILLSHADE_DIR/$BASE.uncomp.tif $HILLSHADE_DIR/$BASE.tif
    echo -n "overview: "
    gdaladdo -r gauss $HILLSHADE_DIR/$BASE.tif 2 4 8 16 32

    # Render and reproject colormap
    echo "Rendering colormap..."
    echo -n "color: "
    gdaldem color-relief $NEDPATH $COLORFILE $COLORMAP_DIR/$BASE.unproj.tif
    echo -n "warp: "
    gdalwarp -t_srs "$SRS900913" -r cubicspline -multi -of GTiff \
       $COLORMAP_DIR/$BASE.unproj.tif $COLORMAP_DIR/$BASE.uncomp.tif
    rm $COLORMAP_DIR/$BASE.unproj.tif
 
    # Render hypsorelief (relief shade with hypsometric tinting)
    # (IM strips any georeferencing in the file, so we have to use
    # a world file to put it back)
    COLORMAP=$COLORMAP_DIR/$BASE.uncomp.tif
    HILLSHADE=$HILLSHADE_DIR/$BASE.uncomp.tif
    # Translate hillshade to PNG+WLD (just to get a valid .wld file)
    echo -n "world: "
    gdal_translate -of PNG -co WORLDFILE=YES \
        $COLORMAP $HYPSORELIEF_DIR/$BASE.png
    # Create the actual PNG file.  We have to force IM to output RGB, since
    # it'll use a palette if it thinks it can, but JPEG compression only works
    # with RGB images.
    echo -n "compose: "
    convert $COLORMAP -modulate 120 \
        \( $HILLSHADE -level 70,95% +level 0%,80% \) \
	-compose screen -composite \
        \( $HILLSHADE -level 0,75% +level 40%,100% \) \
	-compose multiply -composite \
        -modulate 92 -define png:color-type=2 $HYPSORELIEF_DIR/$BASE.png
    echo "."
    # PNG+WLD -> GeoTIFF
    echo -n "compress: "
    gdal_translate -a_srs "$SRS900913" -of GTiff \
       -co COMPRESS=JPEG \
        $HYPSORELIEF_DIR/$BASE.png $HYPSORELIEF_DIR/$BASE.tif
    rm $HYPSORELIEF_DIR/$BASE.png $HYPSORELIEF_DIR/$BASE.png.aux.xml $HYPSORELIEF_DIR/$BASE.wld
    rm $COLORMAP $HILLSHADE
    echo -n "overview: "
    gdaladdo -r gauss $HYPSORELIEF_DIR/$BASE.tif 2 4 8 16 32
    fi
 
    fi;
done; done

# generate overviews and VRT layers
for d in $HILLSHADE_DIR $HYPSORELIEF_DIR ; do
    pushd . ; cd $d
    gdalbuildvrt all.vrt *.tif
    popd
done

# reproject NLCD data, add overviews, and copy VRT file (with
# the custom color table).
gdalwarp -multi -of GTiff \
    -co ZLEVEL=9 -co COMPRESS=DEFLATE -co PREDICTOR=2 -co BIGTIFF=YES \
    -t_srs EPSG:900913 $NLCD_DIR/nlcd2006_landcover_2-14-11.img \
    $NLCD_DIR/nlcd2006.tif
gdaladdo $NLCD_DIR/nlcd2006.tif 2 4 8 16 32 64
cp nlcd2006.vrt $NLCD_DIR