Is a line open to the sea
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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