Date   

Rasterio 1.2b1

Sean Gillies
 

Hi all,

Wheels for Python versions 3.6-3.8 and manylinux1 and macos (built with xcode 9.3) are on PyPI this morning. I'd appreciate it if you could give it a try: python -m pip install rasterio==1.2b1.

At last, these wheels include GDAL 3.2.0 and PROJ 7.2.0. Please note that these wheels include no PROJ data. You must provide your own data or set PROJ_NETWORK=ON in your environment to use the PROJ CDN. See https://proj.org/resource_files.html for details.

I'm particularly interested in feedback about including leaving out the PROJ data. If it's going to impact the way you deploy rasterio, this is the time to let us know.

Yours,

--
Sean Gillies


resampling raster after clipping/masking

Eyal Saiet
 

Hello,
I have found references in rasterio.readthedocs.io both on how to crop/clip a raster and how to resample. I generally understand the steps and sought to combine the two steps into one:

#Clipping Raster
with rasterio.open(dem_raster_pd['rasters'][0]) as raster:
    #create a mask
    out_image,out_transform=rio.mask.mask(raster,shapes,crop=True)
    #adding meta of source
    out_meta=raster.meta


    #updading meta
    out_meta.update({"driver":"GTiff","height":out_image.shape[1],
                                     "width":out_image.shape[2],
                                    "transform":out_transform})

#Because I am resampling next I skipped writing crp_raster
   ##Rescaling to 1 m
    ras_xy=1 #(m)

   # figuring out the factor
    rs_factor=((ras_xy/out_image.shape[1],ras_xy/out_image.shape[2])
    rs_in_factor=(1/rs_factor[0],1/rs_factor[1])

    #But because I don't have a crp_raster raster object I cannot use read and apply yhr out_shape
     data_resample=crp_raster.read(out_shape=(crp_raster.count,int(crp_raster.height*rs_in_factor[0]),int(crp_raster.width*rs_in_factor[1])),resampling=Resampling.bilinear,masked=True)#,masked=True


    #Scaling image trnasform
    transform_resmp=out_transform * out_transform.scale(
       (out_image.shape[2]/data_resample.shape[-1]),
        (out_image.shape[1]/data_resample.shape[-2])
        )

    #updating meta
    out_meta_res=crp_raster.meta
    out_meta_res.update({"height":data_resample.shape[-2],"width":data_resample.shape[-1],"transform":transform_resmp})


    #New filename
    new_ras_nm=products_dir+os.sep+raster.name.split('.')[0]+"_crp_rsmp_1m.tif"

    #Saving
with rasterio.open(new_ras_nm,"w",**out_meta_res) as dest_crp_resamp:
    print('Resampled')
    dest_crp_resamp.write(data_resample)     

I am sure I am stuck over something obvious that currently I cannot see. If I could create a raster object,  I tried crp_raster=rasterio.open("w"..), then I could use the raster_objec.read(out_shape=).




Re: Rasterio 1.2a1

Yann-Sebastien Tremblay-Johnston
 

+1 LGTM. Light testing of new coordinate transformation and warping with rpcs in WSL 2 Ubuntu 20.04.

Sebastien


Rasterio 1.2a1

Sean Gillies
 

Hi all,

Some wheels and an sdist for rasterio 1.2a1 are on PyPI today. I'd appreciate some testing by adventurous folks.

The change log is behind the curve, but the new features are listed in https://github.com/mapbox/rasterio/pulls?q=is%3Apr+milestone%3A1.2.0+is%3Aclosed. There are still a few known bugs to be fixed and features to be finished before 1.2.0. For example, the merge tool will write to a file, but can't yet write datasets that are larger than your computer's memory can store.

Looking forward to feedback,

--
Sean Gillies


Re: masking in rio calc

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


Re: masking in rio calc

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


Re: masking in rio calc

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


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


Re: Sampling from one single point location using sample()

Alexander Jüstel
 

Hey Sean, 

I made it work using your tests! Thanks for that ;)

Cheers Alex

On Mon, Nov 16, 2020 at 4:03 PM Sean Gillies via groups.io <sean=mapbox.com@groups.io> wrote:
Hi Alex,

On Sun, Nov 15, 2020 at 9:06 AM Alexander Jüstel <alexander.juestel@...> wrote:
Hello, 
 
I would like to sample the value of my raster at one single point location. I have passed the coordinates as an iterable to the sample() function. However, I am not able to extract the value from the generator. The sampling works fine for more than one point. However, I need it to work for one point only. I have tried using list(gen) or next(gen) to get the value but nothing works so far.
 
