Calculate gradient, slope and aspect of data grid
[ASPECT, SLOPE, gradN, gradE] = gradientm(Z,
[...] = gradientm(lat, lon, Z)
[...] = gradientm(..., ellipsoid)
[...] = gradientm(lat, lon, Z, ellipsoid, units)
[ASPECT, SLOPE, gradN, gradE] = gradientm(Z, R) computes the slope, aspect and north and east components of the gradient for a regular data grid Z with three-element referencing vector refvec. If the grid contains elevations in meters, the resulting aspect and slope are in units of degrees clockwise from north and up from the horizontal. The north and east gradient components are the change in the map variable per meter of distance in the north and east directions. The computation uses finite differences for the map variable on the default earth ellipsoid.
R can be a geographic raster reference object, a referencing vector, or a referencing matrix.
If R is a geographic raster reference object, its RasterSize property must be consistent with size(Z).
If R is a referencing vector, it must be a 1-by-3 with elements:
[cells/degree northern_latitude_limit western_longitude_limit]
If R is a referencing matrix, it must be 3-by-2 and transform raster row and column indices to or from geographic coordinates according to:
[lon lat] = [row col 1] * R
If R is a referencing matrix, it must define a (non-rotational, non-skewed) relationship in which each column of the data grid falls along a meridian and each row falls along a parallel. Nearest-neighbor interpolation is used by default. NaN is returned for points outside the grid limits or for which lat or lon contain NaN. All angles are in units of degrees.
[...] = gradientm(..., ellipsoid) uses the reference ellipsoid specified by the input ellipsoid, which can be a referenceSphere, referenceEllipsoid, or oblateSpheroid object, or a vector of the form [semimajor_axis eccentricity]. If the map contains elevations in the same units of length as the semimajor axis of the ellipsoid, the slope and aspect are in units of degrees. This calling form is most useful for computations on bodies other than the earth.
[...] = gradientm(lat, lon, Z, ellipsoid, units) specifies the angle units of the latitude and longitude inputs. If omitted, 'degrees' are assumed. For elevation maps in the same units as ellipsoid(1), the resulting slope and aspect are in the specified units. The components of the gradient are the change in the map variable per unit of length, using the same length unit as the semimajor axis of the ellipsoid.
Compute and display the slope for the 30 arc-second (10 km) Korea elevation data. Slopes in the Sea of Japan are up to 8 degrees at this grid resolution.
load korea [aspect, slope, gradN, gradE] = gradientm(map, refvec); worldmap(slope, refvec) geoshow(slope, refvec, 'DisplayType', 'texturemap') cmap = cool(10); demcmap('inc', slope, 1, , cmap) colorbar latlim = getm(gca,'maplatlimit'); lonlim = getm(gca,'maplonlimit'); land = shaperead('landareas',... 'UseGeoCoords', true, 'BoundingBox', [lonlim' latlim']); geoshow(land, 'FaceColor', 'none') set(gca, 'Visible', 'off')
Coarse digital elevation models can considerably underestimate the local slope. For the preceding map, the elevation points are separated by about 10 kilometers. The terrain between two adjacent points is modeled as a linear variation, while actual terrain can vary much more abruptly over such a distance.