diff --git a/docs/changelog.md b/docs/changelog.md index 57a94c4..042ae9e 100644 --- a/docs/changelog.md +++ b/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*) diff --git a/docs/read_geotif.md b/docs/read_geotif.md index 5f9378f..9a05148 100644 --- a/docs/read_geotif.md +++ b/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) ``` diff --git a/docs/transform_points.md b/docs/transform_points.md new file mode 100644 index 0000000..433a9d4 --- /dev/null +++ b/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) diff --git a/mkdocs.yml b/mkdocs.yml index c69c28c..859a6ce 100644 --- a/mkdocs.yml +++ b/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 diff --git a/plantcv/geospatial/__init__.py b/plantcv/geospatial/__init__.py index db53ac4..02b338f 100644 --- a/plantcv/geospatial/__init__.py +++ b/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" ] diff --git a/plantcv/geospatial/read_geotif.py b/plantcv/geospatial/read_geotif.py index 6164d75..14222a5 100644 --- a/plantcv/geospatial/read_geotif.py +++ b/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,7 +62,7 @@ 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 @@ -68,30 +70,55 @@ def read_geotif(filename, bands="B,G,R"): 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 diff --git a/plantcv/geospatial/transform_points.py b/plantcv/geospatial/transform_points.py new file mode 100644 index 0000000..fdb9b6c --- /dev/null +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 0d9f4eb..7617ecb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,8 @@ dependencies = [ "plantcv", "matplotlib", "rasterio", - + "fiona", + "shapely", ] requires-python = ">=3.6" authors = [ diff --git a/tests/conftest.py b/tests/conftest.py index eb6fd03..f2c8e88 100644 --- a/tests/conftest.py +++ b/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(): diff --git a/tests/test_geospatial_read_geotif.py b/tests/test_geospatial_read_geotif.py index ad1d698..7f268b8 100644 --- a/tests/test_geospatial_read_geotif.py +++ b/tests/test_geospatial_read_geotif.py @@ -15,7 +15,7 @@ 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): @@ -23,3 +23,25 @@ def test_geospatial_read_geotif_bad_input(test_data): # 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) diff --git a/tests/test_geospatial_transform_points.py b/tests/test_geospatial_transform_points.py new file mode 100644 index 0000000..cd46133 --- /dev/null +++ b/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 diff --git a/tests/testdata/cropped_empty.tif b/tests/testdata/cropped_empty.tif new file mode 100644 index 0000000..85e00ce Binary files /dev/null and b/tests/testdata/cropped_empty.tif differ diff --git a/tests/testdata/point_crop.geojson b/tests/testdata/point_crop.geojson new file mode 100644 index 0000000..53f469c --- /dev/null +++ b/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 ] } } +] +} diff --git a/tests/testdata/single_test_pts.geojson b/tests/testdata/single_test_pts.geojson new file mode 100644 index 0000000..0eb6509 --- /dev/null +++ b/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 ] } } +] +} diff --git a/tests/testdata/square_crop.geojson b/tests/testdata/square_crop.geojson new file mode 100644 index 0000000..82b8505 --- /dev/null +++ b/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 ] ] ] } } +] +} diff --git a/tests/testdata/test_pts.geojson b/tests/testdata/test_pts.geojson new file mode 100644 index 0000000..3385446 --- /dev/null +++ b/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 ] ] } } +] +}