Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
# Documented at
# First we fill holes in the DEM:
gdal_fillnodata dem30.img dem30_fill.tif
# Now we are going to create a blurred version of the DEM. We'll use this as the generalized DEM
# in our blended DEM
gdalwarp -tr 400 400 -r lanczos -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" dem30_fill.tif dem30_400.tif
gdalwarp -tr 30.929927379942203 30.929927379942203 -r lanczos -co "BLOCKXSIZE=512" -co "BLOCKYSIZE=512" dem30_400.tif dem30_400-30.tif
# We also need a more generalized DEM to for our TPI work.
# See for more info
gdalwarp -tr 5000 5000 -r lanczos dem30_fill.tif dem30_5000.tif
gdalwarp -tr 30.929927379942203 30.929927379942203 -r lanczos dem30_5000.tif dem30_5000-30.tif
# To avoid issues with differently sized images when we use gdal_calc, we'll stack our images into a single image
gdal_merge -separate -of GTiff -o tpistack.tif dem30_5000-30.tif dem30_fill.tif
# Now we calculate a normalized TPI (effectively calculating Standard Deviation over 5000 meter radius)
gdal_calc -A tpistack.tif -B tpistack.tif --A_band=1 --B_band=2 --outfile=tpi.tif --calc="sqrt((B * 1.00000 - A * 1.00000) * (B * 1.00000 - A * 1.00000))"
# Same trick here -- we stack our different images to ensure we have our image dimensions equal
# for gdal_calc
gdal_merge -separate -of GTiff -o stack.tif tpi.tif dem30_400-30.tif dem30_fill.tif
# And now we calculate our weighted DEM
gdal_calc -A stack.tif -B stack.tif -C stack.tif --A_band=2 --B_band=1 --C_band=3 --outfile=result.tif --calc="((A * B) + (C * (288 - B))) / 288"
# We'll convert this to a PNG for use in Povray. See
gdal_translate -ot UInt16 -of PNG result.tif pov-viewshed\result.png
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment