Tweaking illumination/shading

KristofKoch · March 30, 2021, 7:25am

Hi @John_Geodesy, thank you for your explanation – I’ll have a closer look at -Nt. It sounds promising. I was experimenting with -Em135/45.

But first to the Final Frontier that are color spaces and shading. It feels more like the Abyss. I switched to HSV in the CPT and added #COLOR_MODEL=+HSV as per @pwessel’s tip. I observed only minimal differences if at all. Next step was a look at the extremes. How does an illumination at +1 and -1 look like? First no illumination at all:

Just the plain colors from the CPT. Next step: Everything at +1:

Definitely lighter than 0 but also quite colourful. On to the dark side with -1:

Darker, yes, but also very deep red and greens. I expected something like this mockup for -1:

It was created in a graphics editor by using the 0 image and adding a black layer with 33% transparency on top of the land (normal blend mode).

Side by side comparison with 0 superimposed:

For me, the left side appears to “glow” and the right side feels more like “shadow”.

What do you think? Is there a way to circumvent the color shift in GMT and get the result on the right?

Code to play with:

land_polygons_osm_hawaii.bf2.zip (104.8 KB)
faa_hsv.cpt.zip (815 Bytes)

gmt grdcut @srtm_relief_03s -Ghawaii_south.grd -R-156/-155.5/19/19.5
gmt grdgradient hawaii_south.grd -Gshading.grd -Em135/45
gmt grdmath shading.grd 1 GE = shading_0.grd
gmt grdmath shading_0.grd 1 ADD = shading_1.grd
gmt grdmath shading_1.grd -1 MUL = shading_-1.grd



gmt begin sectional_hsv_-1 png
  gmt set PS_LINE_CAP round
  gmt set PS_LINE_JOIN round
  gmt clip -R-156/-155.5/19/19.5 -Jm-155.75/1:500000 land_polygons_osm_hawaii.bf2 -bi2f -B+g232/252/254
  gmt grdimage hawaii_south.grd -Cfaa_hsv.cpt -Ishading_-1.grd
  gmt grdcontour hawaii_south.grd -Cfaa.cpt -Wthinner,204/217/233
  gmt clip -C
  gmt plot land_polygons_osm_hawaii.bf2 -bi2f -Wthinner,7/72/135 -B+t"Shading -1"
gmt end

gmt begin sectional_hsv_0 png
  gmt set PS_LINE_CAP round
  gmt set PS_LINE_JOIN round
  gmt clip -R-156/-155.5/19/19.5 -Jm-155.75/1:500000 land_polygons_osm_hawaii.bf2 -bi2f -B+g232/252/254
  gmt grdimage hawaii_south.grd -Cfaa_hsv.cpt -Ishading_0.grd
  gmt grdcontour hawaii_south.grd -Cfaa.cpt -Wthinner,204/217/233
  gmt clip -C
  gmt plot land_polygons_osm_hawaii.bf2 -bi2f -Wthinner,7/72/135 -B+t"Shading 0"
gmt end

gmt begin sectional_hsv_1 png
  gmt set PS_LINE_CAP round
  gmt set PS_LINE_JOIN round
  gmt clip -R-156/-155.5/19/19.5 -Jm-155.75/1:500000 land_polygons_osm_hawaii.bf2 -bi2f -B+g232/252/254
  gmt grdimage hawaii_south.grd -Cfaa_hsv.cpt -Ishading_1.grd
  gmt grdcontour hawaii_south.grd -Cfaa.cpt -Wthinner,204/217/233
  gmt clip -C
  gmt plot land_polygons_osm_hawaii.bf2 -bi2f -Wthinner,7/72/135 -B+t"Shading +1"
gmt end

pwessel · March 30, 2021, 5:23pm

For your experiments you may need to set
COLOR_HSV_MAX_S = 0
COLOR_HSV_MIN_V = 0

The defaults are 0.1 and 0.3 partly because printers cannot handle the full gamut, but for digital you may use those 0,0. I do that for animations etc.

gmt math RGB2HSV can be used to test color conversions. For a given h,s,v the shading only messes with s,v leaving h along, then converts back to rgb. So the hue should stay unchanged by any shading. I am not sure how that comparsed to “adding black” and what that does to h,s,v.

pwessel · March 30, 2021, 5:51pm

