Does the slope tool's geodesic method have a bug?

[GarrettBradford]

I am working on creating slope surfaces for the entire lower 48 using 1/3 arc-second data from USGS. Recently, I have begun experimenting with running the slope tool using the geodesic method without re-projecting the original rasters (which are in NAD83) to a local coordinate system.

For the most part, the geodesic method appears to work, but it introduces some "noise" in the resulting slope surface when compared slope derived from a planar method and local coordinate system. The major issue I see is that it produces illogical and clearly incorrect values along the 90th west meridian (see maps attached, I apologizes in advance for the ugly graphic). I thought this was user error, or an issue with the input DEM, however I have replicated this using various tiles and input data sources. I have not completed the lower 48, however, everywhere else in the lower 48 where I have completed this appears fine.

What is going on here?

[DanPatterson_Retired]

NAD83 is a datum and not a coordinate system, so can we assume you are working with a Geographic Coordinate system (ie GCS NAD83) with coordinates in decimal degrees?  And does that coordinate system define the z-units?

How Slope works—Help | ArcGIS Desktop 

So the coordinate system isn't critical except for X, Y and Z  being defined

Geodesic method

The geodesic method measures slope in a geocentric 3D coordinate system—also called the Earth Centered, Earth Fixed (ECEF) coordinate system—by considering the shape of the earth as an ellipsoid. The computation result will not be affected by how the dataset is projected. It will use the z-units of the input raster if they are defined in the spatial reference. If the spatial reference of the input does not define the z-units, you will need to do so with the z-unit parameter. The geodesic method produces a more accurate slope than the planar method.

[GarrettBradford]

Hi Dan, thanks for the reply!

Spatial Reference for 1/3 arc-second data: NAD83: EPSG Projection -- Spatial Reference . Sorry if that wasn't clear.

Defining Z-units on the input is not required (see the documentation you linked - you can select them when running the slope tool), however, I did try that as a possible fix and I get the same result.

I have read the tool documentation in some detail and have not found any answers there. As I said, the output appears reasonable everywhere besides the 90th west meridian. Can you think of any reason the slope of the fitted topographic surface relative to the ellipsoid would generate the types of values I have shown in the attachment?

[DanPatterson_Retired]

If it just at that location, and you used the exact same procedures, I would contact the data provider.  They are the most likely group to provide commentary (ie dem cell size difference?)

[GarrettBradford]

Hi Dan, this is not a data issue and it is not because of cell size difference. I believe it is due to a problem with the geodesic method as programed in the slope tool. To confirm this you can create a constant raster (which should by definition has a slope of 0 everywhere).

The code below produces two example using both SRID 4269 ("NAD83") and 26915 ("UTM"). In both cases the error is apparent when slope is calculated using the geodesic method.

import arcpyimport arcpy.sa as SA
#NAD83 example
nad83=arcpy.SpatialReference(4269)
arcpy.env.outputCoordinateSystem = nad83
const_ras = SA.CreateConstantRaster(1.0,"FLOAT",0.0001,"-90.1 29.9 -89.9 30.1")
slope_ras = SA.Slope(const_ras,"DEGREE",'',"GEODESIC","METER")

#UTM example
utm15n = arcpy.SpatialReference(26915)
arcpy.env.outputCoordinateSystem = utm15n
const_ras = SA.CreateConstantRaster(1.0,"FLOAT",10,"779470 3311230 799367 3333922")
slope_ras = SA.Slope(const_ras,"DEGREE",'',"GEODESIC","METER")
slope_ras2 = SA.Slope(const_ras,"DEGREE",'',"PLANAR","METER") #This shows what the slope should look like‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

something is not right...

Look at the median for geodesic... I don't get it, unless the constant raster shouldn't be a flat plane but a curved one conforming to the surface????

And the max... must be some edge thing?

If you are in the U.S, you should file it.  If I report it, Canada has to confirm it, then it goes to Redlands (might be retired by then

 import arcpy
import arcpy.sa as SA
#NAD83 example
nad83=arcpy.SpatialReference(4269)
arcpy.env.outputCoordinateSystem = nad83
const_ras = SA.CreateConstantRaster(1.0,"FLOAT",0.0001,"-90.1 29.9 -89.9 30.1")
slope_ras = SA.Slope(const_ras,"DEGREE",'',"GEODESIC","METER")

#UTM example
utm15n = arcpy.SpatialReference(26915)
arcpy.env.outputCoordinateSystem = utm15n
const_ras = SA.CreateConstantRaster(1.0,"FLOAT",10,"779470 3311230 799367 3333922")
slope_ras = SA.Slope(const_ras,"DEGREE",'',"GEODESIC","METER")
slope_ras2 = SA.Slope(const_ras,"DEGREE",'',"PLANAR","METER") #This shows what the slope should look like



a_planar = arcpy.RasterToNumPyArray(slope_ras2)

a_geod = arcpy.RasterToNumPyArray(slope_ras)

a_planar.max(), a_planar.min()
(0.0, -3.4028235e+38)

a_geod.max(), a_geod.min()
(89.99978, -3.4028235e+38)

np.median(a_planar), np.median(a_geod)
(0.0, 1.2710605)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

even clipping is showing some weirdness in the geodetic

a_p_small = a_planar[1000:1500, 1000:1500]

a_g_small = a_geod[1000:1500, 1000:1500]

a_p_small.min(), a_p_small.max()
(0.0, 0.0)

a_g_small.min(), a_g_small.max()
(0.00015188563, 89.998856)
‍‍‍‍‍‍‍‍‍

Good snag... but if you said everything was good elsewhere, you might want to check some of you other samples using my example.  But do check with the distributing agency and/or file a bug

[GarrettBradford · Accepted solution]

ESRI reviewed this case and started a bug: #BUG-000116642

[DanPatterson_Retired]

Garrett... no one can link to your esri account, so you should just edit out the url in your post and leave the bug #

I will mark your answer correct so that people know that a solution is being worked on