Subscribe to this thread
Home - General / All posts - Computing minimun slope
vincent

1,767 post(s)
#08-Feb-19 16:55

Hi,

I would like to obtain the minimun slope instead of the maximum slope, from a DEM, for all the pixels.

Any idea ?

I have M8 and M9. It can be with SQL, Surface Transform or anything else.

I think to use all 8 "Difference" SQL functions within a window, get the minimun difference out of the results, get the Height value of the center pixel and calculate the slope for the difference between the 2 values and the distance between pixel center.

Dimitri


5,359 post(s)
#08-Feb-19 18:18

I would like to obtain the minimun slope instead of the maximum slope, from a DEM, for all the pixels.

The Slope transform does that now. It generates whatever is the slope at all pixels, the maximum, the minimum and everything in between.

vincent

1,767 post(s)
#08-Feb-19 19:20

Perhaps my need was wrongly explained. I understand that Slope Transform return the slope, which is the maximum one in my idea for an x,y location. I'm after the minimum one, wihich is likely to be perpendicular to the maximum one.

The GIS computed slope is typically

the maximum change in elevation over the distance between the cell and its eight neighbors identifies the steepest downhill descent from the cell.

I need the slope of the minimum change.

Dimitri


5,359 post(s)
#09-Feb-19 05:38

Take a flat plane and tilt it. For each pixel in that plane, the slope is in the direction the plane tilts. The "minimum" slope is 0, because there is no tilt when you walk perpendicular to the slope. It's like a ski run: beginners who cannot prevent themselves from speeding up too much if they point themselves straight down the hill can get off the ski run by shuffling perpendicular to the slope, a traverse.

Strictly speaking, at the purely local level where "slope" is defined, every pixel is at an inclined flat plane, and thus the "minimum slope" for each pixel is zero. (I have worked out a marvelous proof of that proposition but the margin of this forum is too small to write it.)

A surface which is curved in two directions is not a flat plane. See the CalTech article cited for Curvature, Mean in the Transform Images documentation. At any given spot on the surface, you can compute the Normal vector N, that is the direction pointing perpendicular from the surface. What you may be looking for is the inclination calculated in the minimum curvature direction (I forget if it is the x1 or x2 direction in the illustration) perpendicular to N, the normal, at any given pixel.

Another proxy you could use quickly with a point and click would probably be Curvature, Plan, which gives the the curvature perpendicular to the line of maximum slope. See the illustrations in the blog entry cited in the Transform entries for Plan and Profile curvatures. But that assumes a special case where surfaces are curved and not even (a flat plane that is tilted has slope, but not curvature).

Directional Edge detectors (see this link) basically just compute slopes in a given direction. People are usually interested in maximum slopes since that is what happens when values change rapidly, as at edges. But if you wanted minimum slopes you could do something (just guessing here) like adapt the Kirsch example to use TileMin instead of TileMax, to get the minimum value of each of the eight 45 degree directional filters.

All of the above uses 9 or Viewer. If you don't have 9, you can play around with Viewer to try out ideas.

tjhb

8,657 post(s)
#09-Feb-19 06:24

I don't think so.

The "minimum" slope is 0, because there is no tilt when you walk perpendicular to the slope.

This is false, unless we first assume an image having uniform gradient in some direction.

That is far from the usual case for terrain, even locally, and not a good generalization in advance of data.

For any unknown terrain it may well be that the direction of maximum slope (the aspect) is not perpendicular to... anything else. There may in fact be no "traverse"--no pair of neighbours either side of this pixel which are at the same or a similar level--and if there happens to be such a pair, there is no a priori reason to think that they will be at right angles to the steepest descent or the steepest ascent.

(By the same token, there is no reason to suppose that the steepest ascent and the steepest descent "should be" in line. For real terrain they will often not be. We are not just open to peaks and pits, but also to random, noisy, or just "detailed" terrains. )

While it is helpful for many individual purposes to imagine a tangent plane, it is a categorical mistake to address different questions from the same assumed plane without a shared reason. Except in exceptional circumstances, or by accident, the same imaginary plane will not apply.


Vincent is asking a simple question--exactly as simple as asking in which direction the (steepest) slope lies.

In my opinion his suggested answer is roughly correct, but we need to deal with perpendicular neighbours ("rook's move") separately from diagonal neighbours ("bishop's move"), since their radii are different. Their slopes must first be ranked within their own sets, then between the two sets.

