mdbartos / pysheds

:earth_americas: Simple and fast watershed delineation in python.
GNU General Public License v3.0
708 stars 191 forks source link

Grid.from_raster dtype error #195

Open T4mmi opened 2 years ago

T4mmi commented 2 years ago

Hi, I got a weird error trying to read a DEM raster:

In [65]: grid = Grid.from_raster('MNT_25.tif')

File ~\AppData\Roaming\Python\Python39\site-packages\pysheds\sview.py:89, in Raster.__new__(cls, input_array, viewfinder, metadata)
     87     assert np.min_scalar_type(viewfinder.nodata) <= obj.dtype
     88 except:
---> 89     raise TypeError('`nodata` value not representable in dtype of array.')       
     90 # Don't allow original viewfinder and metadata to be modified
     91 viewfinder = viewfinder.copy()

TypeError: `nodata` value not representable in dtype of array.

From what I understand there is something wrong with type checking ... but I cannot say if its a numpy bug or a pysheds "overdoing" since I can reproduce manually the error (note that reading using rasterio works fine):

In [64]: raster = rasterio.open('MNT_25.tif')

In [65]: obj = raster.read()
Out[65]: 
array([[[152, 151, 148, ..., 198, 198, 198],
        [156, 154, 150, ..., 198, 198, 198],
        [160, 157, 152, ..., 198, 198, 198],
        ...,
        [377, 366, 354, ..., 406, 403, 402],
        [382, 370, 359, ..., 402, 400, 398],
        [386, 375, 365, ..., 398, 396, 394]]], dtype=int16)

In [66]: raster.nodata
Out[66]: 32767

In [67]: raster.nodata.dtype
Out[67]: dtype('int16')

In [68]: obj.dtype
Out[68]: dtype('int16')

In [69]: np.min_scalar_type(raster.nodata)
Out[69]: dtype('uint16')

In [70]: np.min_scalar_type(raster.nodata) <= obj.dtype
Out[70]: False

In [71]: np.dtype('uint16') <= np.dtype('int16')
Out[71]: False

https://github.com/mdbartos/pysheds/blob/9d960962bfe742c534da6577c3353169e6b6a501/pysheds/sview.py#L293

meteoDaniel commented 2 years ago

I am having the same issue, any solutions so far ?

T4mmi commented 2 years ago

No solution so far, just a workaround: set an arbitrary value for nodata:

try:
    grid = Grid.from_raster(input_file)
except TypeError:
    grid = Grid.from_raster(input_file, nodata = -9999)
dem = grid.read_raster(input_file, band=band or 1, nodata=grid.nodata)