GRASS GIS producing different slope than GDAL

[mazingaro]

2

When producing a slope in GRASS GIS 7.4 and 7.6 the results are different than with QGIS (GDAL) and ArcMap for the same input data and same parameter (Z ratio is 1.0 and degrees). QGIS and ArcMap generate almost the same output. In GRASS, I tried both for DCELL and FCELL and the results are same. The original layer has a resolution on 1 m and the GRASS region resolution is set to 0.5 m.

Code GRASS: r.slope.aspect elevation=DMR@PERMANENT slope=slope

Code GDAL: gdaldem slope .../DMR.tif .../slope.tif -of GTiff -b 1 -s 1.0

Do I do something wrong or this can be a bug? Left - GRASS calculated slope and right - GDAL calculated one

[mlennert]

What does GRASS GIS give you when you work at the resolution of the data, i.e. 1m, instead of 0.5m ?

[mazingaro]

Replying to mlennert:

What does GRASS GIS give you when you work at the resolution of the data, i.e. 1m, instead of 0.5m ?

The result is the same (with the -a flag)

[mlennert]

Replying to mazingaro:

Replying to mlennert:

What does GRASS GIS give you when you work at the resolution of the data, i.e. 1m, instead of 0.5m ?

The result is the same (with the -a flag)

Weird. Are you sure you are working on the same region extent ?

I just tried with the NC demo data elevation map, setting the region with g.region rast=elevation in GRASS. Here are the univar stats of the difference map between the two:

> r.univar diff
 100%
total null and non-null cells: 8100000
total null cells: 22784

Of the non-null cells:
----------------------
n: 8077216
minimum: -9.67979e-05
maximum: 9.72748e-05
range: 0.000194073
mean: -1.50386e-09
mean of absolute values: 1.30826e-05
standard deviation: 1.69973e-05
variance: 2.88908e-10
variation coefficient: -1.13025e+06 %
sum: -0.0121469677105779

IOW: very small (IMHO negligeable) differences.

However, when I calculate the slope in GRASS with the same region extent, but a resolution set to 5m, I get much more significant differences (diff and univar stats calculated at 5m resolution):

> r.univar diff_restest
 100%
total null and non-null cells: 2025000
total null cells: 5696

Of the non-null cells:
----------------------
n: 2019304
minimum: 0
maximum: 19.4712
range: 19.4712
mean: 3.74369
mean of absolute values: 3.74369
standard deviation: 2.68803
variance: 7.22549
variation coefficient: 71.8015 %
sum: 7559653.93641658

[mazingaro]

Replying to mlennert:

Replying to mazingaro:

Replying to mlennert:

What does GRASS GIS give you when you work at the resolution of the data, i.e. 1m, instead of 0.5m ?

The result is the same (with the -a flag)

Weird. Are you sure you are working on the same region extent ?

I just tried with the NC demo data elevation map, setting the region with g.region rast=elevation in GRASS. Here are the univar stats of the difference map between the two:

> r.univar diff
 100%
total null and non-null cells: 8100000
total null cells: 22784

Of the non-null cells:
----------------------
n: 8077216
minimum: -9.67979e-05
maximum: 9.72748e-05
range: 0.000194073
mean: -1.50386e-09
mean of absolute values: 1.30826e-05
standard deviation: 1.69973e-05
variance: 2.88908e-10
variation coefficient: -1.13025e+06 %
sum: -0.0121469677105779

IOW: very small (IMHO negligeable) differences.

However, when I calculate the slope in GRASS with the same region extent, but a resolution set to 5m, I get much more significant differences (diff and univar stats calculated at 5m resolution):

> r.univar diff_restest
 100%
total null and non-null cells: 2025000
total null cells: 5696

Of the non-null cells:
----------------------
n: 2019304
minimum: 0
maximum: 19.4712
range: 19.4712
mean: 3.74369
mean of absolute values: 3.74369
standard deviation: 2.68803
variance: 7.22549
variation coefficient: 71.8015 %
sum: 7559653.93641658

This is a normal difference, I don't see problems here. The univar stats of the defference map were made for GRASS and GDAL produces slopes?

[mlennert]

The first difference is between gdaldem and r.slope.aspect results at the resolution of the input data.

The second difference is between two runs of r.slope.aspect, one at the resolution of the input data (10m), one at half that resolution (5m). It is meant to show that differences because of different region settings are probably much more significant than differences between the two programs.

I personally would not consider an average difference of almost 4 degrees "normal" if I expected same results...

[mazingaro]

Replying to mlennert:

The first difference is between gdaldem and r.slope.aspect results at the resolution of the input data.

The second difference is between two runs of r.slope.aspect, one at the resolution of the input data (10m), one at half that resolution (5m). It is meant to show that differences because of different region settings are probably much more significant than differences between the two programs.

I personally would not consider an average difference of almost 4 degrees "normal" if I expected same results...

Yes, you are right. Anyway, that means that your slopes are not so different (GRASS-GDAL) while i got them way more different (+10 degrees). Which GRASS are you using?

[mlennert]

Replying to mazingaro:

Replying to mlennert:

