Update: The original error was not from the elevation data, but from the use of "Degree". See below.

An agency interested in purchasing heavily sloped land has decided to use a criterion that "slopes must be ≥ 15 percent on more than 50 percent of the parcel". How can I calculate whether a parcel meets this criterion?

The answer may depend on the distance over which one measures the slopes, and on the data source, but I'd expect that any reasonable choice will approximate the agency's calculations well enough for now.

Unfortunately, I got nonsensical results from using GeoElevationData naively (in Mathematica 11.3, on Windows). For instance, I calculated slopes at twelve angles and three distances from a point on the ski trail called Which Way Glade:

whichwayglade = GeoPosition[{42.203, -74.241}];
slopes[x_, radius_] := 
  Table[(GeoElevationData[GeoDestination[x, {radius, 30 i Degree}]] - 
      GeoElevationData[x]), {i, 12}]/radius;
Round[{slopes[whichwayglade, Quantity[100, "Feet"]], 
  slopes[whichwayglade, Quantity[100, "Yards"]],
  slopes[whichwayglade, Quantity[0.4, "Kilometer"]]}, .001]

The output is

{{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
 {-0.022, -0.022, -0.022, -0.022, -0.022, -0.022,
  -0.022, -0.022, -0.022, -0.022, -0.022, -0.022},
 {-0.128, -0.128, -0.128, -0.128, -0.115, -0.115,
  -0.115, -0.115, -0.115, -0.115, -0.115, -0.115}}

The slope at the point should be roughly the maximum absolute value in one of the lists.

Update: Removing the word "Degree" from the code above leads to reasonable output, shown with pretty spaces as:

{{ 0.,     0.295,  0.295,  0.394,  0.066, 0.066,
  -0.23,  -0.23,  -0.295, -0.295, -0.295, 0.},
 { 0.175,  0.241,  0.295,  0.405,  0.284, 0.098,
  -0.109, -0.131, -0.284, -0.317, -0.23, -0.022},
 {-0.013,  0.187,  0.395,  0.402,  0.3,   0.337,
   0.135,  0.002,  0.032,  0.002, -0.1,  -0.128}}

So the slope at Which Way Glade is about 40% using any of these radii.

But compare these outputs with the topo map, GeoDisk[whichwayglade, Quantity[1, "Kilometer"]] // GeoElevationData // QuantityMagnitude // ListContourPlot

enter image description here

The 100-foot slopes would suggest a flat area. The 100-yard slopes would suggest the top of a perfectly conical mountain. The last outputs would suggest that there is no way to go 1/4-km uphill. So I don't think these results are reliable.

How can I calculate the slope at a point in a reasonable way? Or how can I efficiently calculate the percent of a multi-acre parcel with slopes over 15%?

  • 1
    $\begingroup$ Great question, very interesting. I expect that the method actually used is generating a slope map (gdal.org/programs/gdaldem.html) which you can do in QGIS or similar software (or directly in GDAL). The first method that springs to mind to approximate that would be to subdivide the parcel into small areas and measure the difference between the max and min height in each of those small areas. $\endgroup$
    – Carl Lange
    Commented May 20, 2021 at 13:29
  • 2
    $\begingroup$ I've waited many years for a decent skiing-related question (though not as long as I'd have to wait, three+ eternities at least, for me to become a decent skier). $\endgroup$ Commented May 21, 2021 at 0:34
  • 2
    $\begingroup$ Is it agreed you're speaking in percentages of the footprint of the land, not its surface? Assume a silly parcel with a footprint that's for 50% a uniform 45degree slope, then the rest a flat plateau taking the other 50%. Now the surface of the land is (sqrt[2] + 1)/2 times its footprint; so in one measure it's 50% with the desired slope, in the other measure it's 58.6%? Possibly relevant: Matt Parker, "Does <land area> assume a country is perfectly flat?" youtube.com/watch?v=PtKhbbcc1Rc --- spoilers: It depends on what service gives you their data. $\endgroup$ Commented May 21, 2021 at 16:18
  • $\begingroup$ @user3445853, that is interesting — the agency does not specify footprint vs surface in its public documents. But the reason for focusing on slopy land is its potential for runoff that would dirty nearby rivers and reservoirs. So if they think about the issue, they might note that rainfall is more proportional to the footprint, and focus on that measure. $\endgroup$
    – Matt F.
    Commented May 21, 2021 at 17:36