I have had a go at this today, for 9 (and GPGPU), and I think it works, but so far only on paper. More tomorrow.


By the way it would be fantastic if Manifold could reference its method for calculating the (maximum) slope, in the manual. There are several different methods, with potentially wide variations. (And for the present purpose, it would be great to be able to write a method for minimum slope to match.)

Dimitri


5,359 post(s)
#09-Feb-19 11:47

I don't think so.

The "minimum" slope is 0, because there is no tilt when you walk perpendicular to the slope.

This is false, unless we first assume an image having uniform gradient in some direction.

Just nit-picking, but what I wrote was true. You left off the part of what I wrote where I explicitly stated there was a uniform gradient (the "Take a flat plane" bit). What I wrote:

Take a flat plane and tilt it. For each pixel in that plane, the slope is in the direction the plane tilts. The "minimum" slope is 0, because there is no tilt when you walk perpendicular to the slope.

A lot of these issues have to do with the scale at which you examine a terrain. The Earth, for example, when seen from afar by a godlike giant is smoother than a billiard ball. We all think of a billard ball as smooth, but in reality when you get close in there are all sorts of ridges and valleys, just like with the Earth.

Same with terrain. The calculation of what gives a smooth, or apparent "slope" depends on what scale you look at, over which you average.

adamw


8,447 post(s)
#21-Feb-19 13:01

By the way it would be fantastic if Manifold could reference its method for calculating the (maximum) slope, in the manual.

We compute the rate of change of the raster in horizontal (D) and vertical (E) directions, then compute an arc tangent of sqrt(EE+DD).

For a 3x3 window = radius 1, with heights numbered from z1 to z9 (top down and then left to right), D is computed as (z3 + z6 + z9 - z1 - z4 - z7) / 6, and E is computed as (z1 + z2 + z3 - z7 - z8 - z9) / 6. That's slightly different from what ArcGIS does - they emphasize the center horizontal / vertical by using (z3 + 2 * z6 + z9 - z1 - 2 * z4 - z7) / 8, however what we do generalizes into higher radius values more naturally.

tjhb

8,657 post(s)
#21-Feb-19 15:22

Thank you for this! Very valuable indeed.

(Can a note be added to the manual at some point?)

P.s. Interesting that the notion of a local peak, pit, or saddle simply does not arise on this definition. The centre value z5 is treated as noise for the purpose.

Also that there is no discounting for diagonals. (But circles are not always better than squares.)

P.p.s. For higher radius values, do you continue to take into account the nearer neighbours, or only use the outermost rows?

adamw


8,447 post(s)
#21-Feb-19 15:47

For higher radius values, we compose a best-fit approximation of the surface using all neighbors in the square. The computations coincide exactly with their simplified versions above for radius 1.

This will be in the API doc (being long promised).

tjhb

8,657 post(s)
#21-Feb-19 16:53

Got it, the generalizing makes good sense now. Thanks again.

tjhb

8,657 post(s)
#09-Feb-19 06:43

I other words:

Strictly speaking, at the purely local level where "slope" is defined, every pixel is at an inclined flat plane...

That is false. It is only an inclined flat plane as an approximation for the sole purpose of calculating slope.

The same approximation does not apply to anything else--strictly speaking, not even to calculating aspect.

A surface which is curved in two directions is not a flat plane.

That tautology is (of course) true, but most terrains are locally hyperbolic to some degree, which is often significant--though even more often, local curvature fights with measurement noise. Assuming that every departure from local flatness is noise is at least as bad as assuming that every minor departure is true local form.

Dimitri


5,359 post(s)
#10-Feb-19 10:54

It is only an inclined flat plane as an approximation for the sole purpose of calculating slope.

Correct. Which is the only context in which "slope" is defined. You're basically saying that "slope only means what it means when you calculate it to mean what it means." Well, sure. If you calculated to mean something else, it wouldn't be "slope."

The slope at a point on a curve is the tangent to that point.

The slope at a point on a curved surface is the inclination of the flat plane tangent to that point.

In both cases, you use the classic calculus effect of a limit of the series where shrinking the calculation of points from far apart to infinitely close generates a progression from secant to tangent.

