1

I'm trying to generate a slope raster from a DEM raster. It seems there is a pretty significant error due to the 'z-factor' option.

I'm fairly new to this, but some other posts and responses indicate a general workflow:

1) obtain a DEM. The ones I got (NED) are unprojected: horizontal coords are lat-lon, and vertical coords are meters.

Two possible branches here, A or B:

A2) project the DEM to UTM so that horizontal and vertical coords are in the same units A3) then do slope analysis on the UTM-projected DEM, using z-factor = 1.0

OR

B2) do slope analysis on the unprojected DEM, using z-factor = 111120 (which has been suggested in some posts as a number of meters per degree - see error discussion below) B3) then project the raster to UTM

Problem with option A: the projection process introduces gridded artifacts, and these get amplified in the slope raster:

Original DEM (NAD83, unprojected; zoomed in to one peak; min/max set to show high contrast here):

Original DEM (unprojected)

DEM projected to NAD83 UTM zone 10, note gridded artifacts:

enter image description here

slope analysis output from NAD83-UTMz10-projected DEM, note the same artifacts:

enter image description here

Problem with option B: 111120 seems vastly oversimplified and produces significant errors, see error table below - though, gridded artifacts are greatly reduced since they don't exist in the DEM that you run slope analysis on. First, number of meters per degree of latitude is not the same at all latitudes, due to the earth not being a perfect sphere. Second, and more obvious, number of meters per degree of longitude is a function of latitude, varying from ~111km at the equator to 0 at the poles.

For our area (Nevada County, California) here are some numbers from an online calculator, referenced in other places as a trusted source (but feel free to verify these):

Nevada County S extent ~ 39.0deg (WGS84) Nevada County N extent ~ 39.5deg (WGS84) Nevada City ~ 39.25deg (WGS84)

39.0 deg: 1 deg lat = 111015.45 meters 1 deg lon = 86626.37 meters

39.25 deg: 1 deg lat = 111020.23 meters 1 deg lon = 86320.71 meters

39.5 deg: 1 deg lat = 111025.01 meters (about 9.5 meters longer than at S extent of Nevada County) 1 deg lon = 86013.39 meters (about 613 meters shorter than at S extent of Nevada County)

(at equator, for reference:) 0.0 deg: 1 deg lat = 110574.27 meters 1 deg lon = 111319.46 meters

So, making some straightforward assumptions about how qgis uses the 'Z Factor' argument, here are some resulting errors (note this is in Nevada County but not the same location as the above images), along with Range-Bearing Line from Terrain Navigator Pro which is treated as gospel in this case; also note that this is just for one sample set of from/to; error would vary by compass direction of the from-to line also:

enter image description here

(z-factor is entered by hand for each row; all the other columns are calculated from it)

So, any ideas? Am I calculating something wrong, or is there a better workflow, or is there just no really good way to get around the error?

Tom Grundy
  • 925
  • 1
  • 6
  • 24
  • Look at this answer: http://gis.stackexchange.com/questions/79730/grid-like-lines-appearing-after-using-curvature-tool/79736#79736 You should see an option for bilinear resample under Raster-->Projections-->Warp – Scro Dec 30 '13 at 18:15
  • Your problem in (A) appears to have been diagnosed correctly by @Scro. Your questions related to (B) are answered at http://gis.stackexchange.com/questions/14750/global-dem-to-slope-calculation/40464#40464. – whuber Dec 30 '13 at 20:44
  • it seems to me that there is more than a nearest neighbour resampling issue in your case. Have you tried the reprojection with another software (e.g. gdalwarp)? – radouxju Dec 30 '13 at 20:51
  • @radouxju - gdalwarp is exactly the utility that's used for reprojections in QGIS. – Scro Dec 30 '13 at 21:56
  • The method I used to get the UTM projection in the picture was the 'trivial' projection: right-click the layer in the table of contents, and hit 'Save As...' which lets you pick a CRS. So maybe that uses a cruddy algorithm? I'll try Raster -> Projections -> Warp with the bilinear resample. – Tom Grundy Dec 31 '13 at 00:23
  • I'm trying this (the command generated in the GUI, kindof - the to and from SRS fields don't take correctly so I did them by hand - 4269=NAD83, 29610=NAD83 UTM zone10: gdalwarp -s_srs EPSG:4269 -t_srs EPSG:26910 -r bilinear -of GTiff "D:\GIS layers\NED\40122\float\floatn40w122_13.flt" "D:/GIS layers/NED/projected/floatn40w122_13_utm10n.tif" – Tom Grundy Dec 31 '13 at 15:37
  • The above seems to work nicely. Thanks for the info. – Tom Grundy Dec 31 '13 at 15:55

0 Answers0