Thanks for your support
Cheers
Alex

Is your approach any different than the examples in https://github.com/mapbox/rasterio/blob/master/tests/test_sampling.py?

--
Sean Gillies


Re: Sampling from one single point location using sample()

Sean Gillies
 

Hi Alex,

On Sun, Nov 15, 2020 at 9:06 AM Alexander Jüstel <alexander.juestel@...> wrote:
Hello, 
 
I would like to sample the value of my raster at one single point location. I have passed the coordinates as an iterable to the sample() function. However, I am not able to extract the value from the generator. The sampling works fine for more than one point. However, I need it to work for one point only. I have tried using list(gen) or next(gen) to get the value but nothing works so far.
 
Thanks for your support
Cheers
Alex

Is your approach any different than the examples in https://github.com/mapbox/rasterio/blob/master/tests/test_sampling.py?

--
Sean Gillies


Re: Is it possible to ignore existing overview when performing decimated read?

Sean Gillies
 

Hi Loïc,

A test that requires 3.3 would be welcome. 

The next time we have a rasterio release I will patch GDAL the GDAL library in the wheels we publish on PyPI.


On Mon, Nov 16, 2020, 3:26 AM Loïc Dutrieux <loic.dutrieux@...> wrote:
Many thanks Sean and Even,

So with rasterio compiled against gdal 3.3 I'll be able to force ignore overviews?
If that deserves a new rasterio unit test I'm happy to send a pull request in the coming days.

Cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Sean Gillies via groups.io <sean=mapbox.com@groups.io>
Sent: 13 November 2020 20:57
To: main@rasterio.groups.io
Subject: Re: [rasterio] Is it possible to ignore existing overview when performing decimated read?

Hi Even,

On Fri, Nov 13, 2020 at 10:38 AM Even Rouault <even.rouault@...<mailto:even.rouault@...>> wrote:
On vendredi 13 novembre 2020 09:05:41 CET Sean Gillies via groups.io<https://urldefense.com/v3/__http://groups.io__;!!DOxrgLBm!U8N393pmVlI170WWlMCwvmrIV7A8JSl9WHWOD7E3re7ET8O3tNpqP7DhM2SV30dTeRRzVSk$> wrote:
> HI Loïc,
>
> I don't think we currently have any option other than to remove the
> overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but
> the implementation isn't there. For example, if you used gdalwarp with
> `-ovr NONE` it will set an internal overview level to the value of -1. But
> that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass
> `OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.

Was left as an exercice for a contributor:
https://github.com/OSGeo/gdal/pull/3181<https://urldefense.com/v3/__https://github.com/OSGeo/gdal/pull/3181__;!!DOxrgLBm!U8N393pmVlI170WWlMCwvmrIV7A8JSl9WHWOD7E3re7ET8O3tNpqP7DhM2SV30dTK8LUNCg$>

Even

Thanks


Re: Is it possible to ignore existing overview when performing decimated read?

Loïc Dutrieux
 

Many thanks Sean and Even,

So with rasterio compiled against gdal 3.3 I'll be able to force ignore overviews?
If that deserves a new rasterio unit test I'm happy to send a pull request in the coming days.

Cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Sean Gillies via groups.io <sean=mapbox.com@groups.io>
Sent: 13 November 2020 20:57
To: main@rasterio.groups.io
Subject: Re: [rasterio] Is it possible to ignore existing overview when performing decimated read?

Hi Even,

On Fri, Nov 13, 2020 at 10:38 AM Even Rouault <even.rouault@spatialys.com<mailto:even.rouault@spatialys.com>> wrote:
On vendredi 13 novembre 2020 09:05:41 CET Sean Gillies via groups.io<https://urldefense.com/v3/__http://groups.io__;!!DOxrgLBm!U8N393pmVlI170WWlMCwvmrIV7A8JSl9WHWOD7E3re7ET8O3tNpqP7DhM2SV30dTeRRzVSk$> wrote:
HI Loïc,

I don't think we currently have any option other than to remove the
overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but
the implementation isn't there. For example, if you used gdalwarp with
`-ovr NONE` it will set an internal overview level to the value of -1. But
that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass
`OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.
Was left as an exercice for a contributor:
https://github.com/OSGeo/gdal/pull/3181<https://urldefense.com/v3/__https://github.com/OSGeo/gdal/pull/3181__;!!DOxrgLBm!U8N393pmVlI170WWlMCwvmrIV7A8JSl9WHWOD7E3re7ET8O3tNpqP7DhM2SV30dTK8LUNCg$>