WIth all those COLOR_HSV_* settings pure 0 or 1, your washed out plot should be all white since all hsv values converted with an intensity of +1 will give white. Likewise, your -1 will give pure black. Again, the default limits on these in GMT just puts some brakes on the illumination effects but leaves hue alone.

KristofKoch · March 31, 2021, 7:21pm

Hi @pwessel, thank you for your explanation. After further reading and experimentation it became clear that the interpolation for negative intensities towards Smin and Vmin leads to a saturation increase for colors with S < 1.

(Image taken from the CookBook)

S is altered instead of only changing V to a lower value. This explains why my light green colors become so saturated with negative intensities. Keeping S constant for negative intensities might be the key to my “I just want to add black”. For positive intensities I found interpolating towards Smax and Vmax quite pleasing.

Is there a way to keep S constant for negative intensities while keeping the default behaviour for positive intensities?

KristofKoch · April 1, 2021, 8:55pm

Ok got some steps closer … setting Smin and Vmin to zero, all negative intensities keep a constant saturation during interpolation. By compressing the range of the lighting intensity with grdmdath (in this case from -0.6 to +0.1) I avoid the very dark colors as well as the default increase in saturation.

Not yet where I want to be but coming closer …

Joaquim · April 6, 2021, 4:21pm

Since I’m adding tons of GDAL kg to GMT Julia I implemented also the shading method pointed out by @Andreas in this post. The GDAL method (from the scarce references found in internet that mostly say use a photoshop type to blend color plus shade) does not change the color space and put illumination in Saturation/Value but blends the color image + shade. This blending can apparently take several forms. Anyway, as mentioned I implemented the equation shown by Andreas.

Ofc, in Julia everything seems insultingly simple (the CPT is the one you provided) and the result is interesting

grdcut("@srtm_relief_03s", R=(-156,-155.5,19,19.5), save="hawaii_south.grd");
I = gdalshade("hawaii_south.grd", color="faa.cpt", zfactor=4);
imshow(I, coast=true, fmt=:png)

KristofKoch · April 6, 2021, 5:52pm

Hi @joaquim, this looks very promising! I’ll give it a try later …

On my quest for beautiful hill shading I stumbled over the work of Leland Brown who developed something he calls “texture shading”.

From his 2014 Texture Shading Presentation:

What Texture Shading Does

Two intuitive descriptions:

A type of high-pass filter, or sharpening filter, on the DEM

Enhances narrow terrain features (details)

Attenuates wide features (major landforms)

A measure of relative height compared to nearby terrain

Light areas are higher than nearby terrain (ridges, peaks)

Dark areas are lower than nearby terrain (canyons, valleys)

Mathematical description:

“Fractional Laplacian” operator

Related to fractals, which accounts for the multi-scale properties

At 200%, becomes an ordinary Laplacian (used in image processing)

DEM (left) and texture shading (right) comparison (Image source)

The usual hill shading (left) and hill shading merged with texture shading (right) (Image source)

I found that his texture shading technique solved quite a few of my problems, for example showing terrain features on the sunlit side and removing too heavily accented smaller details present with the usual hill shading.

Combined with the well known hill shading it makes for some very pleasing maps. As the FAA charts I’m trying to emulate feature hand drawn hill shading I certainly can’t emulate them perfectly. But that texture shading technique comes pretty close to the artwork.

Leland Brown made a lot of information – including source code – available for download. He also published a paper on it: Texture Shading: A New Technique For Depicting Terrain Relief (not included in the above-mentioned downloadable archive).

Maybe it is an idea worth a GMT feature request? Your thoughts?

Joaquim · April 6, 2021, 6:20pm

Yes, that looks very nice but probably not obvious to port the code. And since it seems to use FFTs that means the method cannot be applied to grids with NaNs, at least not without some trickery.
Six years ago we had another discussion like this on different shading methods
http://gmt.soest.hawaii.edu/boards/1/topics/1041

KristofKoch · April 6, 2021, 6:38pm

Hi @Joaquim, multi direction hill shading can be done out of the box with GMT – unfortunately it is not really comparable to texture shading. Tried this one already (144 calls to grdgradient for one plot …).

Andreas · April 6, 2021, 7:30pm

@KristofKoch This looks vey interesting. Do you have a handy on-liner how to apply this on an input grid?

I’ve downloaded the code and should browse through it (“RTFM”) - for a quick test though, it would be super nice. If possible.

KristofKoch · April 6, 2021, 7:59pm