2 Answers 2

whichwayglade = GeoPosition[{42.203, -74.241}]

parcel = GeoBoundingBox[whichwayglade, Sqrt[Quantity[100, "Acres"]]]

We'll get our elevation data as GeoPositions so we can accurately measure the distance between the points (GeoDistance doesn't take height into account, according to the documentation). We can modify GeoArraySize here to get as much data as we have computing power - although be aware that the maximum resolution of the SRTM data GeoElevationData is based on is Not That Great)

elev = GeoElevationData[parcel, Automatic, "GeoPosition", 
  GeoArraySize -> {20, 20}

Now we'll split our positions into 2x2 blocks, and get the difference between the minimum and maximum altitude change, and the distance between the points with the minimum and maximum altitude.

mat = Partition[Thread /@ Thread@elev, {2, 2}];

minmax[b_] := SortBy[Flatten@b, #["Elevation"] &][[{1, 4}]]

Now we get the slope percentage by getting the "rise" - the difference in minimum and maximum elevation - and putting it over the "run" - the GeoDistance between those points that have the minimum and maximum elevation.

percents = 
 MapAt[(Quantity[(#[[2]]["Elevation"] - #[[1]]["Elevation"]), 
       "Meters"]/GeoDistance[#])*100 &, 
  MapAt[minmax, mat, {All, All}], {All, All}]

Now plotting these slope percents, we can see which boxes are flatter than others:

enter image description here

And we can trivially determine how many are over a certain percent:

Length@Select[Flatten@percents, GreaterThan[15]]


(of our total Length@Flatten@percents 100).

Side note: we get a really cool animation by incrasing the resolution over time:

slopeMap[bbox_, res_] := Module[{
   mat = Partition[
     Thread /@ 
      Thread@GeoElevationData[bbox, Automatic, "GeoPosition", 
        GeoArraySize -> {res, res}], {2, 2}],
  minmaxes = 
   MapAt[b |-> SortBy[Flatten@b, #["Elevation"] &][[{1, 4}]], 
    mat, {All, All}];
  MapAt[(Quantity[(#[[2]]["Elevation"] - #[[1]]["Elevation"]), 
        "Meters"]/GeoDistance[#])*100 &, 
   MapAt[minmax, mat, {All, All}], {All, All}]]

plots = Table[
    slopeMap[parcel, n]], {n, {5, 8, 10, 15, 20, 30, 40, 50, 100, 


enter image description here

  • $\begingroup$ Yes, GeoZoomLevel should do it for you. $\endgroup$
    – Carl Lange
    Commented May 20, 2021 at 14:30
  • $\begingroup$ Unfortunately I'm not sure, I haven't used 11.3 since 12 came out. Perhaps the GeoElevationData documentation can help you find out how to get different resolutions. I used a smaller resolution just for testing, I'm not sure what the max is. IIRC SRTM data is about a 10m resolution. $\endgroup$
    – Carl Lange
    Commented May 20, 2021 at 14:39
  • $\begingroup$ Let us continue this discussion in chat. $\endgroup$
    – Carl Lange
    Commented May 20, 2021 at 14:43
  • $\begingroup$ For 11.3, I omitted GeoSize (which ends up roughly quadrupling the data points), redefined minmax and percents, and got an answer of 90.9% instead. The new code is: Elevation[b_] := b[[1, 3]]; elevchange[{a_, b_}] := Quantity[b[[1, 3]] - a[[1, 3]], "Meters"]; minmax[b_] := SortBy[Flatten@b, Elevation][[{1, 4}]]; percents = MapAt[100 elevchange[#]/GeoDistance[#] &, MapAt[minmax, mat, {All, All}], {All, All}]; $\endgroup$
    – Matt F.
    Commented May 20, 2021 at 15:38
  • 3
    $\begingroup$ I'm not surprised you got a similar but not quite the same answer - the scale of the boxes is going to have an outcome on the result. Without a defined scale in the regulations you can pretty much massage the scale to get whatever answer you like above a certain baseline. Part of the fractal nature of the earth, I suppose - extremely similar to measuring the length of a coastline. ☺ $\endgroup$
    – Carl Lange
    Commented May 20, 2021 at 16:04

You have to figure out the resolution of the xy coordinates of the data. The docs say at GeoZoomLevel of 12, the resolution is 38 meter at the equator. What it is elsewhere, it does not say. The below assumes that the resolution is in a certain number of degrees in longitute and latitude, and that remains constant. That's the idea. I'm afraid I cannot personally attest to it accuracy.

With that understanding, it's not hard to computing the gradient from an interpolation:

data = GeoDisk[whichwayglade, Quantity[1, "Kilometer"]] // 
    GeoElevationData // QuantityMagnitude;

if = Interpolation[
           42.203` Degree
           ], 1},
        #1] &,
     data, {2}],
myslope[x_, y_] = D[if[y, x], {{x, y}}] // Norm;

Quantity stuff is slow so I precomputed some constants above. Here's a visualization, where the area outside the red lines has slopes over 15%:

  if[y, x],
   Sequence @@ MapThread[Prepend, {Reverse@if@"Domain", {x, y}}]
  myslope[x, y],
   Sequence @@ MapThread[Prepend, {Reverse@if@"Domain", {x, y}}]
  Contours -> {{0.15}}, ContourStyle -> Directive[Red, Thick], 
  ContourShading -> {None, Directive[Opacity[0.3, Red]]}

enter image description here

Added by Matt F:

This method curiously solves the problem without making it easy to calculate the interpolated height or estimated slope at a particular latitude and longitude. However, the same ideas can be used in an elevation function with lat-long inputs, giving the following parallel for the code in the question:

whichwayglade = GeoPosition[{42.203, -74.241}];
sidelength = Quantity[Sqrt[Quantity[100, "Acres"]/
   Quantity[1., "Meters"]^2], "Meters"];
parcel = GeoBoundingBox[GeoDisk[whichwayglade, sidelength]];
data = parcel // GeoElevationData // QuantityMagnitude;
c2 = Quantity[38., "Meter"]/Quantity[1, "Foot"];
c1 = c2 Sin[42.203 Pi/180];
d1 = c1 (Dimensions[data][[1]] - 1)/(parcel[[2, 1, 1]] - parcel[[1, 1, 1]]);
d2 = c2 (Dimensions[data][[2]] - 1)/(parcel[[2, 1, 2]] - parcel[[1, 1, 2]]);
if = Interpolation[Flatten[MapIndexed[N@Append[#2 {c1, c2}, #1] &, data, {2}], 1]];
newelev[geopos_] := Quantity[if[c1 + d1 (geopos[[1, 1]] - parcel[[1, 1, 1]]),
  c2 + d2 (geopos[[1, 2]] - parcel[[1, 1, 2]])], "Feet"];
slopes[x_, radius_] := Table[newelev[GeoDestination[x, {radius, 30 i}]]
  - newelev[x], {i, 12}]/radius;
Round[{slopes[whichwayglade, Quantity[100, "Feet"]], 
  slopes[whichwayglade, Quantity[100, "Yards"]], 
  slopes[whichwayglade, Quantity[0.4, "Kilometer"]]}, .001]
  • $\begingroup$ @MattF. Sorry about data. It's your data, but I forgot you didn't call it data. Carl gets the data with the xy coordinates in real units. I estimate the relative coordinates in real units. I get the gradient from an interpolation. Carl estimates the gradient from the corner points (I think, and I think that could be improved slightly, which I haven't had time to check it out). $\endgroup$
    – Michael E2
    Commented May 20, 2021 at 21:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.