#3860: GRASS GIS producing different slope than GDAL
------------------------+-------------------------
Reporter: mazingaro | Owner: grass-dev@…
Type: defect | Status: reopened
Priority: normal | Milestone: 7.6.2
Component: Raster | Version: unspecified
Resolution: | Keywords: slope
CPU: x86-64 | Platform: Linux
------------------------+-------------------------
Comment (by mankoff):
There may be a documentation bug. There may also be a real bug. According
to my understanding of the manual, slope should not be impacted by region
unless the `-a` flag is set:
To ensure that the raster elevation map is not inappropriately resampled,
the settings for the current region are modified slightly (for the
execution of the program only): the resolution is set to match the
resolution of the elevation raster map and the edges of the region (i.e.
the north, south, east and west) are shifted, if necessary, to line up
along edges of the nearest cells in the elevation map. If the user really
wants the raster elevation map resampled to the current region resolution,
the -a flag should be specified.
Given that, slope calculated in different regions should be the same - if
it is reported correctly from as per the comment above.
{{{
rm -fR nc_spm_08_grass7/DEBUG_SLOPE
grass -c ./nc_spm_08_grass7/DEBUG_SLOPE
export GRASS_OVERWRITE=1
g.region raster=elevation
r.slope.aspect elevation=elevation slope=slope
r.slope.aspect -a elevation=elevation slope=slope_a
g.region res=30 -pa
r.slope.aspect elevation=elevation slope=slope_30
r.slope.aspect -a elevation=elevation slope=slope_30_a
g.region res=100 -pa
r.slope.aspect elevation=elevation slope=slope_100
r.slope.aspect -a elevation=elevation slope=slope_100_a
g.region raster=elevation
r.out.gdal input=elevation output=elevation.tif --o
gdaldem slope elevation.tif slope.tif -of GTiff # -b 1 -s 1.0 not needed
r.in.gdal input=slope.tif output=gdal
g.region raster=elevation # reset region for r.univar
for r in $(g.list mapset=. type=raster); do
(echo -n "${r} "; r.univar ${r} | grep range)
done | sort | tac
diff <(r.univar slope) <(r.univar gdal)
}}}
Has these results:
{{{
slope range: 38.6894
slope_a range: 38.6894
slope_30 range: 14.9465
slope_30_a range: 25.3968
slope_100 range: 4.57874
slope_100_a range: 8.92367
gdal range: 38.6894
15c15
< sum: 7803645.55388512
---
> sum: 7803645.58213263
}}}
As expected `r.aspect.slope` and `r.aspect.slope -a` behave the same when
the region is equal to the raster region. But I would expect `slope_30`
and `slope_100` to both match `slope`, based on the documentation.
The gdal calculated slope does match the GRASS calculated slope.
Summary:
+ The original issue can be solved by setting the region to the raster:
`g.region raster=DEM -a`
+ There appears to be a code, documentation, or my-understanding bug about
the behavior of `r.slope.aspect`. Certainly the used of the word
'slightly' is poorly defined here. If the region resolution and edges are
changed, how more un-slight could it be?
--
Ticket URL: <https://trac.osgeo.org/grass/ticket/3860#comment:14>
GRASS GIS <https://grass.osgeo.org>