The first difference is between gdaldem and r.slope.aspect results at the resolution of the input data.

The second difference is between two runs of r.slope.aspect, one at the resolution of the input data (10m), one at half that resolution (5m). It is meant to show that differences because of different region settings are probably much more significant than differences between the two programs.

I personally would not consider an average difference of almost 4 degrees "normal" if I expected same results...

Yes, you are right. Anyway, that means that your slopes are not so different (GRASS-GDAL) while i got them way more different (+10 degrees). Which GRASS are you using?

Current stable, i.e. 7.6.

Could you make your data available for testing ?

[mazingaro]

Replying to mlennert:

Current stable, i.e. 7.6.

Could you make your data available for testing ?

Dear mlennert, sorry for replying just now - was a wild day at work. Here is the data: ​https://wetransfer.com/downloads/6a96559edb293d40a220c8c65a68e44520190605132711/143eb1 Thanks!

[mankoff]

I think the issue is related to the region.

The r.slope.aspect documentation in the NOTES section implies that the region is adjusted to the raster. I do not find that to be the case. Specifically in the NC test data, I've created the following slope rasters, and then looked at their univariate statistics:

g.region raster=elevation
r.slope.aspect elevation=elevation slope=slope_0
r.slope.aspect -a elevation=elevation slope=slope_1
g.region res=30 -pa
r.slope.aspect elevation=elevation slope=slope_2
r.slope.aspect -a elevation=elevation slope=slope_3

for i in $(seq 0 3); do
  (echo -n "slope_${i} "; r.univar slope_${i} | grep range) | tr '\n' ' '
done

Results are:

slope_0 range: 36.3347 
slope_1 range: 36.3347 
slope_2 range: 13.7754 
slope_3 range: 25.3968

[mazingaro]

Replying to mankoff:

I think the issue is related to the region.

The r.slope.aspect documentation in the NOTES section implies that the region is adjusted to the raster. I do not find that to be the case. Specifically in the NC test data, I've created the following slope rasters, and then looked at their univariate statistics:

g.region raster=elevation
r.slope.aspect elevation=elevation slope=slope_0
r.slope.aspect -a elevation=elevation slope=slope_1
g.region res=30 -pa
r.slope.aspect elevation=elevation slope=slope_2
r.slope.aspect -a elevation=elevation slope=slope_3

for i in $(seq 0 3); do
  (echo -n "slope_${i} "; r.univar slope_${i} | grep range) | tr '\n' ' '
done

Results are:

slope_0 range: 36.3347 
slope_1 range: 36.3347 
slope_2 range: 13.7754 
slope_3 range: 25.3968

Ok, thank you. I am going to close the ticket then.

[neteler]

Replying to mankoff:

I think the issue is related to the region.

The r.slope.aspect documentation in the NOTES section implies that the region is adjusted to the raster. I do not find that to be the case.

Reopening for manual to be corrected (if that's the solution to the reported differences).

[mmetz]

Replying to mankoff:

I think the issue is related to the region.

The r.slope.aspect documentation in the NOTES section implies that the region is adjusted to the raster. I do not find that to be the case. Specifically in the NC test data, I've created the following slope rasters, and then looked at their univariate statistics:

g.region raster=elevation
r.slope.aspect elevation=elevation slope=slope_0
r.slope.aspect -a elevation=elevation slope=slope_1
g.region res=30 -pa
r.slope.aspect elevation=elevation slope=slope_2
r.slope.aspect -a elevation=elevation slope=slope_3

for i in $(seq 0 3); do
  (echo -n "slope_${i} "; r.univar slope_${i} | grep range) | tr '\n' ' '
done

Please also check the extents and resolution of slope_0, slope_1, slope_2, slope_3. You need to set the region to the raster map before running r.univar because r.univar uses the current region. Therefore these results are not what you intend to get because they are all obtained with the same region settings:

Results are:

slope_0 range: 36.3347 
slope_1 range: 36.3347 
slope_2 range: 13.7754 
slope_3 range: 25.3968

Alternatively, try r.info -s which ignores the current region and reports simple raster stats stored in metadata.

Replying to neteler:

Replying to mankoff:

I think the issue is related to the region.

The r.slope.aspect documentation in the NOTES section implies that the region is adjusted to the raster. I do not find that to be the case.

Reopening for manual to be corrected (if that's the solution to the reported differences).

I think the manual is correct.

[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?

[mankoff]

Running r.univar at the region resolution (as above) or at the raster resolution:

for r in $(g.list mapset=. type=raster); do
  g.region raster=${r} -a
  (echo -n "${r} "; r.univar ${r} | grep range)
done | sort | tac

Does not change anything.

[mmetz]

Replying to mankoff:

Running r.univar at the region resolution (as above) or at the raster resolution:

for r in $(g.list mapset=. type=raster); do
  g.region raster=${r} -a
  (echo -n "${r} "; r.univar ${r} | grep range)
done | sort | tac

Does not change anything.

Right, this was a bug in r.slope.aspect, alignment to the input raster did not work (difference between G_get_window() and Rast_get_window()). Fixed in master, relbr76, and relbr74.