You're basically saying "I'm looking at the limit of more microscopic considerations where, in effect, what looks like a curved surface at bigger scales now looks like a flat, inclined plane."

My statement is true: "slope" is an artificial notion, the tangent, that is an approximation to reality at some scale. What is the "slope" of a hyperbolic surface? The surface doesn't have an overall "slope." There is only "slope" at various local spots on that hyperbolic surface, where some approximation equivalent to a tangent plane can be computed.

For any location on a surface, you can imagine a vertical, flat plane, perfectly aligned in the vertical direction with the Normal (z direction) of the coordinate system. As you rotate that plane about the Z vertical axis that passes through the point, you see how where that vertical plane cuts the surface you get a line. Compute the tangent to that line and you get a slope for each such possible line. As you rotate that plane, beginning with the rotation where it gives the maximum slope for the line of intersection, you'll see the slope go down in value until reaches 0.

The proof, by the way, that at every location for which you can compute a non-zero "maximum" slope the "minimum" slope is zero can be visualized without any math, as a simple thought experiment.

Consider a surface covered by a grid of discrete point locations. For each point location, a slope value has been computed. If the surface were a flat XY plane, that is, perpendicular to the Z axis of the Euclidean coordinate system, the slope at all locations would be zero. If the surface has undulations, hills, valleys and saddle points, the slope at various locations might be zero or non-zero.

Suppose one of the undulations of the surface is a hill shaped like a camel's hump. At one of the locations on the side of that hill the slope is 45 (we'll use degrees for conceptual ease) in whatever you consider to be the "straight downhill" direction. Super. Let's teleport ourselves onto the hill so we are standing at that spot looking straight downhill, where we measure a 45 degree angle that is the slope. You've been miniaturized so where you stand looks like a flat, inclined plane (all parts of that surface at a small enough scale look flat).