Even

Thanks!

--
Sean Gillies


Sampling from one single point location using sample()

Alexander Jüstel
 

Hello, 
 
I would like to sample the value of my raster at one single point location. I have passed the coordinates as an iterable to the sample() function. However, I am not able to extract the value from the generator. The sampling works fine for more than one point. However, I need it to work for one point only. I have tried using list(gen) or next(gen) to get the value but nothing works so far.
 
Thanks for your support
Cheers
Alex


Re: Is it possible to ignore existing overview when performing decimated read?

Sean Gillies
 

Hi Even,

On Fri, Nov 13, 2020 at 10:38 AM Even Rouault <even.rouault@...> wrote:
On vendredi 13 novembre 2020 09:05:41 CET Sean Gillies via groups.io wrote:
> HI Loïc,
>
> I don't think we currently have any option other than to remove the
> overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but
> the implementation isn't there. For example, if you used gdalwarp with
> `-ovr NONE` it will set an internal overview level to the value of -1. But
> that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass
> `OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.

Was left as an exercice for a contributor:
https://github.com/OSGeo/gdal/pull/3181

Even

Thanks!

--
Sean Gillies


Re: Is it possible to ignore existing overview when performing decimated read?

Even Rouault
 

On vendredi 13 novembre 2020 09:05:41 CET Sean Gillies via groups.io wrote:
HI Loïc,

I don't think we currently have any option other than to remove the
overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but
the implementation isn't there. For example, if you used gdalwarp with
`-ovr NONE` it will set an internal overview level to the value of -1. But
that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass
`OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.
Was left as an exercice for a contributor:
https://github.com/OSGeo/gdal/pull/3181

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com


Re: Is it possible to ignore existing overview when performing decimated read?

Sean Gillies
 

HI Loïc,

