Reading TIF file

docs/changelog.md

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

docs/read_geotif.md

@@ -2,25 +2,35 @@
 
 Read in data (from tif format, most likely georeferenced image data). 
 
-**plantcv.geospatial.read_geotif**(*filename, mode="rgb"*)
+**plantcv.geospatial.read_geotif**(*filename, bands="B,G,R"*)
 
-**returns** [PlantCV Spectral_data](https://plantcv.readthedocs.io/en/latest/Spectral_data/) object instance 
+**returns** [PlantCV Spectral_data](https://plantcv.readthedocs.io/en/latest/Spectral_data/) object instance.
 
 - **Parameters:**
-    - filename - Filepath to .tif data 
-    - mode - Mode for geotif reading  
+    - 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]`)
+        - Supported symbols and their default wavelengths: 
+            - R (red) = 650nm
+            - G (green) = 560nm
+            - B (blue) = 480nm
+            - RE (rededge) = 717nm
+            - N or NIR (near infrared) = 842nm
 
 - **Example use:**
     - below
 
 
 ```python
-import plantcv.plantcv as pcv 
 import plantcv.geospatial as geo
 
 # Read geotif in
-marker = geo.read_geotif(filename="./data/example_img.tif", mode="rgb")
+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")
 
 ```
 
+![Screenshot](documentation_images/multispec_pseudo_rgb.png)
+
+![Screenshot](documentation_images/rgb.png)
+
 **Source Code:** [Here](https://github.com/danforthcenter/plantcv-geospatial/blob/main/plantcv/geospatial/read_geotif.py)

plantcv/geospatial/__init__.py

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

plantcv/geospatial/read_geotif.py

@@ -0,0 +1,124 @@
+# Read georeferenced TIF files to Spectral Image data
+
+import os
+import cv2
+import rasterio
+import numpy as np
+from plantcv.plantcv import params
+from plantcv.plantcv import fatal_error
+from plantcv.plantcv._debug import _debug
+from plantcv.plantcv.classes import Spectral_data
+
+
+def _find_closest_unsorted(array, target):
+    """Find closest index of array item with smallest distance from the target.
+
+    Parameters
+    ----------
+    array : numpy.ndarray
+        Array of wavelength labels
+    target : int, float
+        Target value
+
+    Returns
+    -------
+    int
+        Index of closest value to the target
+    """
+    return min(range(len(array)), key=lambda i: abs(array[i]-target))
+
+
+def _parse_bands(bands):
+    """Parse bands.
+
+    Parameters
+    ----------
+    bands : str
+        Comma separated string listing the order of bands
+
+    Returns
+    -------
+    list
+        List of bands
+    """
+    # Numeric list of bands
+    band_list = []
+
+    # Parse bands
+    band_strs = bands.split(",")
+
+    # Default values for symbolic bands
+    default_wavelengths = {"R": 650, "G": 560, "B": 480, "RE": 717, "N": 842, "NIR": 842}
+
+    for band in band_strs:
+        # Check if the band symbols are supported
+        if band.upper() not in default_wavelengths:
+            fatal_error(f"Currently {band} is not supported, instead provide list of wavelengths in order.")
+        # Append the default wavelength for each band
+        band_list.append(default_wavelengths[band.upper()])
+
+    return band_list
+
+
+def read_geotif(filename, bands="B,G,R"):
+    """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"
+
+    Returns
+    -------
+    plantcv.plantcv.classes.Spectral_data
+        Orthomosaic image data
+    """
+    img = rasterio.open(filename)
+    img_data = img.read()
+    img_data = img_data.transpose(1, 2, 0)  # reshape such that z-dimension is last
+    height = img.height
+    width = img.width
+    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
+
+    # Mask negative background values
+    img_data[img_data < 0.] = 0
+    # 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
+    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)
+    # Stack bands together, BGR since plot_image will convert BGR2RGB automatically
+    pseudo_rgb = cv2.merge((img_data[:, :, [id_blue]],
+                            img_data[:, :, [id_green]],
+                            img_data[:, :, [id_red]]))
+    # Gamma correction
+    if pseudo_rgb.dtype != 'uint8':
+        pseudo_rgb = pseudo_rgb.astype('float32') ** (1 / 2.2)
+        pseudo_rgb = pseudo_rgb * 255
+        pseudo_rgb = pseudo_rgb.astype('uint8')
+
+    # Make a Spectral_data instance before calculating a pseudo-rgb
+    spectral_array = Spectral_data(array_data=img_data,
+                                   max_wavelength=None,
+                                   min_wavelength=None,
+                                   max_value=np.max(img_data), min_value=np.min(img_data),
+                                   d_type=img.dtypes[0],
+                                   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"))
+
+    return spectral_array

pyproject.toml

@@ -11,6 +11,8 @@ dynamic = ["version"]
 dependencies = [
     "plantcv",
     "matplotlib",
+    "rasterio",
+
 ]
 requires-python = ">=3.6"
 authors = [

tests/conftest.py

@@ -1,4 +1,23 @@
+import os
+import pytest
 import matplotlib
 
 # Disable plotting
 matplotlib.use("Template")
+
+class TestData:
+    def __init__(self):
+        """Initialize simple variables."""
+        # Test data directory
+        self.datadir = os.path.join(os.path.dirname(os.path.abspath(__file__)), "testdata")
+        # Flat image directory
+        self.snapshot_dir = os.path.join(self.datadir, "snapshot_dir")
+        # multispectral image
+        self.cropped_tif = os.path.join(self.datadir, "615.tif")
+        # rgb image
+        self.rgb_tif = os.path.join(self.datadir, "rgb.tif")
+
+@pytest.fixture(scope="session")
+def test_data():
+    """Test data object for the main PlantCV-geospatial package."""
+    return TestData()

tests/test_geospatial_read_geotif.py

@@ -0,0 +1,25 @@
+"""Tests for geospatial.readd_geotif."""
+
+import pytest
+from plantcv.geospatial import read_geotif
+
+
+def test_geospatial_read_geotif(test_data):
+    """Test for plantcv-geospatial."""
+    # read in small 5-band tif image
+    img = read_geotif(filename=test_data.cropped_tif, bands="B,G,R,RE,N")
+    assert img.array_data.shape[2] == 5
+
+
+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
+
+
+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")