Hi @Andreas, I was unable to build from source (most probably due to lack of experience/knowledge on my part). Mr. Brown made-ready-to-use programs for MacOS and Windows available as well, that’s what I’m using.

It’s not that difficult (after all I got it to work). The User-Guide page 2 has all the steps (with commands to copy&paste) to get you up and running. Way easier than me reproducing it here. (Sorry, you have to RTFM )

Basically

gmt grdcut – get the area of interest
gdal_translate – get it into a digestible format
texture – calculate texture data
texture_image – select contrast and generate image

If you used the supplied demo dataset then the gmt grdcut step can be omitted.

Andreas · April 6, 2021, 9:10pm

@KristofKoch Thank you!

I’m on Linux, so I could not used the precompiled binaries. The instructions given in the User Guide didn’t work for me. Some googling and I found out that you have to include -lm (tell gcc to link with the math library) in the two latter commands:

gcc -O2 -funroll-loops -DNOMAIN -c *.c
gcc -O2 -funroll-loops *.o texture.c -o texture -lm
gcc -O2 -funroll-loops *.o texture_image.c -o texture_image -lm

Then I was able to compile and run texture and texture_image.

The results are very interesting. I’ll be playing more with this.

Thanks again for sharing!

KristofKoch · April 6, 2021, 9:20pm

@Andreas thank you for sharing that -lm info – now it runs on my Linux box as well. No more copying files between machines … awesome!

Joaquim · April 7, 2021, 1:56am

can be collapsed in 1 (GMT way) or 1. in 2. (GDAL way)

KristofKoch · April 7, 2021, 7:14am

@Joaquim would you mind elaborating how the two steps can be collapsed into a pure GMT solution?

Joaquim · April 7, 2021, 12:25pm

See https://docs.generic-mapping-tools.org/dev/cookbook/features.html#grid-file-format
It should be something like grdcut ... -Gfile.flt=gd:AIG though I’m not sure on the format driver name. Ans maybe you don’t even need the :AIG bit since GDAL now recognize formats by extension.

PlanetGus · April 7, 2021, 4:35pm

That would be an awesome addition indeed

Andreas · April 7, 2021, 8:14pm

