[GRASS-dev] [GRASS GIS] #3860: GRASS GIS producing different slope than GDAL

[GRASS GIS]

#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>