2

I'm trying do a few different things here..

First, I have an array of values from the watershed polygon 'shape' field in a shapefile.

Rarray = watershed.shape.getPart(0)

I'm making the assumption that the outflow point will be the value with the lowest ZMin. So, ZMin coordinates will be the outflow point (p1).

What I am attempting to do is find the point (p2) within the polygon array that is furthest from this x,y,z outflow point. It should be one of the values within the array

From there, I am trying to calculate the distance between p1 (outflow) and p2(furthest away) so that I can use that value to calculate the relief ratio of the watershed using this formula

ReliefRat = (max elevation - min elevation) / Length of longest axis

So far I have this...

Rarray = watershed.shape.getPart(0)
ReliefRat = (ZMax-ZMin)/(((p2.X-p1.X)**(2.0)) + ((p2.Y-p1.Y)**(2.0)))**(0.5)

...Where p1 is the outflow point. I just don't know how to find p2.

If there is somebody out there that can walk me through this it would be greatly appreciated!

Jen K.
  • 21
  • 2

4 Answers4

0

If p1 is indeed an object containing the x,y values of the outflow point with ZMin, and if p2 contains x,y for the ZMax farthest point, then the expression

(((p2.X-p1.X)**(2.0)) + ((p2.Y-p1.Y)**(2.0)))**(0.5)

is indeed the distance between the two points. This is just an expression of the Pythagorean theorem, because the straight line between two points p1 and p2 is the hypotenuse of a right triangle with base p2.X-p1.X and height p2.Y-p1.Y. Note, see “Which is faster in Python: x**.5 or math.sqrt(x)?” if interested in several ways to compute square roots.

If you haven't already defined a class of point objects with X and Y attributes, ask about that separately or amend the question, and likewise if you don't know how to find the point on a polygon most distant from some other point.

Community
  • 1
  • 1
James Waldby - jwpat7
  • 8,593
  • 2
  • 22
  • 37
  • I have a polygon in a shapefile, listed as 'Polygon ZM' in the shape field. I am able to get an array of points by using...array = watershed.shape.getPart(0)...from there I then have the ZMin and ZMax values from the polygon right? – Jen K. Mar 09 '13 at 00:45
  • I'm using the z-min coordinates as my outflow point, so I need to find the coordinates in the array that are the furthest distance (within the polygon) from the outflow point. I hope I'm asking it the right way – Jen K. Mar 09 '13 at 00:51
  • @JenK., instructions to find `p2` are now in my answer. – askewchan Mar 09 '13 at 05:39
0

It seems that p1 should be the point of outflow, but you still have to somehow determine what p2 is, which is the point in your watershed polygon that's furthest away from p1.

To find p2, I'll assume you know the points that are the vertices of the watershed polygon. Try to make a list of them, which I'll call poly. You want to find the point on the polygon boundary that is furthest from the outflow point, p1:

import numpy as np
poly = np.array([(1.1,5.3),(1.5,2.3),(2.1,3.2),(4.3,4.4)]) # vertices of a polygon
p1 = np.array((1.2,3.5))     # outflow point
ds = map(norm,p1 - poly)     # list of distances between p1 and each vertex of polygon
d_max = max(ds)              # max of those distances
ReliefRat = (ZMax-ZMin)/d_max

If you already know p2, you can find the distance between p1 and p2 quite easily. For better readability, I would recommend using the norm function from numpy, or else write your own:

from numpy import norm
ReliefRat = (ZMax-ZMin)/norm((p2.X-p1.X,p2.Y-p1.Y))

depending on how p1 and p2 are formatted (e.g., as numpy arrays), this could be written as:

ReliefRat = (ZMax-ZMin)/norm(p2 - p1)
askewchan
  • 45,161
  • 17
  • 118
  • 134
  • I'm having trouble finding p2, the furthest point from the outflow point. I made the assumption that the outflow point would be the x,y,z point in the polygon array that had the lowest z-value, thus the minimum elevation. So I have that point. I want to find the point in the polygon array that is the furthest distance from the outflow point. Does that help? – Jen K. Mar 09 '13 at 00:49
  • If you know the vertices of the polygon, you can just check each one, and calculate the distance, and find the one with the largest distance. Then, the point `i` that has the maximum distance from p1 (using the `norm(p_i - p1)`) will be the one you want to use as `p2`. – askewchan Mar 09 '13 at 05:26
0

Not sure I understand the concepts about relief ratio of watersheds and so on, but get a feeling this is about finding distances between points.
I find using complex numbers for this purpose is convenient. Python has default support for it. I understand you have a point (lowest point), and a number of nodes in a x,y coordinate system:

>>> # Lowest point:
... z0 = complex(5, 8)
>>> # A few nodes:
... z1 = complex(5, 10)
>>> z2 = complex(-2, 15)
>>> abs(z1 - z0)
2.0
>>> # Abs give a scalar so order is not important:
... abs(z0 - z1)
2.0
>>> abs(z2 - z0)
9.8994949366116654
>>> # The vector between lowest point z0 and node z2 was the longest axis

By doing z0.real and z0.imag you extract the x and y values respectively. enter image description here

If a lot of nodes like this are stuffed in a list it will be a matter of iterating and searching for the maximum difference.

Tompa
  • 5,131
  • 2
  • 13
  • 13
0

Finding p2

If you also have z-coordinates as elevation, maybe the problem is geometrical. Finding the straight line distance between two points in a 3d space requires a z component. Assume you have some help code like this:

class P():
    """A point in 3d space."""
    def __init__(self, x = 0.0, y = 0.0, z = 0.0):
        self.x = float(x)
        self.y = float(y)
        self.z = float(z)

    def __str__(self):
        return str(self.x) + ', ' + str(self.y) + ', ' + str(self.z)

    def __repr__(self):
        return str(self.x) + ', ' + str(self.y) + ', ' + str(self.z)

def dist(p0, p1):
    """Return the straight distance in 3d space between point p0 and p1
    as an absolute value."""
    x0, y0, z0 = p0.x, p0.y, p0.z
    x1, y1, z1 = p1.x, p1.y, p1.z

    return ((x1 - x0)**2 + (y1 - y0)**2 + (z1 - z0)**2)**(1 / 2.0)

Then I make some points in space:

>>> p1 = P(5, 5, -3) # Outflow point
>>> p2 = P(40, -20, 45)
>>> p3 = P(-5, 30, 40)
>>> p4 = P(-2, 25, 22)
>>> arr = [p2, p3, p4] # An array (list) of 3d points.
>>> p2find = p1 # For the search.

>>> for point in arr:
...    if dist(p1, point) > dist(p1, p2find):
...       p2find = point
... 
>>> p2find
40.0, -20.0, 45.0
>>> p2
40.0, -20.0, 45.0

And so the point p2 being furthest away should be found.

Tompa
  • 5,131
  • 2
  • 13
  • 13