[GRASS-user] r.slope.aspect: Unexpected Results

GRASS GRASS user
1

   I believe that I asked about this before but I don't have the thread
saved, and I don't know where the mail list is archived.

   I have a 10m DEM and tried running r.slope.aspect on it. The output files
(one each for slope and aspect) look like static divided into square cells.
The same problem occurs when I refine the resolution to 1m. I believe that
Markus N. wrote a comment about this but I don't have that message.

   What do I read to learn what the input raster map needs to provide
r.slope.aspect so that the outputs are appropriate?

Thanks,

Rich

2

On Tue, 26 Jan 2010, Rich Shepard wrote:

I have a 10m DEM and tried running r.slope.aspect on it. The output files
(one each for slope and aspect) look like static divided into square cells.
The same problem occurs when I refine the resolution to 1m. I believe that
Markus N. wrote a comment about this but I don't have that message.

   When I run nviz on the same raster map I also get screwey results. I've no
idea what's going on. Here's the results of r.info on that map:

Layer: aberDEM Date: Tue Dec 29 14:05:59 2009
Mapset: beaver_lake Login of Creator: rshepard
Location: Oregon
DataBase: /usr4/grassbase/. Title: Northwest Oregon 10m DEM ( dem )
Timestamp: none

----------------------------------------------------------------------------
Type of Map: raster Number of Categories: 5729
Data Type: CELL
Rows: 7372
Columns: 14816
Total Cells: 109223552
Projection: Lambert Conformal Conic
N: 1334419.15160578 S: 1279151.24118496 Res: 7.49700358
E: 819255.92362222 W: 769192.92822895 Res: 3.37898187
Range of data: min = 7 max = 1476

Data Description:
generated by r.in.gdal

Comments:

r.in.gdal input="/usr4/dem10m-or/hdr.adf" output="dem" title="Northwest Oregon 10m DEM"

Rich

3

Rich Shepard wrote:

   I believe that I asked about this before but I don't have the thread
saved, and I don't know where the mail list is archived.

   I have a 10m DEM and tried running r.slope.aspect on it. The output files
(one each for slope and aspect) look like static divided into square cells.
The same problem occurs when I refine the resolution to 1m. I believe that
Markus N. wrote a comment about this but I don't have that message.

   What do I read to learn what the input raster map needs to provide
r.slope.aspect so that the outputs are appropriate?

First, the current region's grid should at least match that of the
input map, and ideally of the original data.

If the data has been resampled (including projection with r.proj),
this may produce artifacts, particularly if it used method=nearest
(which is the default).

But before you try redoing the projection step, you might want to try
running r.slope.aspect on the "source" data (which may itself be the
result of resampling). If you get similar results, you really need
better source data.

Essentially, slope calculations can be quite sensitive to resampling
artifacts.

