Map water runoff flows near Ellicott City, Maryland

Question

Ellicott City, Maryland was the site of a flash flood several years ago due to the local topography. How can we use the Wolfram Language to show the direction water flows at many points near that city? A good start is code by Roman from the answer here:

center = Entity["City", {"EllicottCity", "Maryland", "UnitedStates"}];
range = Quantity[4, "Miles"];
elevationData = 
  GeoElevationData[
    center, GeoRange -> range, 
    GeoProjection -> Automatic, UnitSystem -> "Imperial"
  ];

I envision using ListVectorPlot to show many vectors that are perpendicular to the equal elevation contours (all pointing downhill). The hard part is making the array of vectors from elevationData.

Answer 1

Just an idea by smoothing the elevation data with MeanFilter and using grad function from @xzczd's answer:

grad[mat_] := grad[mat, ##] & @@ ConstantArray[1, Length@Dimensions@mat]
grad[mat_, dx__] := 
  1/{dx} (NDSolve`FiniteDifferenceDerivative[#, 
       Range@N@Dimensions@mat, mat, DifferenceOrder -> 2] & /@ 
     IdentityMatrix@Length@{dx});
(* Change to control the smoothness *)
smoothingRadius = 16;
elevationDataSmooth = 
  MeanFilter[QuantityMagnitude@elevationData, {smoothingRadius, smoothingRadius}];
vec = Map[{#[[1]], -#[[2]]} &, 
   Reverse[Transpose[-grad[elevationDataSmooth, 100, 100], {3, 2, 1}], 
    {2, 3}], {-2}];
Show[ReliefPlot[elevationData, DataReversed -> True, 
  ColorFunction -> (Opacity[.5, ColorData["GreenBrownTerrain"][#]] &)], 
 ListStreamPlot[vec, StreamPoints -> Fine]]

Answer 2

You could try to fit an interpolation function and then take the negative Gradient of it.

First we must get rid of the unit:

dat=QuantityMagnitude[elevationData];

Then we get an interpolation function:

inter = ListInterpolation[dat]

Now we can take the Gradient:

grad = Grad[inter[x, y], {x, y}]

Then we can plot the gradients:

vecplot = VectorPlot[grad, {x, 1, 438}, {y, 1, 436}]

and make a contour plot of the landscape:

cp = ContourPlot[inter[x, y], {x, 1, 438}, {y, 1, 436}, 
  ColorFunction -> (White &)]

Finally we plot both together:

Show[cp, vecplot]