Merge pull request #21 from danforthcenter/create_circular_geo_rois

Create circular geo rois

docs/changelog.md

@@ -2,6 +2,10 @@
 
 All notable changes to this project will be documented below.
 
+#### geospatial.points2roi_circle
+
+* v0.1dev: rois = **geospatial.points2roi_circle**(*img, geojson, radius*)
+
 #### geospatial.points_to_geojson
 
 * v0.1dev: **geosptial.points_to_geojson**(*img, viewer, out_path*)

docs/points2roi_circle.md

@@ -0,0 +1,33 @@
+## Create Circular ROIs from georeferenced points
+
+Transform the points from a Points-type georeferenced shapefile/GeoJSON into circular Regions of Interest (ROIs). 
+
+**plantcv.geospatial.points2roi_circle**(*img, geojson, radius*)
+
+**returns** list of ROIs (`plantcv.Objects` instance)
+
+- **Parameters:**
+    - img - Spectral image object, likely read in with [`geo.read_geotif`](read_geotif.md)
+    - geojson - Path to the shapefile/GeoJSON containing the points.
+    - radius - Radius of circular ROIs to get created,
+                in units matching the coordinate system of the image
+
+- **Context:**
+    - Directly create ROIs with a consistent georeferenced radius
+- **Example use:**
+    - below
+
+
+```python
+import plantcv.geospatial as geo
+import plantcv.plantcv as pcv
+
+# Read geotif in
+spectral = geo.read_geotif(filename="./data/example_img.tif", bands="b,g,r,RE,NIR")
+rois = geo.points2roi_circle(img=spectral, geojson="./points_example.geojson", radius=1)
+# Segmentation steps here
+labeled_mask = pcv.roi.filter(mask=vegetation_mask, roi=rois, roi_type="partial")
+
+```
+
+**Source Code:** [Here](https://github.com/danforthcenter/plantcv-geospatial/blob/main/plantcv/geospatial/points2roi.py)

mkdocs.yml

@@ -17,6 +17,7 @@ nav:
     - 'Adding/Editing Documentation': documentation.md
   - 'PlantCV Namespace':
     - 'Geopspatial Tools':
+        - Points to ROI (circular): points2roi_circle.md
         - Read Geo-tif Data: read_geotif.md
         - Transform coordinate points: transform_points.md
         - Transform coordinate polygons: transform_polygons.md

plantcv/geospatial/__init__.py

@@ -3,6 +3,7 @@
 from plantcv.geospatial.transform_points import transform_points
 from plantcv.geospatial.transform_polygons import transform_polygons
 from plantcv.geospatial.points_to_geojson import points_to_geojson
+from plantcv.geospatial.points2roi import points2roi_circle
 
 # Auto versioning
 __version__ = version("plantcv-geospatial")
@@ -11,5 +12,6 @@
     "read_geotif",
     "transform_points",
     "transform_polygons",
-    "points_to_geojson"
-    ]
+    "points_to_geojson",
+    "points2roi_circle"
+]

plantcv/geospatial/_helpers.py

@@ -0,0 +1,27 @@
+# PlantCV-geospatial helper functions
+import geopandas
+
+
+def _transform_geojson_crs(img, geojson):
+    """
+    Helper function for opening and transforming Coordinate System
+    of a geojson/shapefile
+
+    Keyword inputs:
+    Inputs:
+    img:        A spectral object from read_geotif.
+    geojson:    Path to the shapefile.
+
+    :param img: [spectral object]
+    :return geojson: str
+    :return gdf: geopandas
+    """
+    gdf = geopandas.read_file(geojson)
+
+    img_crs = img.metadata['crs']
+
+    # Check spectral object and geojson have the same CRS, if not then convert
+    if not gdf.crs == img_crs:
+        gdf = gdf.to_crs(crs=img_crs)
+
+    return gdf

plantcv/geospatial/points2roi.py