Re-projection from lat/lon with equal latitude and longitude
resolutions is problematic as (unless you're close to the equator) a
degree of longitude represents a smaller distance than a degree of
latitude.

Additionally, most of the processes for acquiring geographic data tend
not to produce values on a regular lat/lon grid, so if your data is on
such a grid, it's quite likely to have been resampled already.

Things which may help:

1. Better (i.e. "uncooked") source data. Reprojecting from the
original grid to the desired target grid in one step is better than
reprojecting to lat/lon then again to the target grid.

2. Upsampling the source data with e.g. r.resamp.interp.

3. Filtering the source data with e.g. r.neighbors or r.mfilter.fp.

4. Re-projecting using method=cubic rather than method=nearest.

5. Filtering the projected data with e.g. r.neighbors or r.mfilter.fp.

You wouldn't want to use *all* of these. #1 is nice if you can get it
(and can use it; you need to know the source "projection"), #4 is safe
(it may not be perfect, but it shouldn't be worse than the other
options), but the others could make matters worse if you don't know
what you're doing.

In the worst case, #5 can be used to get rid of the noise at the
expense of resolution (i.e. you'll end up with slope values which
represent an average over a large area).

--
Glynn Clements <glynn@gclements.plus.com>

4

On Tue, 26 Jan 2010, Glynn Clements wrote:

First, the current region's grid should at least match that of the
input map, and ideally of the original data.

Glynn,

   The number of rows and columns has been changed but not the resolution. At
least, not deliberately.

If the data has been resampled (including projection with r.proj), this
may produce artifacts, particularly if it used method=nearest (which is
the default).

   Nope. Source metadata:

Map_Projection:
         Map_Projection_Name: Lambert Conformal Conic
         Lambert_Conformal_Conic:
               Standard_P arallel: 43.000000
               Standard_P arallel: 45.500000
               Longitude_of_Central_Meridian: -120.500000
               Latitude_of_Projection_O rigin: 41.750000
               False_Easting: 1312336.000000
               False_Northing: 0.000000
Planar_Coordinate_Information:
         Planar_Coordinate_Encoding_Method: row and column
         Coordinate_Representation:
               Abscissa_Resolution: 32.828670
               Ordinate_Resolution: 32.828670
         Planar_Distance_Units: User_Defined_Unit
Geodetic_Model:
   Horizontal_Datum_Name: North American Datum of 1983
   Ellipsoid_Name: Geodetic Reference System 80
   Semi-major_Axis: 6378137.000000
   Denominator_of_Flattening_Ratio: 298.257222

   No reprojection necessary.

Essentially, slope calculations can be quite sensitive to resampling
artifacts.

   So I understand.

1. Better (i.e. "uncooked") source data. Reprojecting from the
original grid to the desired target grid in one step is better than
reprojecting to lat/lon then again to the target grid.

   The source originally-imported DEM produces appropriate slope and aspect
maps. The sub-region I cut out with v.in.region (I believe that's the module
I used) doesn't calculate properly. Therefore, the problem is with the
sub-map.

2. Upsampling the source data with e.g. r.resamp.interp.

   I'll try this on the sub-map.

3. Filtering the source data with e.g. r.neighbors or r.mfilter.fp.

   Source data are OK.

4. Re-projecting using method=cubic rather than method=nearest.

   Data came as LCC; no reprojecting needed.

You wouldn't want to use *all* of these. #1 is nice if you can get it
(and can use it; you need to know the source "projection"),

   See above.

Rich

5

On Tue, 26 Jan 2010, Glynn Clements wrote:

2. Upsampling the source data with e.g. r.resamp.interp.

Glynn,

   When would one use r.resamp.interp and when r.resamp.rst? The latter
incorporates functions of r.slope.aspect with the resampling/resolution
change, but what decision criteria guide which module to use?

Thanks,

Rich

6

Rich Shepard wrote:

   The source originally-imported DEM produces appropriate slope and aspect
maps. The sub-region I cut out with v.in.region (I believe that's the module
I used) doesn't calculate properly. Therefore, the problem is with the
sub-map.

It's possible that the cropping may have changed the grid slightly.

What are the grid parameters (bounds, resolution, rows/cols) for the
source map and the cropped version?

You could try repeating the cropping process, using g.region's align=
option to align the grid to the source map.

--
Glynn Clements <glynn@gclements.plus.com>

7

Rich wrote:

Source metadata:

Map_Projection:
        Map_Projection_Name: Lambert Conformal Conic
        Lambert_Conformal_Conic:
              Standard_Parallel: 43.000000
              Standard_Parallel: 45.500000
             
Longitude_of_Central_Meridian: -120.500000
Latitude_of_Projection_Origin: 41.750000
False_Easting: 1312336.000000
False_Northing: 0.000000
Planar_Coordinate_Information:
       
Planar_Coordinate_Encoding_Method: row and column
        Coordinate_Representation:
             
Abscissa_Resolution: 32.828670
Ordinate_Resolution: 32.828670
        Planar_Distance_Units:
User_Defined_Unit
Geodetic_Model:
    Horizontal_Datum_Name: North American Datum of 1983
    Ellipsoid_Name: Geodetic Reference System 80
    Semi-major_Axis: 6378137.000000
    Denominator_of_Flattening_Ratio: 298.257222

  The source originally-imported DEM produces appropriate slope and
aspect maps. The sub-region I cut out with v.in.region (I believe
that's the module I used) doesn't calculate properly. Therefore, the
problem is with the sub-map.

after you do the zoom run:
  g.region -p
and see if the resolution has been modified (bounds take precedence by
default). You can use the -a flag to make the resolution have precedence
at the expense of maintaining the exact bounds:
  g.region -a res=

with whatever the original map resolution was.

* take care if the original map's bounds are not exactly divisible by
the resolution. (e.g. n=125 s=75 res=50 rows=1) then -a does not work
as you might expect.

Hamish

8

On Tue, 26 Jan 2010, Glynn Clements wrote:

It's possible that the cropping may have changed the grid slightly.

Glynn,

   I do hope it's something simple. I can't do more until I have a good DEM
map for the project area.

What are the grid parameters (bounds, resolution, rows/cols) for the
source map and the cropped version?

Source map (demNW raster map):

GRASS 6.4.0svn (Oregon):/usr4/grassbase > g.region rast=demNW -p
projection: 99 (Lambert Conformal Conic)
zone: 0
datum: nad83
ellipsoid: grs80
north: 1735231.43724681
south: 1120657.1324816
west: 143667.18968253
east: 952979.22944945
nsres: 32.82281055
ewres: 32.82281055
rows: 18724
cols: 24657
cells: 461677668

   Resolution in feet.

Cropped map (abernethy vector boundary map):

GRASS 6.4.0svn (Oregon):/usr4/grassbase > g.region vect=abernethy -p
projection: 99 (Lambert Conformal Conic)
zone: 0
datum: nad83
ellipsoid: grs80
north: 1334419.25160578
south: 1279151.24118496
west: 769192.9282895
east: 819255.92362222
nsres: 4.99982001
ewres: 4.99979979
rows: 11054
cols: 10013
cells: 110683702

   Resolution in metres.

You could try repeating the cropping process, using g.region's align=
option to align the grid to the source map.

   I'll try this.

Thanks,

Rich

9

On Tue, 26 Jan 2010, Glynn Clements wrote:

You could try repeating the cropping process, using g.region's align=
option to align the grid to the source map.

   Now I cannot find my notes on how I cropped the larger DEM to the boundary
of the project area. I know that v.overlay works on two vector maps but ...
I've dropped the ball on that operation.

Rich

10

On Tue, 26 Jan 2010, Rich Shepard wrote:

Now I cannot find my notes on how I cropped the larger DEM to the
boundary of the project area. I know that v.overlay works on two vector
maps but ... I've dropped the ball on that operation.

   Found it!

   r.mapcalc "newmap=oldmap"

Rich

11

On Tue, 26 Jan 2010, Glynn Clements wrote:

You could try repeating the cropping process, using g.region's align=
option to align the grid to the source map.

Glynn,

   FIXED the problem.

   I re-ran 'r.mapcalc "aberDEM=demNW"' and it cleaned up all the problems.
r.slope.aspect produces output maps as it should. No more static in
vertical/horizontal striped squares.

   No idea what went wrong the first time, but it's all fixed now.

   Thank you very much for this suggestion. From now on, any problems with
converted data of any type and the first thing I'll do is recreate the map
giving problems.

Rich

12

Rich Shepard wrote:

> It's possible that the cropping may have changed the grid slightly.

Glynn,

   I do hope it's something simple. I can't do more until I have a good DEM
map for the project area.

> What are the grid parameters (bounds, resolution, rows/cols) for the
> source map and the cropped version?

Source map (demNW raster map):

nsres: 32.82281055
ewres: 32.82281055

   Resolution in feet.

Cropped map (abernethy vector boundary map):

nsres: 4.99982001
ewres: 4.99979979

   Resolution in metres.

32.82281055 international feet is 10.00439265564 metres.

32.82281055 US survey feet is 10.00441266446533 metres.

--
Glynn Clements <glynn@gclements.plus.com>

13

Glynn Clements wrote:

Rich Shepard wrote:

Source map (demNW raster map):
    
nsres: 32.82281055
ewres: 32.82281055
    
   Resolution in feet.

Cropped map (abernethy vector boundary map):
    
nsres: 4.99982001
ewres: 4.99979979
    
   Resolution in metres.
    
32.82281055 international feet is 10.00439265564 metres.

32.82281055 US survey feet is 10.00441266446533 metres.
  

I don't think it is possible to have different units within the same location because units are set globally in PERMANENT/PROJ_UNITS. Is it?

14

Markus Metz wrote:

I don't think it is possible to have different units within the same
location because units are set globally in PERMANENT/PROJ_UNITS. Is it?

Correct.

Apparently the issue has been reolved. Which is fortunate; as the
previous post doesn't actually make much sense.

--
Glynn Clements <glynn@gclements.plus.com>