a somewhat refined milti-line example using gdal and gmt tools from CLI:
##Death Valley National Park, CA, USA
##(Mesquite Flats area, -R-117.33/-116.84/36.51/36.84). Steep mountains and large alluvial fans:
#ROI="-R-117.33/-116.84/36.51/36.84"
## Hawaii, -R-156/-155.5/19/19.5
#ROI="-R-156/-155.5/19/19.5"
## Alp mountains -R10.13/10.43/45.46/45.85
#ROI="-R10.13/10.43/45.46/45.85"
## Grand Canyon "-R-112.06/-111.60/35.97/36.4
ROI="-R-112.06/-111.60/35.97/36.4"
gmt makecpt -Cgray -T0/65535/1 -N > textureshade.cpt # texture_image produces a uint16 GeoTIFF image
gmt makecpt -Cgray -T0/255/1 -N > hillshade.cpt # gdaldem hillshade produces a uint8 GeoTIFF image
gmt grdcut @srtm_relief_03s -Ghawaii_south.nc $ROI
gmt grdcut @srtm_relief_03s -Ghawaii_south.flt=gd:EHdr $ROI
DETAIL=1.0
CONTRAST=+2.0
# texture and texture_image must be compiled from source
Source_Code/texture $DETAIL hawaii_south.flt hawaii_south_texture.flt
Source_Code/texture_image $CONTRAST hawaii_south_texture.flt hawaii_south_textureshade.tif
# embedding projection info into the textureshade image
gdal_edit.py -a_srs hawaii_south_textureshade.prj hawaii_south_textureshade.tif
gmt grdcut @srtm_relief_03s -Ghawaii_south.nc $ROI
# creating hillshade GeoTIFF
VERTICAL_ENHANCEMENT=4.0
gdaldem hillshade hawaii_south.flt hawaii_south_hillshade.tif -s 111120 -z $VERTICAL_ENHANCEMENT -compute_edges
# removing NoData value of 0 set by "gdaldem hillshade" by default
# setting NoData makes no sense as texture shading software anyway does not support data with NaNs
# processing NaN values set by gdaldem and gdal_calac.py may create NaN artefacts on the resulting plot by gmt
gdal_edit.py -unsetnodata hawaii_south_hillshade.tif
# blend hillshade with texture shade
A=0.4 # A=0.6 means the result is 60% texture shade value plus 40% hillshade value
gdal_calc.py -A hawaii_south_textureshade.tif -B hawaii_south_hillshade.tif --allBands=B \
--calc="uint8( ( \
$A * (A/65535.) + (1.0 - $A) * (B/255.) \
) * 255 )" --type=Byte --outfile=bat_final.tif --overwrite --hideNoData
# removing NoData value of 255 set by gdal_calc.py
gdal_edit.py -unsetnodata bat_final.tif
# plotting source reilef, shaded relief, texture shade, and hillshade blended with that texture shade
J=-Jx6i
gmt begin texthill jpg -C
gmt set MAP_FRAME_TYPE=plain FORMAT_GEO_MAP="ddd.x" FONT_ANNOT_PRIMARY=auto
gmt subplot begin 2x2 -Fs9c
gmt subplot set 0
gmt grdimage hawaii_south.nc $J -BWSne
gmt colorbar $J -DJBC+h -By+lm -Bxaf
gmt subplot set 1
gmt grdimage hawaii_south_hillshade.tif -Chillshade.cpt $J -BWSne
gmt subplot set 2
gmt grdimage hawaii_south_textureshade.tif -Ctextureshade.cpt $J -BWSne
gmt subplot set 3
gmt grdimage bat_final.tif -Chillshade.cpt $J -BWSne
gmt subplot end
gmt end
result:
now with color blended, color layer produced using gmt grdimage -Aoutput.tif:
##Death Valley National Park, CA, USA
##(Mesquite Flats area, -R-117.33/-116.84/36.51/36.84). Steep mountains and large alluvial fans:
#ROI="-R-117.33/-116.84/36.51/36.84"
#CPT=gray
## Hawaii, -R-156/-155.5/19/19.5
#ROI="-R-156/-155.5/19/19.5"
#CPT=gray
#CPT=faa.cpt
#unzip faa.cpt.zip faa.cpt
### Alp mountains -R10.13/10.43/45.46/45.85
ROI="-R10.13/10.63/45.46/45.85"
CPT=geo
## Grand Canyon
#ROI="-R-112.06/-111.60/35.97/36.25"
#CPT=geo
#CPT=gray
if [[ -z $ROI || -z $CPT ]]; then
printf "ROI or CPT was not indicated, terimnating\n"
exit 0
fi
J=-Jx6i
gmt makecpt -Cgray -T0/65535/1 -N > textureshade.cpt # texture_image produces a uint16 GeoTIFF image
gmt makecpt -Cgray -T0/255/1 -N > hillshade.cpt # gdaldem hillshade produces a uint8 GeoTIFF image
gmt grdcut @srtm_relief_03s -Ghawaii_south.nc $ROI
gmt grdcut @srtm_relief_03s -Ghawaii_south.flt=gd:EHdr $ROI
DETAIL=1.0
CONTRAST=+2.0
# texture and texture_image must be compiled from source
Source_Code/texture $DETAIL hawaii_south.flt hawaii_south_texture.flt
Source_Code/texture_image $CONTRAST hawaii_south_texture.flt hawaii_south_textureshade.tif
# embedding projection info into the textureshade image
gdal_edit.py -a_srs hawaii_south_textureshade.prj hawaii_south_textureshade.tif
# creating hillshade GeoTIFF
#VERTICAL_ENHANCEMENT="4.0 -combined"
#VERTICAL_ENHANCEMENT="4.0 -multidirectional"
VERTICAL_ENHANCEMENT="4.0"
gdaldem hillshade hawaii_south.flt hawaii_south_hillshade.tif -s 111120 -z $VERTICAL_ENHANCEMENT -compute_edges
# removing NoData value of 0 set by "gdaldem hillshade" by default
# setting NoData makes no sense as texture shading software anyway does not support data with NaNs
# processing NaN values set by gdaldem and gdal_calac.py may create NaN artefacts on the resulting plot by gmt
gdal_edit.py -unsetnodata hawaii_south_hillshade.tif
# blend hillshade with texture shade
A=0.4 # A=0.6 means the result is 60% texture shade value plus 40% hillshade value
gdal_calc.py -A hawaii_south_textureshade.tif -B hawaii_south_hillshade.tif --allBands=B \
--calc="uint8( ( \
$A * (A/65535.) + (1.0 - $A) * (B/255.) \
) * 255 )" --type=Byte --outfile=bat_final.tif --overwrite --hideNoData
# removing NoData value of 255 set by gdal_calc.py
gdal_edit.py -unsetnodata bat_final.tif
TEXTHILL="$A*texture shade + (1 - $A)*hillshade"
# generating color relief using grdimage
gmt grdimage hawaii_south.nc -C$CPT -Ahawaii_south_color_relief.tif $J
# blending texture+hillshade with color relief
A=0.6 # A=0.6 means the result is 60% texture+hill shade value plus 40% color relief value
gdal_calc.py -A bat_final.tif -B hawaii_south_color_relief.tif --allBands=B \
--calc="uint8( ( \
$A * (A/255.) + (1.0 - $A) * (B/255.) \
) * 255 )" --type=Byte --outfile=hillshade-overlay.tif --overwrite --hideNoData
COLORTEXTHILL="$A*texture+hillshade + (1 - $A)*color"
# plotting source reilef, shaded relief, texture shade, and
# hillshade (grdimage), hillshade (gdal) blended with that texture shade, and all three blended.
gmt begin texthill jpg -C
gmt set MAP_FRAME_TYPE=plain FORMAT_GEO_MAP="ddd.x" FONT_ANNOT_PRIMARY=auto
gmt set FONT_TITLE 10p
gmt subplot begin 2x3 -Fs9c
gmt subplot set 0
TITLE="source DEM"
gmt grdimage hawaii_south.nc $J -BWSne+t$"$TITLE" -Bxa0.1f0.02 -Bya0.1f0.02
gmt colorbar $J -DJBC+h -By+lm -Bxaf
gmt subplot set 1
TITLE="hillshade: gdaldem hillshade -zfactor=$VERTICAL_ENHANCEMENT"
gmt grdimage hawaii_south_hillshade.tif -Chillshade.cpt $J -BWSne+t"$TITLE" -Bxa0.1f0.02 -Bya0.1f0.02
gmt subplot set 2
TITLE="texture shade, detail=$DETAIL, contrast=$CONTRAST"
gmt grdimage hawaii_south_textureshade.tif -Ctextureshade.cpt $J -BWSne+t"$TITLE" -Bxa0.1f0.02 -Bya0.1f0.02
gmt subplot set 3
TITLE="color: grdimage -C$CPT"
gmt grdimage hawaii_south.nc $J -C$CPT -BWSne+t"$TITLE" -Bxa0.1f0.02 -Bya0.1f0.02
gmt subplot set 4
TITLE=$TEXTHILL
gmt grdimage bat_final.tif -Chillshade.cpt $J -BWSne+t"$TITLE" -Bxa0.1f0.02 -Bya0.1f0.02
gmt subplot set 5
TITLE=$COLORTEXTHILL
gmt grdimage hillshade-overlay.tif $J -BWSne+t"$TITLE" -Bxa0.1f0.02 -Bya0.1f0.02
gmt subplot end
gmt end
rm -f hawaii_south.* hawaii_south_texture.* hawaii_south_textureshade* hawaii_south_hillshade* *.cpt
1 Like
