[GRASS-dev] Re: [gdal-dev] gdaldem hillshade directions / GRASS r.shaded.relief

GRASS GIS GRASS GIS developer
,
1

Vincent,

I'm CC'ing the grass-dev list as gdaldem shares the same formula with the
GRASS r.shaded.relief utility and I also think they are affected (I've only
compared the code, not tested r.shaded.relief, so I could be wrong of course)

I think your analysis is right. To check, I've created an artificial DEM that
is a pyramid with 4 faces oriented to north, east, south and west. Indeed the
azimuth parameter isn't interpreted as documented but as you've noticed. The
documentation is what I would expect 0 = north, 90 = east, 180 = south and
270 = west.

The fix is simple. Instead of the following code in
http://trac.osgeo.org/grass/browser/grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py

#correct azimuth to East (GRASS convention):
# this seems to be backwards, but in fact it works so leave it.
az = float(azimuth) - 90

I'd suggest :

az = 180 - float(azimuth)

The interesting thing is that the issue probably got unnoticed since most
users will use azimuth = 315 (north-west) according to the best practice. But
as 315 - 90 = 225, 180 - 315 = -135 and 225 = - 135 + 360..., the usual value
is the only one for which both formulas are equivalent...

What do the GRASS developers think ?

Best regards,

Even

PS:

Python script to generate the pyramid :

import gdal
import osr

ds = gdal.GetDriverByName('GTiff').Create('testdem.tif', 100, 100, 1)
ds.SetGeoTransform([2,0.01,0,49,0,-0.01])
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326)
ds.SetProjection(sr.ExportToWkt())
for j in range(100):
    data = ''
    for i in range(100):
        val = 255 - 5 * max(abs(50-i),abs(50-j))
        data = data + ('%c' % (val))
    ds.GetRasterBand(1).WriteRaster(0,j,100,1,data)

ds = None

gdaldem commandline :

gdaldem hillshade testdem.tif testdem_shaded.tif -s 111120 -z 100 -az XXXX

Le Tuesday 18 May 2010 14:11:09 Vincent Schut, vous avez écrit :

Is it just me, or is the result/documentation for gdaldem hillshade off
by 90 degrees?

While the docs state "azimuth of the light, in degrees. 0 if it comes
from the top of the raster, 90 from the east, ..." (so: clockwise from
north), I find that a sun azimuth of 0 degrees is west, 90 is south, 180
is east and 270 is north (counter-clockwise from west)...

Vincent.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
gdal-dev Info Page

2

Even wrote:

#correct azimuth to East (GRASS convention):
# this seems to be backwards, but in fact it works so
leave it.
az = float(azimuth) - 90

I'm not into checking-in fudge factors which I don't understand,
and I seem to recall that when I wrote that comment I did
actually investigate the r.mapcalc algorithm and determine that
it was "in fact" doing the right thing, or at least verified the
output was correct.

GRASS has, for decades, used the mathematical convention of
theta measured CCW from the positive x-axis for measuring
aspect, not the nautical convention of of degrees measured CW
from north. Both ways are neither right nor wrong, it's just the
convention used. (a reminder that grass's roots are in modeling
and mathematical analysis, not cartography)

The typical way to convert between the two conventions is:
  naut = 90 - theta
  theta = 90 - naut

There is the possiblity that when gdal copied this code from
GRASS they only got/rewrote half of it, or lost the matching
color table, etc. ??

If you like add it as a wish in the grass trac'er assigned to me,
I could dig into it at some time in the future. sorry I don't
have much time at all right now to look into itfurther. Only to
say that I have looked into this in the past and signed off on
it as ok back then (for use within GRASS, no idea about gdal's
implimentation).

hope that helps a little,
Hamish

3

Hamish,

I've taken some time to experiment a bit with GRASS and GRASS behaviour is
correct of course, so sorry for the noise.

I've finally identified why gdaldem (and initially Matt's hillshade utility)
got it wrong. GRASS r.mapcalc atan(x,y) function has its parameters reversed
in comparison to C atan2() used by GDAL : atan(x,y)=atan2(y,x)....

So the aspect computation was wrong which could also be seen as a wrong
interpretation of the azimuth parameter. But for the particular value of
az=315deg, the beauty of trigonometry makes it such that cos(az-pi/2-atan2
(y,x)) = cos(az-pi/2- atan2(x,y))...

Fixed as ticket #3586 (gdaldem hillshade azimuth parameter wrongly interpreted) – GDAL

Best regards,

Even

Le Wednesday 19 May 2010 03:57:38 Hamish, vous avez écrit :

Even wrote:
> #correct azimuth to East (GRASS convention):
> # this seems to be backwards, but in fact it works so
> leave it.
> az = float(azimuth) - 90

I'm not into checking-in fudge factors which I don't understand,
and I seem to recall that when I wrote that comment I did
actually investigate the r.mapcalc algorithm and determine that
it was "in fact" doing the right thing, or at least verified the
output was correct.

GRASS has, for decades, used the mathematical convention of
theta measured CCW from the positive x-axis for measuring
aspect, not the nautical convention of of degrees measured CW
from north. Both ways are neither right nor wrong, it's just the
convention used. (a reminder that grass's roots are in modeling
and mathematical analysis, not cartography)

The typical way to convert between the two conventions is:
  naut = 90 - theta
  theta = 90 - naut

There is the possiblity that when gdal copied this code from
GRASS they only got/rewrote half of it, or lost the matching
color table, etc. ??

If you like add it as a wish in the grass trac'er assigned to me,
I could dig into it at some time in the future. sorry I don't
have much time at all right now to look into itfurther. Only to
say that I have looked into this in the past and signed off on
it as ok back then (for use within GRASS, no idea about gdal's
implimentation).

hope that helps a little,
Hamish

4

Even wrote:

I've taken some time to experiment a bit with GRASS and
GRASS behaviour is correct of course, so sorry for the noise.

I've finally identified why gdaldem (and initially Matt's hillshade
utility) got it wrong. GRASS r.mapcalc atan(x,y) function has its
parameters reversed in comparison to C atan2() used by GDAL :
atan(x,y)=atan2(y,x)....

So the aspect computation was wrong which could also be seen as a wrong
interpretation of the azimuth parameter. But for the particular value of
az=315deg, the beauty of trigonometry makes it such that
cos(az-pi/2-atan2 (y,x)) = cos(az-pi/2- atan2(x,y))...

Fixed as ticket http://trac.osgeo.org/gdal/ticket/3586

well spotted! Let's just call this one a result of C's math lib being
weird and then none of us have to take the blame. :slight_smile:

[so here's finally a case where lat,long misorder could help?]

cheers,
Hamish