What fraction of a mountain parcel has slopes over 15%?

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:

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

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

88

(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
   },
  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[
   ArrayPlot[
    slopeMap[parcel, n]], {n, {5, 8, 10, 15, 20, 30, 40, 50, 100, 
     200}}];

ListAnimate[plots]

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[
   Flatten[
    MapIndexed[N@Append[#2*
         (*Quantity[38.,"Meters"]/Quantity[1,"Feet"]*)
         38*3.2808398950131235`*
         {Sin[(*whichwayglade//Latitude*)
           42.203` Degree
           ], 1},
        #1] &,
     data, {2}],
    1]
   ];
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%:

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

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]