How is slope calculated by grdgradient?

Hello,
I want to make a slope map from a global geographic digital elevation model in simple cylindrical projection (x/y coordinates are lon/lat and z is the surface height is in kilometers).

I believe grdgradient is the easiest way to do this, but I would like to better understand how is the slope (the “slopegrid” for the -S flag) calculated by grdgradient?

Is it just sqrt( (dz/dx)^2 + (dz/dy)^2 ) where dz/dx and dz/dy are the x and y gradients computed from a finite difference across the central pixel of a 3x3 moving window?

I have read the documentation, and my understanding is that, if the -fg flag is specified, then the dx and dy distances are computed with a flat Earth approximation, correct? What if the -fg flag is not specified?

This is a map of the Moon, so if I specify the radius with --PROJ_ELLIPSOID=1737.4 km, then will grdgradient correctly scale the dx and dy?

Finally, the output slope would be in radians?

Thank you for any help!

Hello mkbarker

I want to make a slope map from a global geographic digital elevation model in simple cylindrical projection (x/y coordinates are lon/lat and z is the surface height is in kilometers).

I believe grdgradient is the easiest way to do this, but I would like to better understand how is the slope (the “slopegrid” for the -S flag) calculated by grdgradient?

Is it just sqrt( (dz/dx)^2 + (dz/dy)^2 ) where dz/dx and dz/dy are the x and y gradients computed from a finite difference across the central pixel of a 3x3 moving window?

Not sure but I think that follow the method described by Horn (1981; http://people.csail.mit.edu/bkph/papers/Hill-Shading.pdf)

I have read the documentation, and my understanding is that, if the -fg flag is specified, then the dx and dy distances are computed with a flat Earth approximation, correct? What if the -fg flag is not specified?

Then it assumes that X,Y,Z are all in the same units.

This is a map of the Moon, so if I specify the radius with --PROJ_ELLIPSOID=1737.4 km, then will grdgradient correctly scale the dx and dy?

Not sure, but I think so.

Finally, the output slope would be in radians?

Mm, no. The result is the magnitud of the gradient vector.

Check these commands (from this script; FCEN-2022/1_Ejercicios/05B_Pendiente/12_Pendiente_Grados.sh at main · Esteban82/FCEN-2022 · GitHub).

#	Calculate gradient from geographic grid (-fg)
	gmt grdgradient $DEM -D -S$GRAD -fg
#	Convert gradient module to inclination (slope)
	gmt grdmath $GRAD ATAN = $RAD  # Convert to randians
	gmt grdmath $GRAD ATAND = $DEG  # Convert to degrees
1 Like

Maybe this will help.

The ‘direction’ grid (as opposed to the magnitude grid), will be in radians I think.

1 Like

Yes, according to this line the planetary body is taken into consideration

dx_grid = GMT->current.proj.DIST_M_PR_DEG * Grid->header->inc[GMT_X] * cosd (lat);
1 Like

But there is nothing of GMT.jl in this question, right?

I wasn’t sure what “GMT.jl” meant…

I made a How-to guide for making a slope map of the moon.

https://www.generic-mapping-tools.org/gmt-examples/gallery/images/slope.html

2 Likes