#!/usr/bin/env python
"""Polygon manipulations
"""
from utilities.numerical_tools import ensure_numeric
def point_on_line(x, y, x0, y0, x1, y1):
"""Determine whether a point is on a line segment
Input: x, y, x0, x0, x1, y1: where
point is given by x, y
line is given by (x0, y0) and (x1, y1)
"""
from Numeric import array, dot, allclose
from math import sqrt
a = array([x - x0, y - y0])
a_normal = array([a[1], -a[0]])
b = array([x1 - x0, y1 - y0])
if dot(a_normal, b) == 0:
#Point is somewhere on the infinite extension of the line
len_a = sqrt(sum(a**2))
len_b = sqrt(sum(b**2))
if dot(a, b) >= 0 and len_a <= len_b:
return True
else:
return False
else:
return False
def inside_polygon(points, polygon, closed = True, verbose = False):
"""Determine points inside a polygon
Functions inside_polygon and outside_polygon have been defined in
terms af separate_by_polygon which will put all inside indices in
the first part of the indices array and outside indices in the last
See separate_points_by_polygon for documentation
"""
from Numeric import array, Float, reshape
if verbose: print 'Checking input to inside_polygon'
try:
points = ensure_numeric(points, Float)
except:
msg = 'Points could not be converted to Numeric array'
raise msg
try:
polygon = ensure_numeric(polygon, Float)
except:
msg = 'Polygon could not be converted to Numeric array'
raise msg
if len(points.shape) == 1:
one_point = True
points = reshape(points, (1,2))
else:
one_point = False
indices, count = separate_points_by_polygon(points, polygon,
closed, verbose)
if one_point:
return count == 1
else:
return indices[:count]
def outside_polygon(points, polygon, closed = True, verbose = False):
"""Determine points outside a polygon
Functions inside_polygon and outside_polygon have been defined in
terms af separate_by_polygon which will put all inside indices in
the first part of the indices array and outside indices in the last
See separate_points_by_polygon for documentation
"""
from Numeric import array, Float, reshape
if verbose: print 'Checking input to outside_polygon'
try:
points = ensure_numeric(points, Float)
except:
msg = 'Points could not be converted to Numeric array'
raise msg
try:
polygon = ensure_numeric(polygon, Float)
except:
msg = 'Polygon could not be converted to Numeric array'
raise msg
if len(points.shape) == 1:
one_point = True
points = reshape(points, (1,2))
else:
one_point = False
indices, count = separate_points_by_polygon(points, polygon,
closed, verbose)
if one_point:
return count != 1
else:
if count == len(indices):
# No points are outside
return []
else:
return indices[count:][::-1] #return reversed
def separate_points_by_polygon(points, polygon,
closed = True, verbose = False):
"""Determine whether points are inside or outside a polygon
Input:
points - Tuple of (x, y) coordinates, or list of tuples
polygon - list of vertices of polygon
closed - (optional) determine whether points on boundary should be
regarded as belonging to the polygon (closed = True)
or not (closed = False)
Outputs:
indices: array of same length as points with indices of points falling
inside the polygon listed from the beginning and indices of points
falling outside listed from the end.
count: count of points falling inside the polygon
The indices of points inside are obtained as indices[:count]
The indices of points outside are obtained as indices[count:]
Examples:
U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U)
will return the indices [0, 2, 1] and count == 2 as only the first
and the last point are inside the unit square
Remarks:
The vertices may be listed clockwise or counterclockwise and
the first point may optionally be repeated.
Polygons do not need to be convex.
Polygons can have holes in them and points inside a hole is
regarded as being outside the polygon.
Algorithm is based on work by Darel Finley,
http://www.alienryderflex.com/polygon/
"""
from Numeric import array, Float, reshape, Int, zeros
#Input checks
try:
points = ensure_numeric(points, Float)
except:
msg = 'Points could not be converted to Numeric array'
raise msg
try:
polygon = ensure_numeric(polygon, Float)
except:
msg = 'Polygon could not be converted to Numeric array'
raise msg
assert len(polygon.shape) == 2,\
'Polygon array must be a 2d array of vertices'
assert polygon.shape[1] == 2,\
'Polygon array must have two columns'
assert len(points.shape) == 2,\
'Points array must be a 2d array'
assert points.shape[1] == 2,\
'Points array must have two columns'
N = polygon.shape[0] #Number of vertices in polygon
M = points.shape[0] #Number of points
px = polygon[:,0]
py = polygon[:,1]
#Used for an optimisation when points are far away from polygon
minpx = min(px); maxpx = max(px)
minpy = min(py); maxpy = max(py)
#Begin main loop
indices = zeros(M, Int)
inside_index = 0 #Keep track of points inside
outside_index = M-1 #Keep track of points outside (starting from end)
for k in range(M):
if verbose:
if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
#for each point
x = points[k, 0]
y = points[k, 1]
inside = False
if not x > maxpx or x < minpx or y > maxpy or y < minpy:
#Check polygon
for i in range(N):
j = (i+1)%N
#Check for case where point is contained in line segment
if point_on_line(x, y, px[i], py[i], px[j], py[j]):
if closed:
inside = True
else:
inside = False
break
else:
#Check if truly inside polygon
if py[i] < y and py[j] >= y or\
py[j] < y and py[i] >= y:
if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
inside = not inside
if inside:
indices[inside_index] = k
inside_index += 1
else:
indices[outside_index] = k
outside_index -= 1
return indices, inside_index
def separate_points_by_polygon_c(points, polygon,
closed = True, verbose = False):
"""Determine whether points are inside or outside a polygon
C-wrapper
"""
from Numeric import array, Float, reshape, zeros, Int
if verbose: print 'Checking input to separate_points_by_polygon'
#Input checks
try:
points = ensure_numeric(points, Float)
except:
msg = 'Points could not be converted to Numeric array'
raise msg
#if verbose: print 'Checking input to separate_points_by_polygon 2'
try:
polygon = ensure_numeric(polygon, Float)
except:
msg = 'Polygon could not be converted to Numeric array'
raise msg
if verbose: print 'check'
assert len(polygon.shape) == 2,\
'Polygon array must be a 2d array of vertices'
assert polygon.shape[1] == 2,\
'Polygon array must have two columns'
assert len(points.shape) == 2,\
'Points array must be a 2d array'
assert points.shape[1] == 2,\
'Points array must have two columns'
N = polygon.shape[0] #Number of vertices in polygon
M = points.shape[0] #Number of points
from polygon_ext import separate_points_by_polygon
if verbose: print 'Allocating array for indices'
indices = zeros( M, Int )
if verbose: print 'Calling C-version of inside poly'
count = separate_points_by_polygon(points, polygon, indices,
int(closed), int(verbose))
return indices, count
class Polygon_function:
"""Create callable object f: x,y -> z, where a,y,z are vectors and
where f will return different values depending on whether x,y belongs
to specified polygons.
To instantiate:
Polygon_function(polygons)
where polygons is a list of tuples of the form
[ (P0, v0), (P1, v1), ...]
with Pi being lists of vertices defining polygons and vi either
constants or functions of x,y to be applied to points with the polygon.
The function takes an optional argument, default which is the value
(or function) to used for points not belonging to any polygon.
For example:
Polygon_function(polygons, default = 0.03)
If omitted the default value will be 0.0
Note: If two polygons overlap, the one last in the list takes precedence
FIXME? : Currently, coordinates specified here are assumed to be relative to the origin (georeference) used by domain. This function is more general than domain so perhaps it'd be an idea to allow the optional argument georeference???
"""
def __init__(self, regions, default = 0.0):
try:
len(regions)
except:
msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
raise msg
T = regions[0]
try:
a = len(T)
except:
msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons
raise msg
assert a == 2, 'Must have two component each: %s' %T
self.regions = regions
self.default = default
def __call__(self, x, y):
from Numeric import ones, Float, concatenate, array, reshape, choose
x = array(x).astype(Float)
y = array(y).astype(Float)
N = len(x)
assert len(y) == N
points = concatenate( (reshape(x, (N, 1)),
reshape(y, (N, 1))), axis=1 )
if callable(self.default):
z = self.default(x,y)
else:
z = ones(N, Float) * self.default
for polygon, value in self.regions:
indices = inside_polygon(points, polygon)
#FIXME: This needs to be vectorised
if callable(value):
for i in indices:
xx = array([x[i]])
yy = array([y[i]])
z[i] = value(xx, yy)[0]
else:
for i in indices:
z[i] = value
return z
def read_polygon(filename):
"""Read points assumed to form a polygon
There must be exactly two numbers in each line
"""
#Get polygon
fid = open(filename)
lines = fid.readlines()
fid.close()
polygon = []
for line in lines:
fields = line.split(',')
polygon.append( [float(fields[0]), float(fields[1])] )
return polygon
def populate_polygon(polygon, number_of_points, seed = None, exclude = None):
"""Populate given polygon with uniformly distributed points.
Input:
polygon - list of vertices of polygon
number_of_points - (optional) number of points
seed - seed for random number generator (default=None)
exclude - list of polygons (inside main polygon) from where points should be excluded
Output:
points - list of points inside polygon
Examples:
populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 )
will return five randomly selected points inside the unit square
"""
from random import uniform, seed as seed_function
seed_function(seed)
points = []
#Find outer extent of polygon
max_x = min_x = polygon[0][0]
max_y = min_y = polygon[0][1]
for point in polygon[1:]:
x = point[0]
if x > max_x: max_x = x
if x < min_x: min_x = x
y = point[1]
if y > max_y: max_y = y
if y < min_y: min_y = y
while len(points) < number_of_points:
x = uniform(min_x, max_x)
y = uniform(min_y, max_y)
append = False
if inside_polygon( [x,y], polygon ):
append = True
#Check exclusions
if exclude is not None:
for ex_poly in exclude:
if inside_polygon( [x,y], ex_poly ):
append = False
if append is True:
points.append([x,y])
return points
def point_in_polygon(polygon, delta=1e-8):
"""Return a point inside a given polygon which will be close to the
polygon edge.
Input:
polygon - list of vertices of polygon
delta - the square root of 2 * delta is the maximum distance from the
polygon points and the returned point.
Output:
points - a point inside polygon
searches in all diagonals and up and down (not left and right)
"""
import exceptions
class Found(exceptions.Exception): pass
point_in = False
while not point_in:
try:
for poly_point in polygon: #[1:]:
for x_mult in range (-1,2):
for y_mult in range (-1,2):
x = poly_point[0]
y = poly_point[1]
if x == 0:
x_delta = x_mult*delta
else:
x_delta = x+x_mult*x*delta
if y == 0:
y_delta = y_mult*delta
else:
y_delta = y+y_mult*y*delta
point = [x_delta, y_delta]
#print "point",point
if inside_polygon( point, polygon, closed = False ):
raise Found
except Found:
point_in = True
else:
delta = delta*0.1
return point
##############################################
#Initialise module
import compile
if compile.can_use_C_extension('polygon_ext.c'):
from polygon_ext import point_on_line
separate_points_by_polygon = separate_points_by_polygon_c
if __name__ == "__main__":
pass