From f6798da2a3b574b29e86befa1148afd7a733b61b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E6=B5=A9=E7=80=9A=E7=8C=AB?= Date: Mon, 21 Oct 2024 17:42:39 +0900 Subject: [PATCH] rough support multi-spectral geotiff Partly support #115 --- easyidp/cvtools.py | 165 +++++++++++++++++++++++------------------- easyidp/data.py | 13 ++++ easyidp/geotiff.py | 40 ++++++++-- tests/test_cvtools.py | 99 +++---------------------- tests/test_geotiff.py | 37 +++++++++- 5 files changed, 183 insertions(+), 171 deletions(-) diff --git a/easyidp/cvtools.py b/easyidp/cvtools.py index d32b811..ea9cd98 100644 --- a/easyidp/cvtools.py +++ b/easyidp/cvtools.py @@ -7,7 +7,7 @@ warnings.filterwarnings("ignore", message="The array interface is deprecated and will no longer work in Shapely 2.0") -def imarray_crop(imarray, polygon_hv, outside_value=0): +def imarray_crop(imarray, polygon_hv, nodata_value=0, transparent_layer=None): """crop a given ndarray image by given polygon pixel positions Parameters @@ -38,7 +38,7 @@ def imarray_crop(imarray, polygon_hv, outside_value=0): >>> band_container = [] >>> for i in range(0, 6): >>> band = multi_spect_imarray[:,:,i] - >>> out, offset = idp.cvtools.imarray_crop(band, polygon_hv, outside_value=your_geotiff.header['nodata']) + >>> out, offset = idp.cvtools.imarray_crop(band, polygon_hv, nodata_value=your_geotiff.header['nodata']) >>> band_container.append(out) >>> final_out = np.dstack(band_container) @@ -50,10 +50,16 @@ def imarray_crop(imarray, polygon_hv, outside_value=0): it is reverted to the numpy imarray axis. horzontal = numpy axis 1, vertical = numpy axis 0. - outside_value: int | float + nodata_value: int | float | specify exact value outside the polgyon, default 0. | But for some DSM geotiff, it could be -10000.0, depends on the geotiff meta infomation + transparent_layer: None | int + | specify one layer to store transparency | outside polygon region, default None + | None: no alpha layer + | 0-x : specific alpha layer + | -1 : the last layer + returns ------- imarray_out : ndarray @@ -117,85 +123,96 @@ def imarray_crop(imarray, polygon_hv, outside_value=0): # remove (160, 160, 1) such fake 3 dimention imarray = np.squeeze(imarray) - dim = len(imarray.shape) - if dim == 2: - # only has 2 dimensions - # e.g. DSM 1 band only, other value outside polygon = empty value - - # here need to reverse - # imarray.shape -> (h, w), but poly2mask need <- (w, h) - roi_cropped = imarray[roi_top_left_offset[1]:roi_max[1], - roi_top_left_offset[0]:roi_max[0]] - rh = roi_cropped.shape[0] - rw = roi_cropped.shape[1] - mask = poly2mask((rw, rh), roi_rm_offset) - - roi_cropped[~mask] = outside_value - imarray_out = roi_cropped - - elif dim == 3: - # has 3 dimensions - # e.g. DOM with RGB or RGBA band, other value outside changed alpha layer to 0 - # coordinate xy reverted between easyidp and numpy - roi_cropped = imarray[roi_top_left_offset[1]:roi_max[1], - roi_top_left_offset[0]:roi_max[0]] - - rh = roi_cropped.shape[0] - rw = roi_cropped.shape[1] - layer_num = roi_cropped.shape[2] - - # here need to reverse - # imarray.shape -> (h, w), but poly2mask need <- (w, h) - mask = poly2mask((rw, rh), roi_rm_offset) - - if layer_num == 3: - # DOM without alpha layer - RGB - # but easyidp will add masked alpha layer to output. - - # change mask data type to fit with the image data type - if np.issubdtype(roi_cropped.dtype, np.integer) and roi_cropped.min() >= 0 and roi_cropped.max() <= 255: - # the image is 0-255 & int type - roi_cropped = roi_cropped.astype(np.uint8) - mask = mask.astype(np.uint8) * 255 - elif np.issubdtype(roi_cropped.dtype, np.floating) and roi_cropped.min() >= 0 and roi_cropped.max() <= 1: - # the image is 0-1 & float type - mask = mask.astype(roi_cropped.dtype) - else: - raise AttributeError(f"Can not handle RGB imarray ranges ({roi_cropped.min()} - {roi_cropped.max()}) with dtype='{roi_cropped.dtype}', " - f"expected (0-1) with dtype='float' or (0-255) with dtype='int'") - - # merge alpha mask with cropped images - imarray_out = np.concatenate([roi_cropped, mask[:, :, None]], axis=2) - - elif layer_num == 4: - # DOM with alpha layer - RGBA + + if dim == 2: + # only has x-y 2 axis, one layer image + layer_num = 1 + elif dim == 3: + # has x-y-z 3 axis, the most common geotiff image + layer_num = imarray.shape[2] + else: + raise ValueError( + f"Current image shape {imarray.shape} is not standard image array, " + f", please check your input image source or whether ROI is smaller than one pixel.") + + ################# + # doing croping # + ################# + # coordinate xy reverted between easyidp and numpy + # imarray.shape -> (h, w), but poly2mask need <- (w, h) + roi_cropped = imarray[roi_top_left_offset[1]:roi_max[1], + roi_top_left_offset[0]:roi_max[0]] + + rh = roi_cropped.shape[0] + rw = roi_cropped.shape[1] + mask = poly2mask((rw, rh), roi_rm_offset) + + im_dtype = roi_cropped.dtype + if np.issubdtype(im_dtype, np.integer): # 整数类型 + mask_max_value = np.iinfo(im_dtype).max + elif np.issubdtype(im_dtype, np.floating): # 浮点类型 + mask_max_value = 1.0 + else: + raise TypeError(f"Unsupported dtype: {im_dtype}") + + ####################################### + # combining image data with mask data # + ####################################### + # using extra layer to store outside + # keep original data, outside using alpha to represent + if transparent_layer is not None: + # ------------- + # Layer 0 -> transparent_layer 0 | Layer_num = 1 + # ------------- + # Layer 1 -> transparent_layer 1 | Layer_num = 2 + # ------------- + # Layer 2 -> transparent_layer 2 | Layer_num = 3 + # ------------- + + # check layers + # input has alpha layer (img.dim == transparent_layer) + # input not has alpha layer, add new layer (img.dim + 1 == transparent_layer) + # raise warning transparent layer not the last layer (img.dim > transparent layer) + if transparent_layer == -1: + transparent_layer = layer_num - 1 + + if transparent_layer == layer_num - 1: + # input already has alpha layer, need to mix two alpha layers + # for example, input is RGBA image # merge orginal mask with polygon_hv mask - original_mask = roi_cropped[:, :, 3].copy() + original_mask = roi_cropped[:, :, transparent_layer].copy() original_mask = original_mask > 0 # change type to bool merged_mask = original_mask * mask # bool = bool * bool - # change mask data type to fit with the image data type - if np.issubdtype(roi_cropped.dtype, np.integer) and roi_cropped.min() >= 0 and roi_cropped.max() <= 255: - # the image is 0-255 & int type - roi_cropped = roi_cropped.astype(np.uint8) - merged_mask = merged_mask.astype(np.uint8) * 255 - elif np.issubdtype(roi_cropped.dtype, np.floating) and roi_cropped.min() >= 0 and roi_cropped.max() <= 1: - # the image is 0-1 & float type - merged_mask = merged_mask.astype(roi_cropped.dtype) - else: - raise AttributeError(f"Can not handle RGB imarray ranges ({roi_cropped.min()} - {roi_cropped.max()}) with dtype='{roi_cropped.dtype}', " - f"expected (0-1) with dtype='float' or (0-255) with dtype='int'") - - imarray_out = np.dstack([roi_cropped[:,:, 0:3], merged_mask]) + mask_converted = merged_mask.astype(im_dtype) * mask_max_value + + imarray_out = roi_cropped.copy() + imarray_out[:,:, transparent_layer] = mask_converted + + + elif transparent_layer == layer_num: # layer_num -1 + 1 + # input not has alpha layer, add new layer + mask_converted = mask.astype(im_dtype) * mask_max_value + + imarray_out = np.dstack([roi_cropped, mask_converted]) + else: - raise TypeError(f'Unable to solve the layer/band number {layer_num}, only one band DSM or 3|4 band RGB|RGBA DOM are acceptable') + raise IndexError(f"Transparent layer (index={transparent_layer}) is not the last layer (layer number: {layer_num})") + + + # not using extra layer to store outside values, + # override each layer with nodata else: - raise ValueError( - f"Only image dimention=2 (mxn) or 3(mxnxd) are accepted, not current" - f"[shape={imarray.shape} dim={dim}], please check whether your ROI " - f"is smaller than one pixel.") + imarray_out = roi_cropped.copy() + # loop each layer, change values + # if layer_num == 1: + # imarray_out[~mask] = nodata_value + # else: + # for i in range(layer_num): + # imarray_out[~mask, i] = nodata_value + imarray_out[~mask] = nodata_value return imarray_out, roi_top_left_offset diff --git a/easyidp/data.py b/easyidp/data.py index 004b181..c618a41 100644 --- a/easyidp/data.py +++ b/easyidp/data.py @@ -574,6 +574,7 @@ def __init__(self, test_out="./tests/out"): * ``.shp.roi_prj`` * ``.shp.testutm_shp`` * ``.shp.testutm_prj`` + * ``.shp.mlayer_shp`` **pcd test module** @@ -588,17 +589,22 @@ def __init__(self, test_out="./tests/out"): * ``.pcd.maize_laz`` * ``.pcd.maize_ply`` + **roi test module** * ``.roi.dxf`` * ``.roi.lxyz_txt`` * ``.roi.xyz_txt`` + **geotiff test module** * ``.tiff.soyweed_part`` + * ``.tiff.mlayer_ndvi `` + * ``.tiff.mlayer_multi`` * ``.tiff.out`` + **metashape test module** * ``.metashape.goya_psx`` @@ -776,6 +782,9 @@ def __init__(self, data_dir, test_out): self.jp_crs_shp = data_dir / "shp_test" / "jp_crs.shp" self.jp_crs_prj = data_dir / "shp_test" / "jp_crs.prj" + # for multi-spectral testing + self.mlayer_shp = data_dir / "shp_test" / "mlayer_roi.shp" + if isinstance(test_out, str): test_out = Path(test_out) self.out = test_out / "shp_test" @@ -828,6 +837,10 @@ def __init__(self, data_dir, test_out): self.soyweed_part = data_dir / "tiff_test" / "2_12.tif" + # for multi-spectral testing + self.mlayer_ndvi = data_dir / "tiff_test" / "mlayer_yamato_ndvi.tif" + self.mlayer_multi = data_dir / "tiff_test" / "mlayer_yamato_multi.tif" + if isinstance(test_out, str): test_out = Path(test_out) self.out = test_out / "tiff_test" diff --git a/easyidp/geotiff.py b/easyidp/geotiff.py index f79ca22..300b8b6 100644 --- a/easyidp/geotiff.py +++ b/easyidp/geotiff.py @@ -14,13 +14,18 @@ class GeoTiff(object): """A GeoTiff class, consisted by header information and file path to raw file. """ - def __init__(self, tif_path=""): + def __init__(self, tif_path="", transparent_layer=None): """The method to initialize the GeoTiff class Parameters ---------- tif_path : str | pathlib.Path, optional the path to geotiff file, by default "" + Transparent_layer: None | int, optional + the transparent | alpha layer of given GeoTiff file, by default None + | None: no alpha layer + | 0-x : specific alpha layer + | -1 : the last layer Example ------- @@ -75,6 +80,9 @@ def __init__(self, tif_path=""): #: The numpy ndarray of GeoTiff images self.imarray = None + #: The layer to represent transparency / alpha + self.transparent_layer = None + if tif_path != "": self.read_geotiff(tif_path) @@ -433,6 +441,7 @@ def crop_polygon(self, polygon_hv, is_geo=True, save_path=None): """ self._not_empty() + if is_geo: poly_px = geo2pixel(polygon_hv, self.header, return_index=True) else: @@ -468,12 +477,15 @@ def crop_polygon(self, polygon_hv, is_geo=True, save_path=None): poly_offseted_px = poly_px - bbox_left_top imarray_out, _ = idp.cvtools.imarray_crop( imarray_bbox, poly_offseted_px, - outside_value=self.header['nodata'] + nodata_value=self.header['nodata'], + transparent_layer=self.transparent_layer ) # check if need save geotiff - if save_path is not None and os.path.splitext(save_path)[-1] == ".tif": - save_geotiff(self.header, imarray_out, bbox_left_top, save_path) + if save_path is not None: + save_path = Path(save_path) + if save_path.suffix == ".tif": + save_geotiff(self.header, imarray_out, bbox_left_top, save_path) return imarray_out @@ -595,8 +607,10 @@ def crop_rectangle(self, left, top, w, h, is_geo=True, save_path=None): out = tifffile_crop(page, top, left, h, w) # check if need save geotiff - if save_path is not None and os.path.splitext(save_path)[-1] == ".tif": - save_geotiff(self.header, out, np.array([left, top]), save_path) + if save_path is not None: + save_path = Path(save_path) + if save_path.suffix == ".tif": + save_geotiff(self.header, out, np.array([left, top]), save_path) return out @@ -891,7 +905,7 @@ def get_header(tif_path): with tf.TiffFile(tif_path) as tif: header = {} # keys: 'width', 'height', 'dim', 'scale', 'tie_point', - # 'nodata', 'crs', 'dtype', 'band_num', + # 'nodata', 'crs', 'dtype', 'band_num' # for export: # 'tags', 'photometric', 'planarconfig', 'compression' page = tif.pages[0] @@ -1637,8 +1651,18 @@ def save_geotiff(header, imarray, left_top_corner, save_path): """ extratags = _offset_geotiff_extratags(header, left_top_corner) - file_ext = os.path.splitext(save_path)[-1] + if not isinstance(save_path, Path): + save_path = Path(save_path) + + # create folder if parent folder not exists + file_ext = save_path.suffix + parent_folder = save_path.parent + if file_ext == ".tif": + # create folder + if not parent_folder.exists(): + parent_folder.mkdir() + # write geotiff with tf.TiffWriter(save_path) as wtif: wtif.write(data=imarray, diff --git a/tests/test_cvtools.py b/tests/test_cvtools.py index c6521bb..7c71cf3 100644 --- a/tests/test_cvtools.py +++ b/tests/test_cvtools.py @@ -136,7 +136,7 @@ def test_imarray_crop_2d_rgb_rgba(): # 3d rgb with int type imarray_rgb_int = plt.imread(photo_path) # imarray_rgb.shape == (3456, 4608, 3) - img_out_rgb_int, offsets_rgb = idp.cvtools.imarray_crop(imarray_rgb_int, roi) + img_out_rgb_int, offsets_rgb = idp.cvtools.imarray_crop(imarray_rgb_int, roi, transparent_layer=3) ax[0, 1].imshow(img_out_rgb_int) ax[0, 1].set_title('rgb(int)') @@ -146,7 +146,7 @@ def test_imarray_crop_2d_rgb_rgba(): # 3d rgb with float type imarray_rgb_float = img_as_float(imarray_rgb_int) - img_out_rgb_float, offsets_rgbf = idp.cvtools.imarray_crop(imarray_rgb_float, roi) + img_out_rgb_float, offsets_rgbf = idp.cvtools.imarray_crop(imarray_rgb_float, roi, transparent_layer=3) ax[1, 1].imshow(img_out_rgb_float) ax[1, 1].set_title('rgb(float)') @@ -179,7 +179,7 @@ def test_imarray_crop_2d_rgb_rgba(): imarray_rgba_int = np.dstack((imarray_rgb_int, np.ones((3456, 4608)) * 255)).astype(np.uint8) # imarray_rgba.shape == (3456, 4608, 4) - im_out_rgba_int, offsets_rgba = idp.cvtools.imarray_crop(imarray_rgba_int, roi) + im_out_rgba_int, offsets_rgba = idp.cvtools.imarray_crop(imarray_rgba_int, roi, transparent_layer=3) ax[0, 2].imshow(im_out_rgba_int) ax[0, 2].set_title('rgba(int)') @@ -188,7 +188,7 @@ def test_imarray_crop_2d_rgb_rgba(): imarray_rgba_float = img_as_float(imarray_rgba_int) - im_out_rgba_float, offsets_rgba = idp.cvtools.imarray_crop(imarray_rgba_float, roi) + im_out_rgba_float, offsets_rgba = idp.cvtools.imarray_crop(imarray_rgba_float, roi, transparent_layer=3) ax[1, 2].imshow(im_out_rgba_float) ax[1, 2].set_title('rgba(float)') @@ -213,11 +213,6 @@ def test_imarray_crop_2d_rgb_rgba_error(): imarray_init= np.random.randint(0,255,(10,10,4), dtype=np.uint8) str_imarray = imarray_init.astype(str) - int_over_255 = imarray_init.astype(np.uint16) + 255 - int_below_0 = imarray_init.astype(np.int16) - 128 - float_over_1 = imarray_init.astype(np.float32) / 255 + 1 - float_below_0 = imarray_init.astype(np.float32) / 255 - 0.5 - img_coord = np.asarray([[0,0],[1,1],[4,4]]) # str type error @@ -229,82 +224,11 @@ def test_imarray_crop_2d_rgb_rgba_error(): ): imarray, offsets = idp.cvtools.imarray_crop(str_imarray, img_coord) - - # int >255 error - with pytest.raises( - AttributeError, - match=re.escape( - "Can not handle RGB imarray ranges (260 - 497) with dtype='uint16', " - ) - ): - # dim = 4 - imarray, offsets = idp.cvtools.imarray_crop(int_over_255, img_coord) - - with pytest.raises( - AttributeError, - match=re.escape( - "Can not handle RGB imarray ranges (260 - 497) with dtype='uint16', " - ) - ): - # dim = 3 - imarray, offsets = idp.cvtools.imarray_crop(int_over_255[:,:,0:3], img_coord) - - # int < 0 error - with pytest.raises( - AttributeError, - match=re.escape( - "Can not handle RGB imarray ranges (-123 - 114) with dtype='int16', " - ) - ): - # dim = 4 - imarray, offsets = idp.cvtools.imarray_crop(int_below_0, img_coord) - - with pytest.raises( - AttributeError, - match=re.escape( - "Can not handle RGB imarray ranges (-123 - 114) with dtype='int16', " - ) - ): - # dim = 3 - imarray, offsets = idp.cvtools.imarray_crop(int_below_0[:,:,0:3], img_coord) - - # float > 1 error - with pytest.raises( - AttributeError, - match=re.escape( - "Can not handle RGB imarray ranges (1.0196079015731812 - 1.9490196704864502) with dtype='float32', " - ) - ): - # dim = 4 - imarray, offsets = idp.cvtools.imarray_crop(float_over_1, img_coord) - - with pytest.raises( - AttributeError, - match=re.escape( - "Can not handle RGB imarray ranges (1.0196079015731812 - 1.9490196704864502) with dtype='float32', " - ) - ): - # dim = 3 - imarray, offsets = idp.cvtools.imarray_crop(float_over_1[:,:,0:3], img_coord) - - # float < 0 error - with pytest.raises( - AttributeError, - match=re.escape( - "Can not handle RGB imarray ranges (-0.4803921580314636 - 0.4490196108818054) with dtype='float32', " - ) - ): - # dim = 4 - imarray, offsets = idp.cvtools.imarray_crop(float_below_0, img_coord) - - with pytest.raises( - AttributeError, - match=re.escape( - "Can not handle RGB imarray ranges (-0.4803921580314636 - 0.4490196108818054) with dtype='float32', " - ) - ): - # dim = 3 - imarray, offsets = idp.cvtools.imarray_crop(float_below_0[:,:,0:3], img_coord) + # These test not compatative to multi-spectal GeoTiffs, deprecated + # int_over_255 = imarray_init.astype(np.uint16) + 255 + # int_below_0 = imarray_init.astype(np.int16) - 128 + # float_over_1 = imarray_init.astype(np.float32) / 255 + 1 + # float_below_0 = imarray_init.astype(np.float32) / 255 - 0.5 def test_imarray_crop_handle_polygon_hv_with_float(): @@ -319,7 +243,8 @@ def test_imarray_crop_handle_polygon_hv_with_float(): np.random.seed(0) img_np = np.random.randint(low=0, high=255, size=(10000,10000,3), dtype=np.uint8) - imarray, offsets = idp.cvtools.imarray_crop(img_np, img_coord) + # the fourth layer (id=3) as new alpha layer + imarray, offsets = idp.cvtools.imarray_crop(img_np, img_coord, transparent_layer=3) assert imarray.shape == (1789, 1811, 4) assert imarray.dtype == np.uint8 @@ -381,7 +306,7 @@ def test_roi_smaller_than_one_pixel_error(): # disucssion #39 with pytest.raises( ValueError, match=re.escape( - "Only image dimention=2 (mxn) or 3(mxnxd) are accepted, not current" + f"Current image shape {one_dim_imarray.shape} is not standard image array" ) ): idp.cvtools.imarray_crop(one_dim_imarray, polygon_hv) \ No newline at end of file diff --git a/tests/test_geotiff.py b/tests/test_geotiff.py index c7efe41..7d5a77f 100644 --- a/tests/test_geotiff.py +++ b/tests/test_geotiff.py @@ -43,7 +43,6 @@ def test_def_get_header(): assert lotus_part["tie_point"][0] == 368024.0839 assert lotus_part["tie_point"][1] == 3955479.7512 - def test_def_get_imarray(): maize_part_np = idp.geotiff.get_imarray(test_data.pix4d.maize_dom) assert maize_part_np.shape == (722, 836, 4) @@ -621,6 +620,40 @@ def test_class_crop_rois(): assert (tif_out_folder / "N1W1.tif").exists() assert out_dict["N2E2"].shape == (320, 320, 4) +def test_class_crop_rois_multispec(): + roi = idp.ROI(test_data.shp.mlayer_shp) + + # 5 layers multispectral with 5th as alpha + multi_tiff = idp.GeoTiff(test_data.tiff.mlayer_multi) + multi_tiff.transparent_layer = -1 + + tif_out_folder = test_data.tiff.out / "multi_crop" + if tif_out_folder.exists(): + shutil.rmtree(tif_out_folder) + tif_out_folder.mkdir() + + out_dict = multi_tiff.crop_rois(roi, save_folder=tif_out_folder) + assert len(out_dict) == 12 + assert (tif_out_folder / "2.tif").exists() + assert out_dict["2"].shape == (180, 180, 5) + + +def test_class_crop_rois_ndvi_special(): + roi = idp.ROI(test_data.shp.mlayer_shp) + # processing multispectral geotiff -> 2 layer ndvi (second as alpha) + ndvi_tiff = idp.GeoTiff(test_data.tiff.mlayer_ndvi) + + tif_out_folder = test_data.tiff.out / "ndvi_crop" + if tif_out_folder.exists(): + shutil.rmtree(tif_out_folder) + tif_out_folder.mkdir() + + out_dict = ndvi_tiff.crop_rois(roi, save_folder=tif_out_folder) + assert len(out_dict) == 12 + assert (tif_out_folder / "2.tif").exists() + assert out_dict["2"].shape == (180, 180, 2) + + def test_class_geo2pixel2geo_executable(): roi = idp.ROI(test_data.shp.lotus_shp, name_field=0) @@ -633,4 +666,4 @@ def test_class_geo2pixel2geo_executable(): roi_test_back = dom.pixel2geo(roi_test_pixel) - np.testing.assert_almost_equal(roi_test, roi_test_back, decimal=5) + np.testing.assert_almost_equal(roi_test, roi_test_back, decimal=5) \ No newline at end of file