Skip to content

Commit

Permalink
rough support multi-spectral geotiff
Browse files Browse the repository at this point in the history 
Partly support #115
  • Loading branch information
@HowcanoeWang
HowcanoeWang committed Oct 21, 2024
1 parent cd0f3d8 commit f6798da
Show file tree
Hide file tree
Showing 5 changed files with 183 additions and 171 deletions.
165 changes: 91 additions & 74 deletions easyidp/cvtools.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down
13 changes: 13 additions & 0 deletions easyidp/data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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**
Expand All @@ -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``
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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"
Expand Down
40 changes: 32 additions & 8 deletions easyidp/geotiff.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
-------
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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,
Expand Down
Loading

0 comments on commit f6798da

Please sign in to comment.