[gdal-dev] gdaldem hillshade directions / GRASS r.shaded.relief
Even Rouault
even.rouault at mines-paris.org
Tue May 18 15:06:06 EDT 2010
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
More information about the gdal-dev
mailing list