diff --git a/docs/changelog.md b/docs/changelog.md index 820a7f1..57a94c4 100644 --- a/docs/changelog.md +++ b/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"*) diff --git a/docs/documentation_images/multispec_pseudo_rgb.png b/docs/documentation_images/multispec_pseudo_rgb.png new file mode 100644 index 0000000..595c903 Binary files /dev/null and b/docs/documentation_images/multispec_pseudo_rgb.png differ diff --git a/docs/documentation_images/rgb.png b/docs/documentation_images/rgb.png new file mode 100644 index 0000000..4c53520 Binary files /dev/null and b/docs/documentation_images/rgb.png differ diff --git a/docs/read_geotif.md b/docs/read_geotif.md index d49898a..5f9378f 100644 --- a/docs/read_geotif.md +++ b/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) diff --git a/plantcv/geospatial/__init__.py b/plantcv/geospatial/__init__.py index 0c878c0..db53ac4 100644 --- a/plantcv/geospatial/__init__.py +++ b/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" + ] diff --git a/plantcv/geospatial/read_geotif.py b/plantcv/geospatial/read_geotif.py new file mode 100644 index 0000000..1e63de4 --- /dev/null +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 06d1726..0d9f4eb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,6 +11,8 @@ dynamic = ["version"] dependencies = [ "plantcv", "matplotlib", + "rasterio", + ] requires-python = ">=3.6" authors = [ diff --git a/tests/conftest.py b/tests/conftest.py index 16b3a25..eb6fd03 100644 --- a/tests/conftest.py +++ b/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() diff --git a/tests/test_geospatial_read_geotif.py b/tests/test_geospatial_read_geotif.py new file mode 100644 index 0000000..ad1d698 --- /dev/null +++ b/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") diff --git a/tests/testdata/615.tif b/tests/testdata/615.tif new file mode 100644 index 0000000..f3d9dba Binary files /dev/null and b/tests/testdata/615.tif differ diff --git a/tests/testdata/rgb.tif b/tests/testdata/rgb.tif new file mode 100644 index 0000000..1846615 Binary files /dev/null and b/tests/testdata/rgb.tif differ