Handle nodata values and masks in read_geotif

docs/changelog.md

@@ -4,7 +4,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"*)
+* v0.1dev: spectral = **geospatial.read_geotif**(*filename, bands="B,G,R", cropto=None*)
 
 #### geospatial.transform_points
 

docs/read_geotif.md

@@ -9,24 +9,29 @@ Read in data (from tif format, most likely georeferenced image data).
 - **Parameters:**
     - filename - Path of the TIF image file.
     - 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
+        - Supported symbols (case insensitive) and their default wavelengths: 
+            - "R" (red) = 650nm
+            - "G" (green) = 560nm
+            - "B" (blue) = 480nm
+            - "RE" (rededge) = 717nm
+            - "N" or "NIR" (near infrared) = 842nm
+            - "mask" (a binary mask to represent regions of no data in the image) = 0nm
     - cropto - A geoJSON-type shapefile to crop the input image as it is read in. Default is None. 
 
+- **Context:**
+    - This function aims to handle variability in data type, depth, and common "No-Data" values of Geo-tifs. There is some flexibility in formats supported but we encourage people to reach out on [GitHub](https://github.com/danforthcenter/plantcv-geospatial/issues) and collaborate with the PlantCV community to expand our support.
+    - Negative values are masked to a value of 0 to account for common no data values, and for errant negative values that can result from calibration since reflectance is bounded 0-1.
+    - Utilizing `cropto` can significantly reduce the memory needed to run a geospatial workflow. 
+
 - **Example use:**
     - below
 
-
 ```python
 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="R,G,B",
+ortho2 = geo.read_geotif(filename="./data/example_rgb_img.tif", bands="R,G,B,mask",
                          cropto="./shapefiles/experimental_bounds.geojson)
 
 ```

plantcv/geospatial/read_geotif.py

@@ -62,49 +62,76 @@ def _parse_bands(bands):
     return band_list
 
 
-def read_geotif(filename, bands="R,G,B", cropto=None):
-    """Read Georeferenced TIF image from file.
+def _read_geotif_and_shapefile(filename, cropto):
+    """Read Georeferenced TIF image from file and shapefile for cropping.
 
     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 "R,G,B"
+    cropto : str
+        Path of the shapefile to crop the image
 
     Returns
     -------
-    plantcv.plantcv.classes.Spectral_data
-        Orthomosaic image data in a Spectral_data class instance
+    tuple
+        Tuple of image data, geotransform, data type, and crs
     """
     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:
+            if len(shapefile) != 1:
                 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)
+            img_data, _ = mask(src, shapes, crop=True)
+            metadata = src.meta.copy()
             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
+        metadata = img.meta.copy()
+
+    return img_data, d_type, metadata
+
+
+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 "R,G,B"
+
+    Returns
+    -------
+    plantcv.plantcv.classes.Spectral_data
+        Orthomosaic image data in a Spectral_data class instance
+    """
+    # Read the geotif image and shapefile for cropping
+    img_data, d_type, metadata = _read_geotif_and_shapefile(filename, cropto)
 
     img_data = img_data.transpose(1, 2, 0)  # reshape such that z-dimension is last
     height, width, depth = img_data.shape
     wavelengths = {}
 
+    # Check for mask
+    mask_layer = None
+    for i in range(depth):
+        if len(np.unique(img_data[:, :, [i]])) == 2:
+            mask_layer = img_data[:, :, [i]]
+            img_data = np.delete(img_data, i, 2)
+
     # Parse bands if input is a string
     if isinstance(bands, str):
         bands = _parse_bands(bands)
@@ -114,17 +141,19 @@ def read_geotif(filename, bands="R,G,B", cropto=None):
     # 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}")
+    if depth < len(bands):
+        fatal_error("your image depth is less than the specified number of 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()
+    if mask_layer is not None:
+        img_data = np.where(mask_layer == 0, 0, img_data)
     # Find which bands to use for red, green, and blue bands of the pseudo_rgb image
-    id_red = _find_closest_unsorted(array=np.array([float(i) for i in wl_keys]), target=630)
-    id_green = _find_closest_unsorted(array=np.array([float(i) for i in wl_keys]), target=540)
-    id_blue = _find_closest_unsorted(array=np.array([float(i) for i in wl_keys]), target=480)
+    id_red = _find_closest_unsorted(array=np.array([float(i) for i in wavelengths]), target=630)
+    id_green = _find_closest_unsorted(array=np.array([float(i) for i in wavelengths]), target=540)
+    id_blue = _find_closest_unsorted(array=np.array([float(i) for i in wavelengths]), target=480)
     # Stack bands together, BGR since plot_image will convert BGR2RGB automatically
     pseudo_rgb = cv2.merge((img_data[:, :, [id_blue]],
                             img_data[:, :, [id_green]],
@@ -146,9 +175,7 @@ def read_geotif(filename, bands="R,G,B", cropto=None):
                                    wavelength_units="nm", array_type="datacube",
                                    pseudo_rgb=pseudo_rgb, filename=filename,
                                    default_bands=[480, 540, 630],
-                                   geo_transform=geo_transform,
-                                   geo_crs=geo_crs)
+                                   metadata=metadata)
 
-    _debug(visual=pseudo_rgb, filename=os.path.join(params.debug_outdir,
-                                                    f"{params.device}_pseudo_rgb.png"))
+    _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

@@ -15,7 +15,7 @@ def transform_points(img, geojson):
     :param geojson: str
     :return coord: list
     """
-    geo_transform = img.geo_transform
+    geo_transform = img.metadata["transform"]
 
     coord = []
     with fiona.open(geojson, 'r') as shapefile:

plantcv/geospatial/transform_polygons.py

@@ -15,7 +15,7 @@ def transform_polygons(img, geojson):
     :param geojson: str
     :return coord: list
     """
-    geo_transform = img.geo_transform
+    geo_transform = img.metadata["transform"]
     coord = []
     with fiona.open(geojson, 'r') as shapefile:
         for row in shapefile:

tests/test_geospatial_read_geotif.py

@@ -13,7 +13,7 @@ def test_geospatial_read_geotif(test_data):
 
 def test_geospatial_read_geotif_rgb(test_data):
     """Test for plantcv-geospatial."""
-    # read in small 5-band tif image
+    # read in small tif image
     img = read_geotif(filename=test_data.rgb_tif, bands="R,G,B")
     assert img.pseudo_rgb.shape == (284, 261, 3)
 
@@ -30,6 +30,13 @@ def test_geospatial_read_geotif_bad_crop(test_data):
     # 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_bad_bands(test_data):
+    """Test for plantcv-geospatial."""
+    # read in small 5-band tif image
+    with pytest.raises(RuntimeError):
+        _ = read_geotif(filename=test_data.rgb_tif, bands="B,G,R,RE,N")
 
 
 def test_geospatial_read_geotif_polygon_crop(test_data):