# Map water runoff flows near Ellicott City, Maryland

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.

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]
1/{dx} (NDSolveFiniteDifferenceDerivative[#,
Range@N@Dimensions@mat, mat, DifferenceOrder -> 2] & /@
IdentityMatrix@Length@{dx});

(* Change to control the smoothness *)

elevationDataSmooth =

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]]


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]
`