Merge branch 'main' into add_min_max_wavelength

docs/changelog.md

@@ -5,3 +5,7 @@ All notable changes to this project will be documented below.
 #### geospatial.read_geotif
 
 * v0.1dev: spectral = **geospatial.read_geotif**(*filename, bands="B,G,R"*)
+
+#### geospatial.transform_points
+
+* v0.1dev: coord = **geospatial.transform_points**(*img, geojson*)

docs/read_geotif.md

@@ -2,19 +2,20 @@
 
 Read in data (from tif format, most likely georeferenced image data). 
 
-**plantcv.geospatial.read_geotif**(*filename, bands="B,G,R"*)
+**plantcv.geospatial.read_geotif**(*filename, bands="R,G,B", cropto=None*)
 
 **returns** [PlantCV Spectral_data](https://plantcv.readthedocs.io/en/latest/Spectral_data/) object instance.
 
 - **Parameters:**
     - filename - Path of the TIF image file.
-    - bands - Comma separated string representing the order of image bands (e.g., `bands="B,G,R"`), or a list of wavelengths (e.g., `bands=[650, 560, 480]`)
+    - bands - Comma separated string representing the order of image bands (e.g., `bands="R,G,B"`), or a list of wavelengths (e.g., `bands=[650, 560, 480]`)
         - Supported symbols and their default wavelengths: 
             - R (red) = 650nm
             - G (green) = 560nm
             - B (blue) = 480nm
             - RE (rededge) = 717nm
             - N or NIR (near infrared) = 842nm
+    - cropto - A geoJSON-type shapefile to crop the input image as it is read in. Default is None. 
 
 - **Example use:**
     - below
@@ -25,7 +26,8 @@ import plantcv.geospatial as geo
 
 # Read geotif in
 ortho1 = geo.read_geotif(filename="./data/example_img.tif", bands="b,g,r,RE,NIR")
-ortho2 = geo.read_geotif(filename="./data/example_rgb_img.tif", bands="B,G,R")
+ortho2 = geo.read_geotif(filename="./data/example_rgb_img.tif", bands="R,G,B",
+                         cropto="./shapefiles/experimental_bounds.geojson)
 
 ```
 

docs/transform_points.md

@@ -0,0 +1,29 @@
+## Transform Points
+
+Transform the points from a georeferenced shapefile/GeoJSON based on an image size. 
+
+**plantcv.geospatial.transform_points**(*img, geojson*)
+
+**returns** list of transformed coordinates
+
+- **Parameters:**
+    - img - Spectral image object, likely read in with [`geo.read_geotif`](read_geotif.md)
+    - geojson - Path to the shapefile/GeoJSON containing the points. Can be Point or MultiPoint geometry.
+
+- **Context:**
+    - Transformed points can be used downstream for Python analysis, such as defining ROIs. 
+- **Example use:**
+    - below to define plot boundaries
+
+
+```python
+import plantcv.geospatial as geo
+
+# Read geotif in
+spectral = geo.read_geotif(filename="./data/example_img.tif", bands="b,g,r,RE,NIR")
+coords = geo.transform_points(img=spectral, geojson="./experimental_bounds_2024.shp")
+roi_objects = pcv.roi.custom(img=spectral.pseudo_rgb, vertices=coords)
+
+```
+
+**Source Code:** [Here](https://github.com/danforthcenter/plantcv-geospatial/blob/main/plantcv/geospatial/transform_points.py)

mkdocs.yml

@@ -18,6 +18,7 @@ nav:
   - 'PlantCV Namespace':
     - 'Geopspatial Tools':
         - Read Geo-tif Data: read_geotif.md
+        - Transform coordinate points: transform_points.md
 markdown_extensions:
   - toc:
       permalink: True

plantcv/geospatial/__init__.py

@@ -1,9 +1,11 @@
 from importlib.metadata import version
 from plantcv.geospatial.read_geotif import read_geotif
+from plantcv.geospatial.transform_points import transform_points
 
 # Auto versioning
 __version__ = version("plantcv-geospatial")
 
 __all__ = [
-    "read_geotif"
+    "read_geotif",
+    "transform_points"
     ]

plantcv/geospatial/read_geotif.py

@@ -4,14 +4,16 @@
 import cv2
 import rasterio
 import numpy as np
-from plantcv.plantcv import params
-from plantcv.plantcv import fatal_error
+import fiona
+from rasterio.mask import mask
+from plantcv.plantcv import warn, params, fatal_error
 from plantcv.plantcv._debug import _debug
 from plantcv.plantcv.classes import Spectral_data
+from shapely.geometry import shape, MultiPoint, mapping
 
 
 def _find_closest_unsorted(array, target):
-    """Find closest index of array item with smallest distance from the target.
+    """Find closest index of array item with smallest distance from
 
     Parameters
     ----------
@@ -60,38 +62,63 @@ def _parse_bands(bands):
     return band_list
 
 
-def read_geotif(filename, bands="B,G,R"):
+def read_geotif(filename, bands="R,G,B", cropto=None):
     """Read Georeferenced TIF image from file.
 
     Parameters
     ----------
     filename : str
         Path of the TIF image file.
     bands : str, list, optional
-        Comma separated string listing the order of bands or a list of wavelengths, by default "B,G,R"
+        Comma separated string listing the order of bands or a list of wavelengths, by default "R,G,B"
 
     Returns
     -------
     plantcv.plantcv.classes.Spectral_data
-        Orthomosaic image data
+        Orthomosaic image data in a Spectral_data class instance
     """
-    img = rasterio.open(filename)
-    img_data = img.read()
+    if cropto:
+        with fiona.open(cropto, 'r') as shapefile:
+            # polygon-type shapefile
+            if len(shapefile) == 1:
+                shapes = [feature['geometry'] for feature in shapefile]
+            # points-type shapefile
+            else:
+                points = [shape(feature["geometry"]) for feature in shapefile]
+                multi_point = MultiPoint(points)
+                convex_hull = multi_point.convex_hull
+                shapes = [mapping(convex_hull)]
+        # rasterio does the cropping within open
+        with rasterio.open(filename, 'r') as src:
+            img_data, geo_transform = mask(src, shapes, crop=True)
+            d_type = src.dtypes[0]
+            geo_crs = src.crs.wkt
+
+    else:
+        img = rasterio.open(filename)
+        img_data = img.read()
+        d_type = img.dtypes[0]
+        geo_transform = img.transform
+        geo_crs = img.crs.wkt
+
     img_data = img_data.transpose(1, 2, 0)  # reshape such that z-dimension is last
-    height = img.height
-    width = img.width
+    height, width, depth = img_data.shape
     wavelengths = {}
 
     # Parse bands if input is a string
     if isinstance(bands, str):
         bands = _parse_bands(bands)
-
     # Create a dictionary of wavelengths and their indices
     for i, wl in enumerate(bands):
         wavelengths[wl] = i
-
+    # Check if user input matches image dimension in z direction
+    if depth != len(bands):
+        warn(f"{depth} bands found in the image data but {filename} was provided with {bands}")
     # Mask negative background values
     img_data[img_data < 0.] = 0
+    if np.sum(img_data) == 0:
+        fatal_error(f"your image is empty, are the crop-to bounds outside of the {filename} image area?")
+
     # Make a list of wavelength keys
     wl_keys = wavelengths.keys()
     # Find which bands to use for red, green, and blue bands of the pseudo_rgb image
@@ -113,12 +140,15 @@ def read_geotif(filename, bands="B,G,R"):
                                    max_wavelength=max(wavelengths, key=wavelengths.get),
                                    min_wavelength=min(wavelengths, key=wavelengths.get),
                                    max_value=np.max(img_data), min_value=np.min(img_data),
-                                   d_type=img.dtypes[0],
+                                   d_type=d_type,
                                    wavelength_dict=wavelengths, samples=int(width),
                                    lines=int(height), interleave=None,
                                    wavelength_units="nm", array_type="datacube",
-                                   pseudo_rgb=pseudo_rgb, filename=filename, default_bands=[480, 540, 630])
-
-    _debug(visual=pseudo_rgb, filename=os.path.join(params.debug_outdir, f"{params.device}_pseudo_rgb.png"))
+                                   pseudo_rgb=pseudo_rgb, filename=filename,
+                                   default_bands=[480, 540, 630],
+                                   geo_transform=geo_transform,
+                                   geo_crs=geo_crs)
 
+    _debug(visual=pseudo_rgb, filename=os.path.join(params.debug_outdir,
+                                                    f"{params.device}_pseudo_rgb.png"))
     return spectral_array

plantcv/geospatial/transform_points.py

@@ -0,0 +1,30 @@
+# Transform georeferenced GeoJSON/shapefile points into python coordinates
+import fiona
+
+
+def transform_points(img, geojson):
+    """Takes a points-type shapefile/GeoJSON and transforms to numpy coordinates
+    Inputs:
+    img:        A spectral object from read_geotif.
+    geojson:    Path to the shape file containing the points.
+
+    Returns:
+    coord:      Transformed points as a list of numpy coordinates
+
+    :param img: [spectral object]
+    :param geojson: str
+    :return coord: list
+    """
+    geo_transform = img.geo_transform
+
+    coord = []
+    with fiona.open(geojson, 'r') as shapefile:
+        for i, _ in enumerate(shapefile):
+            if type((shapefile[i].geometry["coordinates"])) is list:
+                pixel_point = ~(geo_transform) * (shapefile[i].geometry["coordinates"][0])
+            if type((shapefile[i].geometry["coordinates"])) is tuple:
+                pixel_point = ~(geo_transform) * (shapefile[i].geometry["coordinates"])
+            rounded = (int(pixel_point[0]), int(pixel_point[1]))
+            coord.append(rounded)
+
+    return coord

pyproject.toml

@@ -12,7 +12,8 @@ dependencies = [
     "plantcv",
     "matplotlib",
     "rasterio",
-
+    "fiona",
+    "shapely",
 ]
 requires-python = ">=3.6"
 authors = [

tests/conftest.py

@@ -14,8 +14,18 @@ def __init__(self):
         self.snapshot_dir = os.path.join(self.datadir, "snapshot_dir")
         # multispectral image
         self.cropped_tif = os.path.join(self.datadir, "615.tif")
+        # empty image
+        self.empty_tif = os.path.join(self.datadir, "cropped_empty.tif")
         # rgb image
         self.rgb_tif = os.path.join(self.datadir, "rgb.tif")
+        # multiPoints shapefilex
+        self.pts_geojson = os.path.join(self.datadir, "test_pts.geojson")
+        # points shapefilex
+        self.single_pts_geojson = os.path.join(self.datadir, "single_test_pts.geojson")
+        # polygon shapefile
+        self.square_crop = os.path.join(self.datadir, "square_crop.geojson")
+        # points shapefile
+        self.point_crop = os.path.join(self.datadir, "point_crop.geojson")
 
 @pytest.fixture(scope="session")
 def test_data():

tests/test_geospatial_read_geotif.py

@@ -15,11 +15,33 @@ def test_geospatial_read_geotif_rgb(test_data):
     """Test for plantcv-geospatial."""
     # read in small 5-band tif image
     img = read_geotif(filename=test_data.rgb_tif, bands="R,G,B")
-    assert img.array_data.shape[2] == 4
+    assert img.pseudo_rgb.shape == (284, 261, 3)
 
 
 def test_geospatial_read_geotif_bad_input(test_data):
     """Test for plantcv-geospatial."""
     # read in small 5-band tif image
     with pytest.raises(RuntimeError):
         _ = read_geotif(filename=test_data.cropped_tif, bands="p")
+
+
+def test_geospatial_read_geotif_bad_crop(test_data):
+    """Test for plantcv-geospatial."""
+    # read in small 5-band tif image
+    with pytest.raises(RuntimeError):
+        _ = read_geotif(filename=test_data.empty_tif, bands="B,G,R,RE,N")
+
+
+def test_geospatial_read_geotif_polygon_crop(test_data):
+    """Test for plantcv-geospatial."""
+    # read in rgb image with a polygon-type shapefile
+    img = read_geotif(filename=test_data.rgb_tif, bands=[650, 560, 480],
+                      cropto=test_data.square_crop)
+    assert img.pseudo_rgb.shape == (80, 83, 3)
+
+
+def test_geospatial_read_geotif_point_crop(test_data):
+    """Test for plantcv-geospatial."""
+    # read in rgb image with a polygon-type shapefile
+    img = read_geotif(filename=test_data.rgb_tif, bands="R,G,B", cropto=test_data.point_crop)
+    assert img.pseudo_rgb.shape == (41, 46, 3)

tests/test_geospatial_transform_points.py

@@ -0,0 +1,19 @@
+"""Tests for geospatial.readd_geotif."""
+
+from plantcv.geospatial import read_geotif, transform_points
+
+
+def test_geospatial_transform_points(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")
+    coords = transform_points(img=img, geojson=test_data.pts_geojson)
+    assert len(coords) == 4
+
+
+def test_geospatial_transform_single_points(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")
+    coords = transform_points(img=img, geojson=test_data.single_pts_geojson)
+    assert len(coords) == 8

tests/testdata/point_crop.geojson

@@ -0,0 +1,11 @@
+{
+"type": "FeatureCollection",
+"name": "point_crop",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32615" } },
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 720198.797875224030577, 4302928.175875684246421 ] } },
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 720199.770441700122319, 4302928.186563228257 ] } },
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 720199.781129243900068, 4302927.310184645466506 ] } },
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ 720198.851312942570075, 4302927.352934820577502 ] } }
+]
+}

tests/testdata/single_test_pts.geojson

@@ -0,0 +1,15 @@
+{
+"type": "FeatureCollection",
+"name": "single_test_pts",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32615" } },
+"features": [
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ 720201.455548222991638, 4302929.009466916322708 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ 720200.939437506720424, 4302928.515762113966048 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ 720200.359096014057286, 4302927.93337285425514 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ 720199.889628155739047, 4302927.349608448334038 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ 720199.406284277676605, 4302926.797856044024229 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ 720198.922733629238792, 4302926.185171222314239 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ 720198.442067798110656, 4302925.539512793533504 ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "Point", "coordinates": [ 720197.996919097611681, 4302925.057554029859602 ] } }
+]
+}

tests/testdata/square_crop.geojson

@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"name": "square_crop",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32615" } },
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 720200.657507826690562, 4302928.218625859357417 ], [ 720199.535315738874488, 4302929.041566723957658 ], [ 720198.855781324207783, 4302928.114928886294365 ], [ 720199.977973412023857, 4302927.291988021694124 ], [ 720200.657507826690562, 4302928.218625859357417 ] ] ] } }
+]
+}

tests/testdata/test_pts.geojson

@@ -0,0 +1,11 @@
+{
+"type": "FeatureCollection",
+"name": "test_pts",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32615" } },
+"features": [
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "MultiPoint", "coordinates": [ [ 720198.48000389453955, 4302927.004775654524565 ] ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "MultiPoint", "coordinates": [ [ 720199.890186188276857, 4302928.102779876440763 ] ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "MultiPoint", "coordinates": [ [ 720200.501309029874392, 4302929.21947811357677 ] ] } },
+{ "type": "Feature", "properties": { "id": null }, "geometry": { "type": "MultiPoint", "coordinates": [ [ 720197.757872894289903, 4302926.144173440523446 ] ] } }
+]
+}