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

[Even Rouault]

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 at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev