GFZ / arosics

AROSICS - Automated and Robust Open-Source Image Co-Registration Software
https://git.gfz-potsdam.de/danschef/arosics
Apache License 2.0
137 stars 26 forks source link

'TypeError: ufunc 'invert' not supported for the input types' when using input mask #33

Closed himal-ella closed 7 months ago

himal-ella commented 7 months ago

Description Hi, I have been using arosics to look at displacements between images but have been running into some problems when using a mask layer. The same input images always works without a mask and works from some mask layers, however seemingly some mask shapes throw the TypeError below:

_TypeError: ufunc 'invert' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''_

I have checked that the input masks always have the same dimensions and coordinate systems as the input images and are binary 0,1 tif images; the error is also not related to whether the mask area touches the edge of the image. The masks are generated from polygons in QGIS and for some masks the code runs fine and for others it throws the error.

I am not sure what is going on and can't work out the pattern as to why certain masks are failing. Many thanks, Ella

What I Did arosics code:

kwargs = {
    'grid_res'     : 10,                          
    'window_size'  : (30,30),                  
    'nodata'       : (0,0),                                                                                         
    'mask_baddata_tgt' : mask,                     
}

FT = COREG_LOCAL(image_1,image_2,**kwargs)
table = FT.CoRegPoints_table

Mask Generation

binary_mask = features.rasterize(
        [(geom, 0) for geom in polygon.geometry],
        out_shape=out_shape,
        transform=transform,
        fill=1,
        all_touched=True,
        dtype=np.uint16,
    )

Error message:

