General functionality for handling vector data moved to vectorlib and extended...

General functionality for handling vector data moved to vectorlib and extended with the help of external package delaynay-triangulation

lmgeo/toolbox/vcliplib.py

-# Copyright (c) 2004-2020 WUR, Wageningen
 """VclipLib - a Python library for clipping from vector based spatial data"""
 from array import array
 from enum import Enum
-
-__author__ = "Steven B. Hoek"
+from .vectorlib import ShapeRecords, SimpleShape
 
 class HLocation(Enum):
     HBETWEEN = 1
@@ -24,26 +22,6 @@ def Locate(x, y, bbox):
     elif (y>=bbox[3]): result[1] = VLocation.ABOVE
     return result
 
-class ShapeRecords(list):
-    def __init__(self):
-        pass
-
-class SimpleShape(object):
-    shape = None
-    def __init__(self):
-        self.shape = InnerShape()
-    
-class InnerShape(object):
-    # List called parts holds index which indicates the point of next part
-    # List called points holds list / array with only 2 places
-    bbox = None
-    parts = None
-    points = None
-    def __init__(self):
-        self.bbox = 4*[0.0]
-        self.parts = []
-        self.points = []
-        
 def get_clip_result(shpRecs, bbox):
     # Get a simplified copy of each record, with its shape and the parts and points thereof
     simpleRecords = ShapeRecords()

lmgeo/toolbox/vectorlib.py

+"""VectorLib - a Python library for representing and handling vector based spatial data"""
+from math import sqrt
+from delaunay_triangulation.triangulate import delaunay
+from delaunay_triangulation import typing
+
+def distance(xy1, xy2):
+    # Be prepared to handle sequenced coordinates as well as x & y attributes
+    if hasattr(xy1, "__len__") and hasattr(xy2, "__len__"):
+        xdist, ydist = xy1[0] - xy2[0], xy1[1] - xy2[1]
+    elif hasattr(xy1, "x") and hasattr(xy1, "y") and hasattr(xy2, "x") and hasattr(xy2, "y"):
+        xdist, ydist = xy1.x - xy2.x, xy1.y - xy2.y
+    result = sqrt(xdist**2 + ydist**2)
+    return result
+
+class Triangle(object):
+    __points = []
+    
+    def __init__(self, points):
+        # Check input - be prepared to handle input in different forms
+        if hasattr(points, "__len__"):
+            if len(points) != 3: raise ValueError("Sequence of length 3 expected!")
+        else:
+            if not hasattr(points, "__getitem__"): ValueError("Indexed object expected!")
+        try:
+            for i in range(3):
+                if (not hasattr(points[i], "x")) or (not hasattr(points[i], "y")): 
+                    raise ValueError("Attributes x and y expected!")
+        except Exception as e:
+            print(e)
+        finally:
+            self.__points = points
+
+    def area(self):
+        # Calculate the length of the sides
+        a = distance(self.__points[0], self.__points[1])
+        b = distance(self.__points[1], self.__points[2])
+        c = distance(self.__points[0], self.__points[2])
+        
+        # Now we'll use Heron's formula
+        s = (a + b + c) / 2 # semi-perimeter
+        result = sqrt(s * (s - a) * (s - b) * (s - c))
+        return result
+        
+    def centroid(self):
+        # Validate input
+        if hasattr(self, "__len__"): 
+            assert len(self.__points) == 3, "Triangle must have 3 corners!" 
+            
+        # Calculate the centroid by averaging coordinates
+        xc = [pt[0] for pt in self.__points] / 3
+        yc = [pt[1] for pt in self.__points] / 3
+        return xc, yc
+
+# TODO: develop this class 
+class Polygon(object):
+    __vertices = []
+     
+    def __init__(self, vertices):
+        pass
+         
+    def area(self):
+        pass
+     
+    def contains_point(self, point):
+        # For the meantime, use method contains_points from class Path found in matplotlib.path
+        pass
+        
+    def triangles(self):
+        # Carry out a Delaunay triangulation to split up the polgon into triangles
+        pass
+        
+    def centroid(self):
+        pass
+    
+class ShapeRecords(list):
+    def __init__(self):
+        pass
+
+class SimplePoint(object):
+    x = None
+    y = None
+    
+    def __init__(self, x, y):
+        self.x = x
+        self.y = y   
+
+class SimpleShape(object):
+    shape = None
+    def __init__(self):
+        self.shape = InnerShape()
+        
+    def get_point_series(self):
+        result = []
+        for i in range(len(self.shape.parts)):
+            points = []
+            i_start = self.shape.parts[i]
+            if i==len(self.shape.parts)-1:
+                i_end = len(self.shape.points)
+            else:
+                i_end = self.shape.parts[i+1]
+            for i in self.shape.points[i_start:i_end]:
+                pt  = SimplePoint(i[0], i[1])
+                points.append(pt)
+            result.append(points)
+        return result 
+        
+    # Simple method for estimating centroid    
+    def get_dirty_centroid(self):
+        points = sum(self.get_point_series(), [])
+        x_coords = [p.x for p in points]
+        y_coords = [p.y for p in points]
+        _len = len(x_coords)
+        return [sum(x_coords)/_len, sum(y_coords)/_len]
+    
+    def get_centroid(self):
+        result = (-1.0, -1.0)
+        try:
+            # Get all the points and then do the triangulation
+            points = sum(self.get_point_series(), [])
+            triangles = delaunay([typing.Vertex(p.x, p.y) for p in points])
+    
+            # Now loop over the triangles
+            total_area = total_xweighted = total_yweighted = 0.0
+            for tr in triangles:
+                mytr = Triangle(tr)
+                total_area += mytr.area()
+                total_xweighted += mytr.centroid()[0] * mytr.area() 
+                total_yweighted += mytr.centroid()[1] * mytr.area()
+            result = (total_xweighted / total_area, total_yweighted / total_area)
+            
+        except Exception as e:
+            print(e)
+        finally:
+            return result
+        
+    def get_nearest_point(self, otherpoint):
+        points = sum(self.get_point_series(), [])
+        mydicts = [{"point":(p.x, p.y), "distance":distance((p.x, p.y), otherpoint)} for p in points]
+        minval = min([d["distance"] for d in mydicts])
+        result = list(filter(lambda d: abs(d["distance"] - minval) < 0.000000001, mydicts))[0]
+        return result
+    
+    def get_area(self):
+        # Get all the points and then do the triangulation
+        points = sum(self.get_point_series(), [])
+        triangles = delaunay([typing.Vertex(p.x, p.y) for p in points])
+
+        # Now loop over the triangles
+        result = 0.0
+        for tr in triangles:
+            mytr = Triangle(tr)
+            result += mytr.area()
+        return result
+        
+class InnerShape(object):
+    # A list called parts holds index which indicates the point of next part
+    # A list called points holds list / array with only 2 places
+    bbox = None
+    parts = None
+    points = None
+    def __init__(self):
+        self.bbox = 4*[0.0]
+        self.parts = []
+        self.points = []
+        
\ No newline at end of file

requirements.txt

@@ -5,4 +5,4 @@ tifffile==2019.3.18
 tables==3.5.2
 netCDF4==1.5.1.2
 libtiff==0.4.2
-
+delaunay-triangulation=1.0.3