feat: adds features module with 2D spatial index and basic property utils

doc/index.rst

@@ -13,6 +13,7 @@ Imagery Transmission Format (NITF) standards.
    _apidoc/aws.osml.photogrammetry
    _apidoc/aws.osml.gdal
    _apidoc/aws.osml.image_processing
+   _apidoc/aws.osml.features
 
 
 Indices and tables

setup.cfg

@@ -43,6 +43,7 @@ install_requires =
     scikit-optimize>=0.9.0
     cachetools>=5.3.0
     geojson>=3.0.0
+    shapely>=2.0.2
     pyproj>=3.6.0
     omegaconf==2.3.0;python_version<'3.10.0'
     xsdata>=23.8

src/aws/osml/features/__init__.py

@@ -0,0 +1,8 @@
+from .feature_index import Feature2DSpatialIndex, STRFeature2DSpatialIndex
+from .imaged_feature_property_accessor import ImagedFeaturePropertyAccessor
+
+__all__ = [
+    "Feature2DSpatialIndex",
+    "ImagedFeaturePropertyAccessor",
+    "STRFeature2DSpatialIndex",
+]

src/aws/osml/features/feature_index.py

@@ -0,0 +1,67 @@
+from abc import ABC, abstractmethod
+from typing import Iterable, Optional
+
+import geojson
+import shapely
+
+from .imaged_feature_property_accessor import ImagedFeaturePropertyAccessor
+
+
+class Feature2DSpatialIndex(ABC):
+    """
+    A query-only spatial index allowing clients to lookup features using 2D geometries
+    """
+
+    @abstractmethod
+    def find_intersects(self, geometry: shapely.Geometry) -> Iterable[geojson.Feature]:
+        """
+        Return the features intersecting the input geometry.
+
+        :param geometry: geometry to query the index
+        :return: the features
+        """
+
+    @abstractmethod
+    def find_nearest(self, geometry: shapely.Geometry, max_distance: Optional[float] = None) -> Iterable[geojson.Feature]:
+        """
+        Return the nearest feature for the input geometry based on distance within two-dimensional Cartesian space.
+
+        :param geometry: geometry to query the index
+        :param max_distance: maximum distance
+        :return: the nearest features
+        """
+
+
+class STRFeature2DSpatialIndex(Feature2DSpatialIndex):
+    """
+    Implementation of the 2D spatial index for GeoJSON features using Shapely's Sort-Tile-Recursive (STR)
+    tree datastructure.
+    """
+
+    def __init__(
+        self,
+        feature_collection: geojson.FeatureCollection,
+        use_image_geometries: bool = True,
+        property_accessor: ImagedFeaturePropertyAccessor = ImagedFeaturePropertyAccessor(),
+    ) -> None:
+        self.use_image_geometries = use_image_geometries
+        self.features = feature_collection.features
+        if use_image_geometries and property_accessor is not None:
+            geometries = [property_accessor.find_image_geometry(feature) for feature in self.features]
+        else:
+            geometries = [(shapely.shape(feature.geometry), feature) for feature in self.features]
+
+        self.index = shapely.STRtree(geometries)
+
+    def find_intersects(self, geometry: shapely.Geometry) -> Iterable[geojson.Feature]:
+        result_indexes = self.index.query(geometry, predicate="intersects")
+        return [self.features[i] for i in result_indexes]
+
+    def find_nearest(self, geometry: shapely.Geometry, max_distance: Optional[float] = None) -> Iterable[geojson.Feature]:
+        if max_distance is None:
+            if self.use_image_geometries:
+                max_distance = 50
+            else:
+                max_distance = 1.0
+        result_indexes = self.index.query_nearest(geometry, max_distance=max_distance)
+        return [self.features[i] for i in result_indexes]

src/aws/osml/features/imaged_feature_property_accessor.py