I don't think we currently have any option other than to remove the overviews. In theory GDAL's OVERVIEW_LEVEL open option could help us, but the implementation isn't there. For example, if you used gdalwarp with `-ovr NONE` it will set an internal overview level to the value of -1. But that bypasses GDALOpenEx (used by rasterio). GDALOpenEx won't let us pass `OVERVIEW_LEVEL=-1. It fails with a "Cannot open overview level -1" error.

To illustrate

>>> rasterio.open("tests/data/green.tif").shape
(64, 64)
>>> rasterio.open("tests/data/green.tif", OVERVIEW_LEVEL=0).shape
(32, 32)
>>> rasterio.open("tests/data/green.tif", OVERVIEW_LEVEL="NONE").shape
(32, 32)
>>> rasterio.open("tests/data/green.tif", OVERVIEW_LEVEL=-1).shape
Traceback (most recent call last):
  File "rasterio/_base.pyx", line 261, in rasterio._base.DatasetBase.__init__
    self._hds = open_dataset(filename, flags, driver, kwargs, None)
  File "rasterio/_shim.pyx", line 78, in rasterio._shim.open_dataset
    return exc_wrap_pointer(hds)
  File "rasterio/_err.pyx", line 213, in rasterio._err.exc_wrap_pointer
    raise exc
rasterio._err.CPLE_OpenFailedError: Cannot open overview level -1 of tests/data/green.tif

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/sean/projects/rasterio/rasterio/env.py", line 433, in wrapper
    return f(*args, **kwds)
  File "/home/sean/projects/rasterio/rasterio/__init__.py", line 221, in open
    s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
  File "rasterio/_base.pyx", line 263, in rasterio._base.DatasetBase.__init__
    raise RasterioIOError(str(err))
rasterio.errors.RasterioIOError: Cannot open overview level -1 of tests/data/green.tif


On Fri, Nov 13, 2020 at 4:33 AM Loïc Dutrieux <loic.dutrieux@...> wrote:
Hi Denis,

Thanks for pointing to that issue. I may have misunderstood it but I don't think it helps answering my question.

Cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Denis Rykov <rykovd@...>
Sent: 10 November 2020 21:36:13
To: main@rasterio.groups.io
Subject: Re: [rasterio] Is it possible to ignore existing overview when performing decimated read?

You can find an example how to do that here: https://github.com/mapbox/rasterio/issues/1929<https://urldefense.com/v3/__https://github.com/mapbox/rasterio/issues/1929__;!!DOxrgLBm!XaInc3_J0VCLXQo6gYXo0rBWHs4b0p4ZkMuimU2g9HtJy0eRb55BN-J7JDPt-MB3Dhg5K3E$>.

On Tue, Nov 10, 2020, 8:12 PM Loïc Dutrieux <loic.dutrieux@...<mailto:loic.dutrieux@...>> wrote:
Hi everyone,

I'm trying to perform a decimated read of a dataset that already contains overviews. Regardless of which value I pass to the resampling= argument of the read method, it seems that the already existing overview is used. Is there any way to ignore it? See the reproducible example below.

Cheers,
Loïc

---

import tempfile
import os

import numpy as np
import rasterio
from rasterio.enums import Resampling


filename = os.path.join(tempfile.gettempdir(), 'overview_test.tif')

# Create random int array
shape = (100, 100)
arr = np.random.randint(0, 9999, size=shape, dtype=np.uint16)

meta = {'height': 100,
        'width': 100,
        'driver': 'GTiff',
        'dtype': np.uint16,
        'count': 1}

# Write random array to file and compute first overview
with rasterio.open(filename, 'w', **meta) as dst:
    dst.write(arr, 1)
    dst.build_overviews([2], Resampling.nearest)

# Read with downsample
with rasterio.open(filename) as src:
    arr_avg = src.read(1, out_shape=(1,50,50), resampling=Resampling.average,
                                  out_dtype=np.float)
    arr_nrt = src.read(1, out_shape=(1,50,50), resampling=Resampling.nearest)

print(np.max(arr_avg - arr_nrt))
# when the source file does not contain overviews, the max of the difference array
# is > 0


--
Sean Gillies


Re: Is it possible to ignore existing overview when performing decimated read?

Loïc Dutrieux
 

Hi Denis,

Thanks for pointing to that issue. I may have misunderstood it but I don't think it helps answering my question.

Cheers,
Loïc
________________________________________
From: main@rasterio.groups.io <main@rasterio.groups.io> on behalf of Denis Rykov <rykovd@gmail.com>
Sent: 10 November 2020 21:36:13
To: main@rasterio.groups.io
Subject: Re: [rasterio] Is it possible to ignore existing overview when performing decimated read?

You can find an example how to do that here: https://github.com/mapbox/rasterio/issues/1929<https://urldefense.com/v3/__https://github.com/mapbox/rasterio/issues/1929__;!!DOxrgLBm!XaInc3_J0VCLXQo6gYXo0rBWHs4b0p4ZkMuimU2g9HtJy0eRb55BN-J7JDPt-MB3Dhg5K3E$>.

On Tue, Nov 10, 2020, 8:12 PM Loïc Dutrieux <loic.dutrieux@ec.europa.eu<mailto:loic.dutrieux@ec.europa.eu>> wrote:
Hi everyone,

I'm trying to perform a decimated read of a dataset that already contains overviews. Regardless of which value I pass to the resampling= argument of the read method, it seems that the already existing overview is used. Is there any way to ignore it? See the reproducible example below.

Cheers,
Loïc

---

import tempfile
import os

import numpy as np
import rasterio
from rasterio.enums import Resampling


filename = os.path.join(tempfile.gettempdir(), 'overview_test.tif')

# Create random int array
shape = (100, 100)
arr = np.random.randint(0, 9999, size=shape, dtype=np.uint16)

meta = {'height': 100,
'width': 100,
'driver': 'GTiff',
'dtype': np.uint16,
'count': 1}

# Write random array to file and compute first overview
with rasterio.open(filename, 'w', **meta) as dst:
dst.write(arr, 1)
dst.build_overviews([2], Resampling.nearest)

# Read with downsample
with rasterio.open(filename) as src:
arr_avg = src.read(1, out_shape=(1,50,50), resampling=Resampling.average,
out_dtype=np.float)
arr_nrt = src.read(1, out_shape=(1,50,50), resampling=Resampling.nearest)

print(np.max(arr_avg - arr_nrt))
# when the source file does not contain overviews, the max of the difference array
# is > 0


Re: Rasterizing polygon from numpy arrays

alexisshakas@...
 

Works like a charm, thank you!


Re: Rasterizing polygon from numpy arrays

Brendan Ward
 

Yes, you should be able to pass a list / tuple (these are iterables) of Shapely geometries as input to geometry_mask.


Re: Rasterizing polygon from numpy arrays

alexisshakas@...
 

Hi Brendan and thnk you for the quick response!

The coordinates should not be an issue, since the polygon coordinates are in the same frame as the image.

Can i directly input a Shapely polygon to the geometry_mask class? I was a bit lost on whether this needs to be in some other format, like a sort of iterable dictionary.

121 - 140 of 780