TypeError                                 Traceback (most recent call last)
(file:///C:/Code/)...FT_script.ipynb Cell 9 line 1
      [2](vscode-notebook-cell:/...line=1) kwargs = {
      [3](vscode-notebook-cell:/c...line=2)     'grid_res'     : 10,                         
      [4](vscode-notebook-cell:/c...line=3)     'window_size'  : (30,30),                       
   (...)
     [10](vscode-notebook-cell:/...line=9)     'mask_baddata_tgt' : mask,          
     [11](vscode-notebook-cell:/c...line=10) }
     [13](vscode-notebook-cell:/c...line=12) FT = COREG_LOCAL(image_1,image_2,**kwargs)
---> [14](vscode-notebook-cell:/c...line=13) table = FT.CoRegPoints_table

F(file:///C:/.....Lib/site-packages/arosics/CoReg_local.py:466), in COREG_LOCAL.CoRegPoints_table(self)
    [459](file:///C:/....arosics/CoReg_local.py:459) @property
    [460](file:///C:/....arosics/CoReg_local.py:460) def CoRegPoints_table(self) -> GeoDataFrame:
    [461](file:///C/....arosics/CoReg_local.py:461)     
"""Return a GeoDataFrame containing all the results from coregistration for all points in the tie point grid.
    [462](file:///C:/....arosics/CoReg_local.py:462) 
    [463](file:///C:/....arosics/CoReg_local.py:463)     Columns of the GeoDataFrame: 
'geometry','POINT_ID','X_IM','Y_IM','X_MAP','Y_MAP','X_WIN_SIZE', 'Y_WIN_SIZE',
    [464](file:///C:/...arosics/CoReg_local.py:464)                                 
 'X_SHIFT_PX','Y_SHIFT_PX', 'X_SHIFT_M', 'Y_SHIFT_M', 'ABS_SHIFT' and 'ANGLE'
    [465](file:///C:/...arosics/CoReg_local.py:465)     """

--> [466](file:///C:.../arosics/CoReg_local.py:466)     return self.tiepoint_grid.CoRegPoints_table

(file:///C:..../arosics/CoReg_local.py:456), in COREG_LOCAL.tiepoint_grid(self)
    [454](file:///C:/...arosics/CoReg_local.py:454)     return self._tiepoint_grid
    [455](file:///C:/...arosics/CoReg_local.py:455) else:
...
--> [366](file:///C:/...Lib/site-packages/pandas/core/internals/blocks.py:366)     result = func(self.values, **kwargs)
    [368](file:///C:/...site-packages/pandas/core/internals/blocks.py:368)     result = maybe_coerce_values(result)
    [369](file:///C:/...site-packages/pandas/core/internals/blocks.py:369)     return self._split_op_result(result)

TypeError: ufunc 'invert' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
danschef commented 7 months ago

What is the data type of the mask you provide? Your call of features.rasterize looks like it produces a mask in unsigned integer 16 bit. I suggest to pass boolean masks to AROSICS, i.e., the data type should be bool.

If that fixes your issue, I will adapt AROSICS to return a more useful error message if an unexpected data type if given.

himal-ella commented 7 months ago

Thanks, yes my data is Int8 with 0 1 values as the rasterio features.rasterize function doesn't support saving as a bool data type. I tried loading the mask and then converting it to a boolean, as below, however when I use this as an input to arosics it still throws the same error as before.

with rasterio.open(mask_path) as src: mask_raster = src.read(1)

crs = src.crs
transform = src.transform
width = src.width
height = src.height

mask_input = maskraster.astype(np.bool) #Convert raster to bool

danschef commented 7 months ago

Sorry, I cannot reproduce it.

from geoarray import GeoArray
from arosics import COREG_LOCAL
import numpy as np

# some sample data
ref = GeoArray('/path/to/reference_image.tif')
ref = GeoArray('/path/to/target_image.tif')

# create an empty bad data mask with the same dimension like the target image
mask = np.zeros((tgt.rows, tgt.columns), dtype=bool)

# fill in some bad data
mask[1000:2000, 1000:2000] = True

# make the mask a GeoArray instance (you can also pass a the file path of you bad data mask to COREG_LOCAL)
gA_mask = GeoArray(
    mask,
    geotransform=tgt.gt,
    projection=tgt.prj
)

# call COREG_LOCAL
CRL = COREG_LOCAL(ref, tgt, gri_res=100, mask_baddata_tgt=gA_mask)
CRL.calculate_spatial_shifts()
CRL.view_CoRegPoints()

The mask looks like this: grafik

And as expected, no shift vectors are computed at the bad-data-area, so everything works as it should: grafik

Could you provide a simple code snippet that allows me to reproduce the issue? Maybe I can also reproduce it if you would upload your mask file.

himal-ella commented 7 months ago

Hi, Thanks for your help so far, apologies for the delayed reply. I have attached some example data which should allow you to reproduce the issue.

Here is the code for the function I used to create the mask and I have also uploaded the mask itself here as well as the original input polygon.

Many thanks, Ella

mask_problem_data.zip

## Mask code ##
def create_mask_layer(im_target, shapefile_path, MASK_export_path):
    with rasterio.open(im_target) as src:
        raster = src.read(1)
        transform = src.transform
        crs = src.crs

    polygon = gpd.read_file(shapefile_path)

    out_shape = raster.shape
    binary_mask = features.rasterize(
        [(geom, 0) for geom in polygon.geometry],
        out_shape=out_shape,
        transform=transform,
        fill=1,
        all_touched=True,
        dtype=np.int8,
    )

    profile = {
        'driver': 'GTiff',
        'width': out_shape[1],
        'height': out_shape[0],
        'count': 1,
        'dtype': np.int8,
        'crs': crs,
        'transform': transform
    }

    with rasterio.open(MASK_export_path, 'w', **profile) as dst:
        dst.write(binary_mask, 1)

create_mask_layer(im_target, shapefile_path, MASK_export_path)

## Arosics code ##

kwargs = {
    'grid_res'     : 10,    
    'window_size'  : (30,30),    
    'nodata'       : (0,0), 
    'min_reliability': 0,    
    'max_shift': 10,         
    'max_iter' : 10,       
    'tieP_filter_level' : 0,         
    'mask_baddata_tgt' : mask,                      
}

FT = COREG_LOCAL(image_1,image_2,**kwargs)
table = FT.CoRegPoints_table
danschef commented 7 months ago

Thanks, now I can reproduce it. The actual traceback is:

../arosics/CoReg_local.py:490: in calculate_spatial_shifts
    self._tiepoint_grid.get_CoRegPoints_table()
../arosics/Tie_Point_Grid.py:352: in get_CoRegPoints_table
    GDF = self._exclude_bad_XYpos(GDF)
../arosics/Tie_Point_Grid.py:275: in _exclude_bad_XYpos
    GDF = GDF[(~GDF['REF_BADDATA']) & (~GDF['TGT_BADDATA'])]
../../../../mambaforge/envs/arosics/lib/python3.11/site-packages/pandas/core/generic.py:1513: in __invert__
    new_data = self._mgr.apply(operator.invert)
../../../../mambaforge/envs/arosics/lib/python3.11/site-packages/pandas/core/internals/managers.py:352: in apply
    applied = b.apply(f, **kwargs)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

self = NumpyBlock: 1 dtype: float64, func = <built-in function invert>
kwargs = {}

    @final
    def apply(self, func, **kwargs) -> list[Block]:
        """
        apply the function to my values; return a block if we are not
        one
        """
>       result = func(self.values, **kwargs)
E       TypeError: ufunc 'invert' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

../../../../mambaforge/envs/arosics/lib/python3.11/site-packages/pandas/core/internals/blocks.py:366: TypeError

I am investigating...

danschef commented 7 months ago

Merge request in progress: https://git.gfz-potsdam.de/danschef/arosics/-/merge_requests/42

danschef commented 7 months ago

Ok, thanks for reporting, the issue is fixed in AROSICS v1.9.4 which will become available very soon.

himal-ella commented 7 months ago

Okay, thanks so much for the help!