Skip to content

Instantly share code, notes, and snippets.

@wassname
Last active August 29, 2015 14:24
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save wassname/7bb22aa7e36dd8a53a0f to your computer and use it in GitHub Desktop.
Is a line open to the sea
"""
Determine if a line is open to sea by using aspect derived from a digital
elevation model, and comparing to line angle and the change in line angle.
It works like this, if we are turning in the same direction as downhill then it is concave.
The magnitude of concavity is the component of the downhill direction that is perpendicular to direction
Then we do an average of concavity over the whole line, weighted by segment length
"""
import numpy as np
import pylab
from matplotlib import pyplot as plt
from numpy import *
from pylab import *
from scipy.stats import circmean
pi = np.pi
sin = np.sin
cos = np.cos
ar = np.array
def to360(angle):
"convert angle range of -179 to 180 to 0 to 360"
return (angle+360)%360
def from360(angle):
"convert angleset range of 0 to 360 to -179 to 180"
angle[angle>180]-=360
return angle
def angles(line):
"""get angles between points in line in degrees"""
x,y=line.T
ang=np.arctan2(np.diff(x),np.diff(y))
ang=to360(np.rad2deg(np.unwrap(ang)))
return np.append(ang,0)
def line_bins(line):
"""Work out the bin length for a line. The bin is the distance from halfway
to the previous point to halfway to the next point"""
seg_l=pylab.distances_along_curve(line)
seg_l=np.hstack((0,seg_l))
left_bin=seg_l/2
right_bin=np.concatenate([seg_l[1:],seg_l[:1]])
bin_l=left_bin+right_bin
# normalise
bin_l=bin_l/sum(bin_l)
return bin_l
def weighted_line_average(line,record):
"""
Do a weighted average for an attribute on a line. The average is weighted by bin length
where the bin is the distance from halfway to the previous point to halfway to the next point
"""
# get bin length
bin_l=line_bins(line)
# now do weighted mean
return np.average(record,weights=bin_l)
def howOpen(line,downhill):
"""
Determine if a line is open to sea by using aspect derived from a digital
elevation model, and comparing to line angle and the change in line angle.
It works like this, if we are turning in the same direction as downhill then it is concave.
The magnitude of concavity is the component of the dowhill direction that is perpendicular to direction
Then we do a average of concavity over the whole line, weighted by segment length
"""
downhill=to360(downhill)
# get the angle for each segment
segment_angle=angles(line)[:-1]
segment_angle=to360(segment_angle)
# get the angle for each node as the average of the segment before and after
node_angle=ar([circmean(ar([a,b]),high=360,low=0) for a,b in zip(segment_angle[1:],segment_angle[:-1])])
node_angle=to360(node_angle)
# is the line turning left or right at each node?
turnRight=np.diff(segment_angle)/np.abs(np.diff(segment_angle))
# Get the component of the downhill direction that is orthogonal to the node direction
weight=-sin(np.deg2rad(downhill[1:-1]-node_angle)) # how much is the line facing downhill, 0-1
# Now get average of weight*turnDirection and weight the average by length of segments on each side of a node
nodeOpen=-weight*turnRight
lineOpen=np.concatenate(([0],nodeOpen,[0]))
openness=weighted_line_average(line,lineOpen)
return openness
import unittest
class TestHelperMethods(unittest.TestCase):
def test_open_to_sea(self):
# Make some test data
circle = lambda t,a,r: np.array([a+r*cos(t), a+r*sin(t)]).T
r=2
frac=0.5
t=np.arange(0,frac*2*pi+pi/10,pi/10)[::-1]
ccircle_cd=circle(t,0,r) # semicircle concave/facing south
t=np.arange(-frac*pi,frac*pi+pi/10,pi/10)[::-1]
ccircle_cl=circle(t,0,r) # a semi circle facing west,
t=np.arange(-frac*pi,frac*pi+pi/10,pi/10)
ccircle_cl2=circle(t,0,r) # a semi circle facing west, but with reversed points
t=np.arange(-frac*2*pi-pi/10,0,pi/10)[::-1]
ccircle_cnne=circle(t,0,r) # semicircle facing North North East
# test a semicircle facing south
for i in np.arange(-89,90,1):
o=howOpen(ccircle_cd,np.ones(ccircle_cd.shape[0])*i)
self.assertFalse(o>0.0,msg=str(i))
for i in np.arange(91,270,1):
o=howOpen(ccircle_cd,np.ones(ccircle_cd.shape[0])*i)
self.assertTrue(o>0.0,msg="{},{}".format(i,o))
# test a semi circle facing west
for i in np.arange(0,180,1):
o=howOpen(ccircle_cl,np.ones(ccircle_cl.shape[0])*i)
self.assertFalse(o>0.0,msg=str(i))
for i in np.arange(180,360,1):
o=howOpen(ccircle_cl,np.ones(ccircle_cl.shape[0])*i)
self.assertTrue(o>0.0,msg="{},{}".format(i,o))
# test a semi circle facing west, but with reversed points
for i in np.arange(0,180,1):
o=howOpen(ccircle_cl,np.ones(ccircle_cl2.shape[0])*i)
self.assertFalse(o>0.0,msg=str(i))
for i in np.arange(180,360,1):
o=howOpen(ccircle_cl,np.ones(ccircle_cl2.shape[0])*i)
self.assertTrue(o>0.0,msg="{},{}".format(i,o))
# test a semi circle facing North North East
o=howOpen(ccircle_cnne,np.ones(ccircle_cnne.shape[0])*10)
self.assertTrue(o>0.0,msg="{},{}".format(i,o))
# and again backwards
o=howOpen(ccircle_cnne,np.ones(ccircle_cnne.shape[0])*-190)
self.assertFalse(o>0.0,msg="{},{}".format(i,o))
if __name__=="__main__":
unittest.main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment