diff --git a/docs/changelog.md b/docs/changelog.md index 9689b17..f5bf5be 100644 --- a/docs/changelog.md +++ b/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 diff --git a/docs/read_geotif.md b/docs/read_geotif.md index 9a05148..eea7edb 100644 --- a/docs/read_geotif.md +++ b/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) ``` diff --git a/plantcv/geospatial/read_geotif.py b/plantcv/geospatial/read_geotif.py index 14222a5..4d33190 100644 --- a/plantcv/geospatial/read_geotif.py +++ b/plantcv/geospatial/read_geotif.py @@ -62,20 +62,20 @@ 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: @@ -83,28 +83,55 @@ def read_geotif(filename, bands="R,G,B", cropto=None): 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 diff --git a/plantcv/geospatial/transform_points.py b/plantcv/geospatial/transform_points.py index fdb9b6c..afb6eee 100644 --- a/plantcv/geospatial/transform_points.py +++ b/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: diff --git a/plantcv/geospatial/transform_polygons.py b/plantcv/geospatial/transform_polygons.py index 7818f9c..3289834 100644 --- a/plantcv/geospatial/transform_polygons.py +++ b/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: diff --git a/tests/test_geospatial_read_geotif.py b/tests/test_geospatial_read_geotif.py index 7f268b8..a5e96cb 100644 --- a/tests/test_geospatial_read_geotif.py +++ b/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):