masking in rio calc


Gregory, Matthew
 

Hi all,

I'm wondering if I'm using rio-calc incorrectly to achieve what I'd like to do. I have a floating-point raster that ranges from 0.0-1.0 with areas in the raster that are masked. I'd like to simply create a 0|1 threshold raster based on 0.5 as a threshold value *and* retain the mask which I'd like to set to 255 in a uint8 output. My command is:

rio calc
-f GTiff
-t uint8
--overwrite
--masked
--profile nodata=255
"(>= (read 1) 0.5)"
raw.tif
threshold.tif

When I run this, I believe the ">=" in the expression gives me a np.bool type which then gets filled here (https://github.com/mapbox/rasterio/blob/master/rasterio/rio/calc.py#L189) before being cast to the output dtype. When I swap these two code chunks: type conversion (lines 194-198) before the masking (lines 189-192), I get the expected result. But I'm guessing I'm introducing some unintended consequences.

One other minor point: In order to use my specified profile nodata value, I had to first cast it to an int (line 192):

res = res.filled(int(kwargs['nodata'])))

Thanks for any advice,
matt


Sean Gillies
 

Hi Matt,

On Tue, Nov 17, 2020 at 5:57 PM Gregory, Matthew <matt.gregory@...> wrote:
Hi all,

I'm wondering if I'm using rio-calc incorrectly to achieve what I'd like to do.  I have a floating-point raster that ranges from 0.0-1.0 with areas in the raster that are masked.  I'd like to simply create a 0|1 threshold raster based on 0.5 as a threshold value *and* retain the mask which I'd like to set to 255 in a uint8 output.  My command is:

  rio calc
    -f GTiff
    -t uint8
    --overwrite
    --masked
    --profile nodata=255
    "(>= (read 1) 0.5)"
    raw.tif
    threshold.tif

When I run this, I believe the ">=" in the expression gives me a np.bool type which then gets filled here (https://github.com/mapbox/rasterio/blob/master/rasterio/rio/calc.py#L189) before being cast to the output dtype.  When I swap these two code chunks: type conversion (lines 194-198) before the masking (lines 189-192), I get the expected result.  But I'm guessing I'm introducing some unintended consequences.

One other minor point: In order to use my specified profile nodata value, I had to first cast it to an int (line 192):

  res = res.filled(int(kwargs['nodata'])))

Thanks for any advice,
matt

You may get better results using numpy.where, like "(where (>= (read 1) 0.5) 1 0)".

I'll dig into the filling code, that looks like a bug to me.

--
Sean Gillies


Sean Gillies
 

Update: https://github.com/mapbox/rasterio/issues/2041. Thanks for bringing this up, Matt.

On Wed, Nov 18, 2020 at 1:54 PM Sean Gillies <sean.gillies@...> wrote:
Hi Matt,

On Tue, Nov 17, 2020 at 5:57 PM Gregory, Matthew <matt.gregory@...> wrote:
Hi all,

I'm wondering if I'm using rio-calc incorrectly to achieve what I'd like to do.  I have a floating-point raster that ranges from 0.0-1.0 with areas in the raster that are masked.  I'd like to simply create a 0|1 threshold raster based on 0.5 as a threshold value *and* retain the mask which I'd like to set to 255 in a uint8 output.  My command is:

  rio calc
    -f GTiff
    -t uint8
    --overwrite
    --masked
    --profile nodata=255
    "(>= (read 1) 0.5)"
    raw.tif
    threshold.tif

When I run this, I believe the ">=" in the expression gives me a np.bool type which then gets filled here (https://github.com/mapbox/rasterio/blob/master/rasterio/rio/calc.py#L189) before being cast to the output dtype.  When I swap these two code chunks: type conversion (lines 194-198) before the masking (lines 189-192), I get the expected result.  But I'm guessing I'm introducing some unintended consequences.

One other minor point: In order to use my specified profile nodata value, I had to first cast it to an int (line 192):

  res = res.filled(int(kwargs['nodata'])))

Thanks for any advice,
matt

You may get better results using numpy.where, like "(where (>= (read 1) 0.5) 1 0)".

I'll dig into the filling code, that looks like a bug to me.

--
Sean Gillies


Gregory, Matthew
 

Great, thanks Sean.

 

As I’m sure you figured out from my original example, but I wasn’t explicit about, was that filling with 255 when the numpy dtype was bool resulted in True, which then gets set to 1 on the cast to uint8.  Therefore, all input mask pixels become 1 in the output.

 

matt

 

-----------------------

 

Update: https://github.com/mapbox/rasterio/issues/2041. Thanks for bringing this up, Matt.

 

On Wed, Nov 18, 2020 at 1:54 PM Sean Gillies <sean.gillies@...> wrote:

Hi Matt,

 

On Tue, Nov 17, 2020 at 5:57 PM Gregory, Matthew <matt.gregory@...> wrote:

Hi all,

I'm wondering if I'm using rio-calc incorrectly to achieve what I'd like to do.  I have a floating-point raster that ranges from 0.0-1.0 with areas in the raster that are masked.  I'd like to simply create a 0|1 threshold raster based on 0.5 as a threshold value *and* retain the mask which I'd like to set to 255 in a uint8 output.  My command is:

  rio calc
    -f GTiff
    -t uint8
    --overwrite
    --masked
    --profile nodata=255
    "(>= (read 1) 0.5)"
    raw.tif
    threshold.tif

When I run this, I believe the ">=" in the expression gives me a np.bool type which then gets filled here (https://github.com/mapbox/rasterio/blob/master/rasterio/rio/calc.py#L189) before being cast to the output dtype.  When I swap these two code chunks: type conversion (lines 194-198) before the masking (lines 189-192), I get the expected result.  But I'm guessing I'm introducing some unintended consequences.

One other minor point: In order to use my specified profile nodata value, I had to first cast it to an int (line 192):

  res = res.filled(int(kwargs['nodata'])))

Thanks for any advice,
matt

 

You may get better results using numpy.where, like "(where (>= (read 1) 0.5) 1 0)".

 

I'll dig into the filling code, that looks like a bug to me.

 

--

Sean Gillies