From the gdal command given in the user guide, gdal_translate -of EHdr [..], it seems to be the EHdr driver (https://gdal.org/drivers/raster/ehdr.html). The AIG driver in gdal unfortunately seems not to support creation of datasets (see https://gdal.org/drivers/raster/index.html).

I’m trying this, but getting bad grids:

gmt grdcut -Gno.nc -R0/10/50/60 @earth_relief_05m

gdal_translate -of EHdr -ot Float32 no.nc no.flt

gmt grdinfo no.flt
no.flt: Title: no.hdr
no.flt: Command:
no.flt: Remark:
no.flt: Gridline node registration used [Cartesian grid]
no.flt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
no.flt: x_min: NaN x_max: NaN x_inc: NaN name: x n_columns: 0
no.flt: y_min: NaN y_max: NaN y_inc: NaN name: y n_rows: 0
no.flt: v_min: NaN v_max: NaN name: z
no.flt: scale_factor: 1 add_offset: 0

Done via gmt:

gmt grdcut -Gno.flt=gd:EHdr -R0/10/50/60 @earth_relief_05m

gmt grdinfo no.flt
no.flt: Title: no.hdr
no.flt: Command:
no.flt: Remark:
no.flt: Gridline node registration used [Cartesian grid]
no.flt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
no.flt: x_min: NaN x_max: NaN x_inc: NaN name: x n_columns: 0
no.flt: y_min: NaN y_max: NaN y_inc: NaN name: y n_rows: 0
no.flt: v_min: NaN v_max: NaN name: z
no.flt: scale_factor: 1 add_offset: 0

Guess it’s a gdal issue?

I’m getting side tracked now but;

Note format:

gmt grdcut -Gno.asc=ei -R0/10/50/60 @earth_relief_05m
grdcut [WARNING]: ESRI Arc/Info ASCII Interchange file must use proxy for NaN; default to -9999

gmt grdinfo no.asc
no.asc: Title:
no.asc: Command:
no.asc: Remark:
no.asc: Pixel node registration used [Cartesian grid]
no.asc: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
no.asc: x_min: 0 x_max: 10 x_inc: 0.0833333333333 name: x n_columns: 120
no.asc: y_min: 50 y_max: 60 y_inc: 0.0833333333333 name: y n_rows: 120
no.asc: v_min: NaN v_max: NaN name: z
no.asc: scale_factor: 1 add_offset: 0

v_min and v_max is both reported as NaN, but if you cat no.asc it contains only non-NaN’s.

Joaquim · April 7, 2021, 10:00pm

Worked fine for me

grdcut @srtm_relief_03s -Ghawaii_south.flt=gd:EHdr -R-156/-155.5/19/19.5

gdalinfo hawaii_south.flt
Driver: EHdr/ESRI .hdr Labelled
Files: hawaii_south.flt
       hawaii_south.flt.aux.xml
       hawaii_south.hdr
       hawaii_south.prj
Size is 601, 601
Coordinate System is:
...

But on a second look, it clipped the negative depths and replaced them by 0!!!

Origin = (-156.000416666666666,19.500416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left  (-156.0004167,  19.5004167) (156d 0' 1.50"W, 19d30' 1.50"N)
Lower Left  (-156.0004167,  18.9995833) (156d 0' 1.50"W, 18d59'58.50"N)
Upper Right (-155.4995833,  19.5004167) (155d29'58.50"W, 19d30' 1.50"N)
Lower Right (-155.4995833,  18.9995833) (155d29'58.50"W, 18d59'58.50"N)
Center      (-155.7500000,  19.2500000) (155d45' 0.00"W, 19d15' 0.00"N)
Band 1 Block=601x1 Type=Float32, ColorInterp=Undefined
  Min=-1.000 Max=4160.000
  Minimum=-1.000, Maximum=4160.000, Mean=1172.830, StdDev=1136.810
  NoData Value=nan
  Metadata:
    STATISTICS_MAXIMUM=4160
    STATISTICS_MEAN=1172.8304905025
    STATISTICS_MINIMUM=-1
    STATISTICS_STDDEV=1136.809938332
    STATISTICS_VALID_PERCENT=100

EDIT: Nope, it’s correct. The remote a SRTM has no negative nodes. We get the same if grdcut writes a nc file.

Andreas · April 8, 2021, 6:37am

Ok, so the file is correct, but gmt is not able to report the ranges(?)

gmt grdcut -Gno.flt=gd:EHdr -R0/10/50/60 @earth_relief_05m

gmt grdinfo no.flt
no.flt: Title: no.hdr
no.flt: Command: 
no.flt: Remark: 
no.flt: Gridline node registration used [Cartesian grid]
no.flt: Grid file format: ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)
no.flt: x_min: NaN x_max: NaN x_inc: NaN name: x n_columns: 0
no.flt: y_min: NaN y_max: NaN y_inc: NaN name: y n_rows: 0
no.flt: v_min: NaN v_max: NaN name: z
no.flt: scale_factor: 1 add_offset: 0

gdalinfo no.flt
Driver: EHdr/ESRI .hdr Labelled
Files: no.flt
       no.flt.aux.xml
       no.hdr
       no.prj
Size is 120, 120
Coordinate System is:
GEOGCRS["GCS_unknown",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["Degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["Degree",0.0174532925199433]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["Degree",0.0174532925199433]]]
Data axis to CRS axis mapping: 1,2
Origin = (0.000000000000000,59.999999999999964)
Pixel Size = (0.083333333333333,-0.083333333333333)
Corner Coordinates:
Upper Left  (   0.0000000,  60.0000000) (  0d 0' 0.00"E, 60d 0' 0.00"N)
Lower Left  (   0.0000000,  50.0000000) (  0d 0' 0.00"E, 50d 0' 0.00"N)
Upper Right (  10.0000000,  60.0000000) ( 10d 0' 0.00"E, 60d 0' 0.00"N)
Lower Right (  10.0000000,  50.0000000) ( 10d 0' 0.00"E, 50d 0' 0.00"N)
Center      (   5.0000000,  55.0000000) (  5d 0' 0.00"E, 55d 0' 0.00"N)
Band 1 Block=120x1 Type=Float32, ColorInterp=Undefined
  Min=-676.500 Max=1497.000 
  Minimum=-676.500, Maximum=1497.000, Mean=20.463, StdDev=217.651
  NoData Value=nan
  Metadata:
    STATISTICS_MAXIMUM=1497
    STATISTICS_MEAN=20.463229166667
    STATISTICS_MINIMUM=-676.5
    STATISTICS_STDDEV=217.65120068201
    STATISTICS_VALID_PERCENT=100