Suppose you now aim yourself to the left a bit, instead of straight downhill but more in the way of a traverse. The slope angle will be less than 45. As you continue moving your gaze to the left and up, counter-clockwise, the slope angle will be less and less until at some point you reach 0. Keep on going and you end up with a repetition in mirror image, just as when you stand on that microscopic spot and either look directly downhill to get 45 or directly uphill to also get 45 (it's not said to be a -45 slope when you look uphill or a +45 slope when you look downhill).

That zero slope direction, by the way, is the direction of the perpendicular to the normal of the coordinate system, where a flat XY plane with zero slope (another way of saying it is perpendicular to the normal) would cut the surface of the camel's hump.

Where people get confused by seemingly tricky surfaces like concave, convex or saddle points is they forget slope is computed to get the tangent plane, in effect looking at such small regions where seeing the surface as a tangent flat plane makes sense. However you calculate, the ultimate notion of "slope" as the tangent flat plane remains the same in all of them.

adamw


8,447 post(s)
#21-Feb-19 13:14

It is only an inclined flat plane as an approximation for the sole purpose of calculating slope.

Agree, the plane is just a local approximation used to compute things like maximum slope = slope. But how to define "minimum" slope? I am just interested what it is we are after here.

adamw


8,447 post(s)
#21-Feb-19 13:07

Regarding minimum slope, I am not sure how to define it. "Maximum" slope is computed from the surface being approximated locally by a plane. If we define "minimum" slope as "take all possible directions, compute a slope in each, then take the minimum value", then, with the approximation being a plane there will always be a direction in which the slope is zero.

What is it specifically that you want to compute? A measure of how rocky the surface is in terms of "whatever direction you take, it is going to be at least this rocky, so if this minimum value exceeds some threshold [do or don't do something here]"? If so, that seems like something resembling curvature.

tjhb

8,657 post(s)
#21-Feb-19 15:07

I suspect you have the right hunch about curvature (it could even be our remarkable friend Gaussian, which holds all the trees up by their branches).

But is it true that there is always some local plane for which slope is zero? Yes for a saddle; yes for a symmetrical peak or pit (the tangent plane right here)—but for an asymmetrical peak or pit?

OK yes again. As Dimitri has said (paraphrasing), we can’t define a slope unless we draw a tangent plane; we can’t draw a tangent plane unless we (at least implicitly) draw a normal; if we draw a normal, we also define (implicitly) an orthogonal direction in which slope is zero.

But there is one much simpler approach as well. How about if (for some reason) we just want the minimum slope in a D8 sense? That is, no fitting, no concept of continuous direction, just looking at the 8 immediate neighbours. There is always one who has the nicest garden, one who just seems to love weeds.

That could be useful in a context where D8 was already being used (or by extension D-infinity). For example if asking something like, how much rain before this soil becomes waterlogged?

[Added. But in that case profile curvature may be a better measure.]

I don’t know the purpose here though.

adamw


8,447 post(s)
#21-Feb-19 15:54

If we define minimum slope as the minimum of 8 slopes from center to its neighbors, then sure, that's easy to compute and perhaps useful. We don't have a built-in function to compute that, but writing a custom function that does it isn't difficult. Regular slope wouldn't be the maximum slope in this sense, but perhaps that isn't important.

adamw


8,447 post(s)
#21-Feb-19 16:23

Example code (the biggest issue is that it assumes all pixels are visible - this has to be adjusted).

Script (C#):

// C#

 

using M = Manifold;

 

class Script

{

 

static Manifold.Context Manifold;

static void Main()

{

}

 

// assumes float64 pixels

// assumes all pixels are visible

// computes min of 8 slopes from center to neighbors

//

 

// slopemin for pixel

static double SlopeMinValue(double beg, double end, double dx)

{

  return (end - beg) / dx;

}

static double SlopeMinValueInc(double slope, double beg, double end, double dx)

{

  return System.Math.Min(slope, (end - beg) / dx);

}

static double SlopeMinPixel(M.Tile.PixelSet<double> pixels, int x, int y)

{

  double center = pixels[x, y];

  double dx = 1;

  double dxdiag = 1.4142; // sqrt(2)

  double value;

 

  // use horizontal / vertical neighbors

  value = SlopeMinValue(center, pixels[x-1, y], dx);

  value = SlopeMinValueInc(value, center, pixels[x+1, y], dx);

  value = SlopeMinValueInc(value, center, pixels[x, y-1], dx);

  value = SlopeMinValueInc(value, center, pixels[x, y+1], dx);

 

  // use diagonal neighbors

  value = SlopeMinValueInc(value, center, pixels[x-1, y-1], dxdiag);

  value = SlopeMinValueInc(value, center, pixels[x-1, y+1], dxdiag);

  value = SlopeMinValueInc(value, center, pixels[x+1, y-1], dxdiag);

  value = SlopeMinValueInc(value, center, pixels[x+1, y+1], dxdiag);

 

  return value;

}

 

// slopemin for tile

static public M.Tile SlopeMin(M.Tile t, int radius)

{

  M.Application app = Manifold.Application;

  M.TileBuilder builder = app.CreateTileBuilder();

 

  // start with blank tile of the same size as source

  builder.StartTile(t.Width, t.Height, t.Type);

 

  // set pixels

  M.Tile.PixelSet<double> pixelsOld =

    (M.Tile.PixelSet<double>)t.Pixels;

  M.TileBuilder.PixelSet<double> pixelsNew =

    (M.TileBuilder.PixelSet<double>)builder.Pixels;

  for (int y = radius; y < builder.Height - radius; ++y)

  {

    for (int x = radius; x < builder.Width - radius; ++x)

      pixelsNew[x, y] = SlopeMinPixel(pixelsOld, x, y);

  }

 

  // return composed tile

  return builder.EndTile();

}

 

}

Query:

--SQL9

FUNCTION SlopeMin(@t TILE, @r INT32TILE AS SCRIPT [SlopeMin]

  ENTRY 'Script.SlopeMin';

 

VALUE @t TILE = StringJsonTile('[2,2,3,2,1,2,3,3,3]', 3, 3, 1, false);

VALUES (TileJson(SlopeMin(@t, 1)));

MXB attached. Hope this is useful as a starting point.

Attachments:
slopemin.mxb

tjhb

8,657 post(s)
#21-Feb-19 16:40

That is a great start not just for this but for heaps of other stuff.

The simple split into SlopeMinValue and SlopeMinValueInc (probably second nature to you) is a helpful design lesson in itself.

Thanks!

vincent

1,767 post(s)
#25-Feb-19 18:46

Thank you Adam. I will look if that can be useful. You were right, I turned my mind to various curvature after posting this at first. They are really useful.

Manifold User Community Use Agreement Copyright (C) 2007-2017 Manifold Software Limited. All rights reserved.