@@ -0,0 +1,136 @@
+import json
+from typing import Optional
+
+import geojson
+import shapely
+
+
+class ImagedFeaturePropertyAccessor:
+    """
+    This class contains utility functions that ensure the property names / values for features derived from imagery
+    are consistently implemented. These specifications are still evolving so the intent is to encapsulate all of the
+    names in this one class so that changes do not ripple through the rest of the software baseline.
+    """
+
+    IMAGE_GEOMETRY = "imageGeometry"
+    IMAGE_BBOX = "imageBBox"
+
+    BOUNDS_IMCORDS = "bounds_imcoords"
+    GEOM_IMCOORDS = "geom_imcoords"
+    DETECTION = "detection"
+    TYPE = "type"
+    COORDINATES = "coordinates"
+    PIXEL_COORDINATES = "pixelCoordinates"
+
+    def __init__(self, allow_deprecated: bool = True):
+        """
+        Construct an instance of the property accessor with configuration options.
+
+        :param allow_deprecated: if true the accessor will work with deprecated property names.
+        """
+        self.allow_deprecated = allow_deprecated
+        pass
+
+    def find_image_geometry(self, feature: geojson.Feature) -> Optional[shapely.Geometry]:
+        """
+        This function searches through the properties of a GeoJSON feature that are known to contain the geometry
+        of the feature in image coordinates. If found an appropriate 2D shape is constructed and returned. Note that
+        this search is conducted in priority order giving preference to the current preferred "imageGeometry" and
+        "bboxGeometry" properties. If neither of those is available and the accessor has been configured to search
+        deprecated properties then the "geom_imcoords", "detection", and "bounds_imcoords" properties are searched
+        in that order.
+
+        :param feature: a GeoJSON feature that might contain an image geometry property
+        :return: a 2D shape representing the image geometry or None
+        """
+        # The "imageGeometry" property is the current preferred encoding of image geometries for these
+        # features. The format follows the same type and coordinates structure used by shapely so we can
+        # construct the geometry directly from these values.
+        if self.IMAGE_GEOMETRY in feature.properties:
+            return shapely.geometry.shape(feature.properties[self.IMAGE_GEOMETRY])
+
+        # If a full image geometry is not provided we might be able to construct a Polygon boundary from the
+        # "imageBBox" property. The property contains a [minx, miny, maxx, maxy] bounding box. If available we
+        # can construct a Polygon boundary from those 4 corners.
+        if self.IMAGE_BBOX in feature.properties:
+            bbox = feature.properties[self.IMAGE_BBOX]
+            return shapely.geometry.box(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3])
+
+        # !!!!! ALL PROPERTIES BELOW THIS LINE ARE DEPRECATED !!!!!
+        if self.allow_deprecated:
+            # The current convention for the "geom_imcoords" allows a single external ring for a Polygon boundary to be
+            # captured as a list of coordinates.
+            if self.GEOM_IMCOORDS in feature.properties:
+                return shapely.geometry.Polygon(shell=feature.properties[self.GEOM_IMCOORDS])
+
+            # Some inputs may have a "detection" property with child "type" and "pixelCoordinates" properties. If these
+            # are found we can construct the appropriate shape.
+            if self.DETECTION in feature.properties and self.PIXEL_COORDINATES in feature.properties[self.DETECTION]:
+                temp_geom = {
+                    self.TYPE: feature.properties[self.DETECTION][self.TYPE],
+                    self.COORDINATES: feature.properties[self.DETECTION][self.PIXEL_COORDINATES],
+                }
+                return shapely.geometry.shape(temp_geom)
+
+            # The current convention for "bounds_imcoords" is a [minx, miny, maxx, maxy] bounding box. If available we
+            # can construct a Polygon boundary from those 4 corners.
+            if self.BOUNDS_IMCORDS in feature.properties:
+                bbox = feature.properties[self.BOUNDS_IMCORDS]
+                return shapely.geometry.box(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3])
+
+        # All properties that might contain the image geometry are missing. This feature does not have image
+        # coordinates.
+        return None
+
+    def update_existing_image_geometries(self, feature: geojson.Feature, geometry: shapely.Geometry) -> None:
+        """
+        This function searches through the properties of a GeoJSON feature that are known to contain the geometry
+        of the feature in image coordinates. If found each property is overwritten with information from the
+        geometry provided. Note that for bounding box properties the bounds of the input geometry are used.
+
+        :param feature: a GeoJSON feature that might contain an image geometry property
+        :param geometry: the geometry to set property values for.
+        """
+        if self.IMAGE_GEOMETRY in feature.properties:
+            ImagedFeaturePropertyAccessor.set_image_geometry(feature, geometry)
+
+        if self.IMAGE_BBOX in feature.properties:
+            ImagedFeaturePropertyAccessor.set_image_bbox(feature, geometry)
+
+        # !!!!! ALL PROPERTIES BELOW THIS LINE ARE DEPRECATED !!!!!
+        if self.allow_deprecated:
+            if self.GEOM_IMCOORDS in feature.properties:
+                coordinates = shapely.geometry.mapping(geometry)[self.COORDINATES]
+                if isinstance(geometry, shapely.geometry.Polygon):
+                    feature.properties[self.GEOM_IMCOORDS] = coordinates[0]
+                else:
+                    feature.properties[self.GEOM_IMCOORDS] = coordinates
+
+            if self.DETECTION in feature.properties and self.PIXEL_COORDINATES in feature.properties[self.DETECTION]:
+                geometry_mapping = shapely.geometry.mapping(geometry)
+                feature.properties[self.DETECTION][self.TYPE] = geometry_mapping[self.TYPE]
+                feature.properties[self.DETECTION][self.PIXEL_COORDINATES] = geometry_mapping[self.COORDINATES]
+
+            if self.BOUNDS_IMCORDS in feature.properties:
+                feature.properties[self.BOUNDS_IMCORDS] = list(geometry.bounds)
+
+    @classmethod
+    def set_image_geometry(cls, feature: geojson.Feature, geometry: shapely.Geometry) -> None:
+        """
+        Add or set the "imageGeometry" property for a feature. This is a 2D geometry that supports a variety of
+        types (points, lines, polygons, etc.)
+
+        :param feature: a GeoJSON feature that will contain the property
+        :param geometry: the geometry value
+        """
+        feature.properties[cls.IMAGE_GEOMETRY] = json.loads(shapely.to_geojson(geometry))
+
+    @classmethod
+    def set_image_bbox(cls, feature: geojson.Feature, geometry: shapely.Geometry) -> None:
+        """
+        Add or set the "imageBBox" property for a feature. this is a [minx, miny, maxx, maxy] bounds for this object.
+
+        :param feature: a GeoJSON feature that will contain the property
+        :param geometry: the geometry value
+        """
+        feature.properties[cls.IMAGE_BBOX] = list(geometry.bounds)

test/aws/osml/features/test_feature_index.py

@@ -0,0 +1,29 @@
+import unittest
+
+import geojson
+import shapely
+
+
+class TestFeatureIndex(unittest.TestCase):
+    def setUp(self):
+        from aws.osml.features import STRFeature2DSpatialIndex
+
+        test_features = []
+        for r in range(0, 30, 10):
+            for c in range(0, 30, 10):
+                test_features.append(geojson.Feature(geometry=None, properties={"imageBBox": [c, r, c + 5, r + 5]}))
+        test_fc = geojson.FeatureCollection(features=test_features)
+
+        self.index = STRFeature2DSpatialIndex(test_fc, use_image_geometries=True)
+
+    def test_find_partial_intersects(self):
+        results = self.index.find_intersects(shapely.box(-1, -1, 11, 11))
+        assert len(list(results)) == 4
+
+    def test_find_contains(self):
+        results = self.index.find_intersects(shapely.box(-1, -1, 31, 31))
+        assert len(list(results)) == 9
+
+    def test_find_nearest(self):
+        results = self.index.find_nearest(shapely.Point(1, 1), max_distance=5)
+        assert len(list(results)) == 1