Merge pull request #3 from marcbrittain/dev/v0.0.4

update main to v0.0.4

Author: marcbrittain

Examples/Example - GeoLineStrings.ipynb

@@ -77,8 +77,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[[ 3.  3.  3.]\n",
-      " [ 6. 10.  2.]]\n"
+      "[[3. 3. 3.]]\n"
      ]
     }
    ],
@@ -365,8 +364,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[[2.30084261e+06 1.35167639e+05 1.00000000e+02]\n",
-      " [2.31208049e+06 7.91798511e+04 1.00000000e+02]]\n"
+      "[[2.30084261e+06 1.35167639e+05 1.00000000e+02]]\n"
      ]
     }
    ],
@@ -384,8 +382,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[[-71.2  42.5 100. ]\n",
-      " [-71.3  42.  100. ]]\n"
+      "[[-71.2  42.5 100. ]]\n"
      ]
     }
    ],
@@ -449,70 +446,6 @@
    "source": [
     "gls3.geolinestring_length()"
    ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "4c5eb5a0",
-   "metadata": {},
-   "source": [
-    "# 4. Known Bugs\n",
-    "\n",
-    "There is an edge case where if there are multiple coordinates with (x,y) overlap, but different z values, an intersection point may be ignored. In most applications, this may not arise, but in general needs to be addressed and is being actively worked on. The below code demonstrates the bug."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 18,
-   "id": "9675a260",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Intersection point at [1,1,0] and [4,3,0]\n",
-    "# There are multiple (x,y) overlap between the coords:\n",
-    "#   1. [1,1,0]\n",
-    "#   2. [1,1,5]\n",
-    "coords5 = [[1,1,0],[4,3,0]]\n",
-    "coords6 = [[1,1,5],[1,1,0],[4,3,0]]\n",
-    "\n",
-    "gls5 = GeoLineString(coords5,is_xy=True)\n",
-    "gls6 = GeoLineString(coords6,is_xy=True)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 19,
-   "id": "7bf7d05f",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "True\n"
-     ]
-    }
-   ],
-   "source": [
-    "print(gls5.intersects(gls6)) # this is okay "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 20,
-   "id": "9dbb69e6",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "[[4. 3. 0.]]\n"
-     ]
-    }
-   ],
-   "source": [
-    "print(gls5.intersection(gls6)) # this is not correct, [1,1,0] is not identified"
-   ]
   }
  ],
  "metadata": {

README.md

@@ -1,6 +1,6 @@
 ![downloads](https://img.shields.io/pypi/dm/pygeoshape) ![license](https://img.shields.io/pypi/l/pygeoshape?color=blue) ![pypi](https://img.shields.io/pypi/v/pygeoshape)
 
-# PyGeoShape v0.0.3
+# PyGeoShape v0.0.4
 3D extension to [shapely](https://github.com/shapely/shapely/tree/main) and [pyproj](https://github.com/pyproj4/pyproj) to make working with geospatial/trajectory data easier in python. If you found PyGeoShape helpful, please consider adding a star to this repository. Thanks!
 
 ## Getting Started

pygeoshape/geolinestring.py

@@ -2,7 +2,9 @@
 from shapely.geometry import LineString
 from pyproj import Transformer
 import numpy as np
+import numba as nb
 import matplotlib.pyplot as plt
+from pygeoshape.utils import fast_intersection_append
 
 
 class GeoLineString:
@@ -45,18 +47,18 @@ def __init__(
 
         self.np_coords = np.array(self.xy_coords)
         self.length = self.geolinestring_length()
-        dimensions = len(self.xy_coords[0])
+        self.dimensions = len(self.xy_coords[0])
 
         # need linestring for each dimension
         # x, y, [z]
         # x, z, [y]
         # y, z, [x]
 
-        if dimensions == 2:
+        if self.dimensions == 2:
             self.xy = LineString([(coord[0], coord[1])
                                  for coord in self.xy_coords])
 
-        if dimensions == 3:
+        if self.dimensions == 3:
 
             self.xy = LineString(
                 [(coord[0], coord[1], coord[2]) for coord in self.xy_coords]
@@ -144,7 +146,7 @@ def intersects(self, geo_obj):
 
         Args:
         ----
-        geo_obj (GeoLineString): Second GeoLineString for comparision
+        geo_obj (GeoLineString): Second GeoLineString for comparison
 
         Returns:
         -------
@@ -164,77 +166,63 @@ def intersection(self, geo_obj, lonlat=False):
 
         Args:
         ----
-        geo_obj (GeoLineString): Second GeoLineString for comparision
+        geo_obj (GeoLineString): Second GeoLineString for comparison
         lonlat (bool): Format of output coordinates. True returns the coordinates in longitude, latitude. False returns coordinates in x, y
 
         Returns:
         -------
         intersection (List): Intersection coordinates
 
         """
-        xy_check = self.xy.intersects(geo_obj.xy)
-        xz_check = self.xz.intersects(geo_obj.xz)
-        yz_check = self.yz.intersects(geo_obj.yz)
 
-        intersection = []
+        intersection = nb.typed.List.empty_list(
+            nb.types.UniTuple(nb.float64, self.dimensions))
 
-        if all([xy_check, xz_check, yz_check]):
+        inter1 = self.xy.intersection(geo_obj.xy)
+        inter2 = self.xz.intersection(geo_obj.xz)
+        inter3 = self.yz.intersection(geo_obj.yz)
 
-            inter1 = self.xy.intersection(geo_obj.xy)
-            inter2 = self.xz.intersection(geo_obj.xz)
-            inter3 = self.yz.intersection(geo_obj.yz)
+        if not hasattr(inter1, "geoms"):
+            xy = np.array(inter1.coords)[0]
+            n_xy_intersections = 1
 
-            if not hasattr(inter1, "geoms"):
-                xy = list(inter1.coords)[0]
-                n_xy_intersections = 1
-            else:
-                n_xy_intersections = len(inter1.geoms)
-
-            if not hasattr(inter2, "geoms"):
-                xz = list(inter2.coords)[0]
-                n_xz_intersections = 1
-            else:
-                n_xz_intersections = len(inter2.geoms)
-
-            if not hasattr(inter3, "geoms"):
-                yz = list(inter3.coords)[0]
-                n_yz_intersections = 1
-            else:
-                n_yz_intersections = len(inter3.geoms)
-
-            for i in range(n_xy_intersections):
-
-                for j in range(n_xz_intersections):
-
-                    for k in range(n_yz_intersections):
+        else:
+            n_xy_intersections = len(inter1.geoms)
 
-                        if hasattr(inter1, "geoms"):
-                            xy = list(inter1.geoms[i].coords)
-                        else:
-                            xy = list(inter1.coords)
+        if not hasattr(inter2, "geoms"):
+            xz = np.array(inter2.coords)[0]
+            n_xz_intersections = 1
+        else:
+            n_xz_intersections = len(inter2.geoms)
 
-                        if hasattr(inter2, "geoms"):
-                            xz = list(inter2.geoms[j].coords)
-                        else:
-                            xz = list(inter2.coords)
+        if not hasattr(inter3, "geoms"):
+            yz = np.array(inter3.coords)[0]
+            n_yz_intersections = 1
+        else:
+            n_yz_intersections = len(inter3.geoms)
 
-                        if hasattr(inter3, "geoms"):
-                            yz = list(inter3.geoms[k].coords)
-                        else:
-                            yz = list(inter3.coords)
+        for i in range(n_xy_intersections):
 
-                        for ii in range(len(xy)):
+            for j in range(n_xz_intersections):
 
-                            x_y_z = xy[ii]
-                            x_z_y = (x_y_z[0], x_y_z[2], x_y_z[1])
+                for k in range(n_yz_intersections):
 
-                            if x_z_y in xz:
+                    if hasattr(inter1, "geoms"):
+                        xy = np.array(inter1.geoms[i].coords)
+                    else:
+                        xy = np.array(inter1.coords)
 
-                                y_z_x = (x_y_z[1], x_y_z[2], x_y_z[0])
+                    if hasattr(inter2, "geoms"):
+                        xz = np.array(inter2.geoms[j].coords)
+                    else:
+                        xz = np.array(inter2.coords)
 
-                                if y_z_x in yz:
+                    if hasattr(inter3, "geoms"):
+                        yz = np.array(inter3.geoms[k].coords)
+                    else:
+                        yz = np.array(inter3.coords)
 
-                                    intersection.append(x_y_z)
+                    fast_intersection_append(xy, xz, yz, intersection)
 
             intersection = np.unique(intersection, axis=0)