1
$\begingroup$

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.

$\endgroup$

2 Answers 2

5
$\begingroup$

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

enter image description here

$\endgroup$
3
$\begingroup$

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

![enter image description here

and make a contour plot of the landscape:

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

enter image description here

Finally we plot both together:

Show[cp, vecplot]

![enter image description here

$\endgroup$

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.