Date   

Re: Rasterio 1.2b1

Alan Snow
 

Is that because on the other hand GDAL 3.2.0 and/or PROJ 7.2.0 are heavier than their counterparts from release 1.1.8?
The new wheel now includes the proj.db, which is close to the size of grids that were removed. Also, PROJ requires more dependencies, so that is also a likely factor.

would you expect there to be problems if I were to have pyproj<3.0.0 and rasterio>=1.2 installed in my Python environment, due to their different PROJ versions, or would that still work OK?
It might work, but I wouldn't recommend it. At the very least, there are different versions proj.db in each one. If you do go that route, I would recommend you proceed with caution.


Re: Rasterio 1.2b1

Guillaume Lostis <g.lostis@...>
 

Hi Sean,

Thanks a lot for this upcoming release, I am really looking forward to using some of its new features!

Out of interest, the 1.2b1 wheels do not seem to be lighter than the 1.1.8 wheels even though they do not contain PROJ data. Is that because on the other hand GDAL 3.2.0 and/or PROJ 7.2.0 are heavier than their counterparts from release 1.1.8?

Also Sean and Alan, would you expect there to be problems if I were to have pyproj<3.0.0 and rasterio>=1.2 installed in my Python environment, due to their different PROJ versions, or would that still work OK?

Best,

Guillaume

On Tue, Dec 15, 2020 at 6:59 PM Alan Snow <alansnow21@...> wrote:
Thanks for working on this Sean!

I verified that PROJ_NETWORK=ON works with a basic test:

rio_geom.py:
from rasterio.warp import transform_geom

geometry = [
{
"type": "Polygon",
"coordinates": [
[
[-94.07955380199459, 41.69085871273774],
[-94.06082436942204, 41.69103313774798],
[-94.06063203899649, 41.67932439500822],
[-94.07935807746362, 41.679150041277325],
[-94.07955380199459, 41.69085871273774],
]
],
}
]


print(transform_geom("+proj=latlon", "+proj=utm +zone=10 +datum=NAD27", geometry))

$ python rio_geom.py
[{'type': 'Polygon', 'coordinates': [[(2914912.501340564, 5039833.18814224), (2916474.3810722474, 5040429.817378172), (2916971.9920050735, 5039127.068439782), (2915410.0220987815, 5038530.542053506), (2914912.501340564, 5039833.18814224)]]}]
$ PROJ_NETWORK=ON python rio_geom.py
[{'type': 'Polygon', 'coordinates': [[(2914912.2992929644, 5039834.208342302), (2916474.1765534407, 5040430.837641376), (2916971.7831669645, 5039128.087552843), (2915409.815547156, 5038531.561302773), (2914912.2992929644, 5039834.208342302)]]}]


Re: rasterio.features.shapes with holes in polygons

Sean Gillies
 

Hi Paolo,

Welcome! I ran the your file through rio-shapes and jq

rio shapes --bidx 1 ~/Downloads/sample_cover.tif | jq -c 'select(.properties.val == 3.0)' | fio collect

and then uploaded to a gist:


It looks like the holes aren't lost there. How did you make the second image?

On Tue, Dec 15, 2020 at 3:47 PM Paolo Corti <pcorti@...> wrote:
Hello users and devs

First: thanks a lot for this wonderful project, I am really enjoying using it for my geospatial needs.

For a specific process I am using rasterio.features.shapes and it is working greatly for my purpose. However I have a few cases where I am experiencing strange results.
For example, see this 1 band raster: 

Screen Shot 2020-12-15 at 5.33.13 PM.png
I would like to extract the polygons for the yellow area (value: 3).

When running the following code:

with rasterio.open(cover_path) as src:
cover_arr = src.read(1)
mask = (cover_arr == 3)
polygons = shapes(cover_arr, mask=mask, transform=src.transform)

I get these results (in blue):

Screen Shot 2020-12-15 at 5.35.37 PM.png

Apparently holes are not correctly parsed.

It would be great to know if I am doing something wrong here or if it is something which is not working properly.

For your reference, I am attaching a copy of the raster I am using here so you can reproduce this (sample_cover.tif)

Thanks in advance!
Paolo

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1



--
Sean Gillies


rasterio.features.shapes with holes in polygons

Paolo Corti
 

Hello users and devs

First: thanks a lot for this wonderful project, I am really enjoying using it for my geospatial needs.

For a specific process I am using rasterio.features.shapes and it is working greatly for my purpose. However I have a few cases where I am experiencing strange results.
For example, see this 1 band raster: 

Screen Shot 2020-12-15 at 5.33.13 PM.png
I would like to extract the polygons for the yellow area (value: 3).

When running the following code:

with rasterio.open(cover_path) as src:
cover_arr = src.read(1)
mask = (cover_arr == 3)
polygons = shapes(cover_arr, mask=mask, transform=src.transform)

I get these results (in blue):

Screen Shot 2020-12-15 at 5.35.37 PM.png

Apparently holes are not correctly parsed.

It would be great to know if I am doing something wrong here or if it is something which is not working properly.

For your reference, I am attaching a copy of the raster I am using here so you can reproduce this (sample_cover.tif)

Thanks in advance!
Paolo

--
Paolo Corti
Geospatial software developer
web: http://www.paolocorti.net
twitter: @capooti
skype: capooti
#drt3jc1


Re: Rasterio 1.2b1

Alan Snow
 

Thanks for working on this Sean!

I verified that PROJ_NETWORK=ON works with a basic test:

rio_geom.py:
from rasterio.warp import transform_geom

geometry = [
{
"type": "Polygon",
"coordinates": [
[
[-94.07955380199459, 41.69085871273774],
[-94.06082436942204, 41.69103313774798],
[-94.06063203899649, 41.67932439500822],
[-94.07935807746362, 41.679150041277325],
[-94.07955380199459, 41.69085871273774],
]
],
}
]


print(transform_geom("+proj=latlon", "+proj=utm +zone=10 +datum=NAD27", geometry))

$ python rio_geom.py
[{'type': 'Polygon', 'coordinates': [[(2914912.501340564, 5039833.18814224), (2916474.3810722474, 5040429.817378172), (2916971.9920050735, 5039127.068439782), (2915410.0220987815, 5038530.542053506), (2914912.501340564, 5039833.18814224)]]}]
$ PROJ_NETWORK=ON python rio_geom.py
[{'type': 'Polygon', 'coordinates': [[(2914912.2992929644, 5039834.208342302), (2916474.1765534407, 5040430.837641376), (2916971.7831669645, 5039128.087552843), (2915409.815547156, 5038531.561302773), (2914912.2992929644, 5039834.208342302)]]}]


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

281 - 300 of 945