@@ -0,0 +1,56 @@
+# Transform georeferenced GeoJSON/shapefile points into python coordinates
+import numpy as np
+from plantcv.geospatial.transform_polygons import transform_polygons
+from plantcv.geospatial._helpers import _transform_geojson_crs
+from plantcv.plantcv import Objects
+
+
+def points2roi_circle(img, geojson, radius):
+    """Takes a points-type shapefile/GeoJSON and transforms circular ROIs,
+    saves these out to a new geoJSON file and creates ROI Objects instances
+    Inputs:
+    img:        A spectral object from read_geotif.
+    geojson:    Path to the shape file containing the points.
+    radius:     Radius of circular ROIs to get created,
+                in units matching the coordinate system (CRS) of the image
+                e.g. meters
+
+    Returns:
+    rois:       List of circular ROIs (plantcv Objects class instances)
+
+    :param img: [spectral object]
+    :param geojson: str
+    :param radius: float
+    :return rois: list
+    """
+    gdf = _transform_geojson_crs(img=img, geojson=geojson)
+
+    gdf['geometry'] = gdf.geometry.buffer(radius)
+
+    buffered_geojson = geojson + '_circles.geojson'
+    gdf.to_file(buffered_geojson, driver='GeoJSON')
+
+    geo_rois = transform_polygons(img=img, geojson=buffered_geojson)
+
+    return _points2roi(geo_rois)
+
+
+def _points2roi(roi_list):
+    """
+    Helper that takes ROI contour coordinates and populates
+    a plantcv Objects class instance
+
+    Inputs:
+    roi_list  = List of ROI contours from georeferenced origin
+
+    Returns:
+    rois    = grouped contours list
+
+    :param roi_list: list
+    :return rois: plantcv.plantcv.classes.Objects
+    """
+    rois = Objects()
+    for roi in roi_list:
+        rois.append(contour=[np.array(roi)], h=[])
+
+    return rois

pyproject.toml

@@ -14,6 +14,7 @@ dependencies = [
     "rasterio",
     "fiona",
     "shapely",
+    "geopandas",
     "geojson",
     "napari",
 ]

tests/conftest.py

@@ -28,6 +28,8 @@ def __init__(self):
         self.point_crop = os.path.join(self.datadir, "point_crop.geojson")
         # multi polygon shapefile
         self.multipolygon = os.path.join(self.datadir, "multipolygon_fortests.geojson")
+        # epsg4326 points shapefile
+        self.epsg4326_geojson = os.path.join(self.datadir, "epsg4326points.geojson")
 
 @pytest.fixture(scope="session")
 def test_data():

tests/test_geospatial_helpers.py

@@ -0,0 +1,11 @@
+"""Tests for geospatial._helpers"""
+
+from plantcv.geospatial._helpers import _transform_geojson_crs
+from plantcv.geospatial import read_geotif
+
+def test_geospatial_helpers_transform_geojson_crs(test_data):
+    """Test for plantcv-geospatial."""
+    # read in small 5-band tif image
+    img = read_geotif(filename=test_data.rgb_tif, bands="B,G,R")
+    gdf = _transform_geojson_crs(img=img, geojson=test_data.epsg4326_geojson)
+    assert gdf.crs == img.metadata['crs']

tests/test_geospatial_points2roicircle.py

@@ -0,0 +1,12 @@
+"""Tests for geospatial.points2roi_circle"""
+
+from plantcv.geospatial import read_geotif, points2roi_circle
+import numpy as np
+
+
+def test_geospatial_points2roi_circle(test_data):
+    """Test for plantcv-geospatial."""
+    # read in small 3-band tif image
+    img = read_geotif(filename=test_data.rgb_tif, bands="B,G,R")
+    rois = points2roi_circle(img=img, geojson=test_data.pts_geojson, radius=0.5)
+    assert np.all(rois.contours[0][0][0] == np.array([119, 170]))

tests/testdata/epsg4326points.geojson

@@ -0,0 +1,12 @@
+{
+"type": "FeatureCollection",
+"name": "epsg4326points",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
+"features": [
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ -90.457889274175201, 38.847795218128709 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ -90.457842139564363, 38.847839210432156 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ -90.457876704945647, 38.847903627733636 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ -90.458036962622501, 38.847880060428217 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ -90.457790291492444, 38.847996325801624 ] } }
+]
+}