Date   
Re: issue with opening/closing datasets

Amine Aboufirass
 

Hi Sean, the issue is that I am writing functions where the output is often a rasterio dataset. I don't know if this can be accomplished using a with statement:
function open_raster(filename):
    rasterio_dataset_object = rasterio.open(filename)
    return rasterio_dataset_object

function do_stuff_with_raster(rasterio_dataset_object):
    ###do stuff with raster
    return rasterio_dataset_object

dataset = open_raster('raster.tif')
new_raster = do_stuff_with_raster(dataset)
new_raster.close() 

Thanks,
Amine

On Fri, Mar 22, 2019 at 5:25 PM Sean Gillies <sean.gillies@...> wrote:
Hi Amine,

I think you have made in error in pasting code into the GitHub issue. The code you've given will fail at dataset = memfile.open because you haven't assigned memfile yet.

The message you see printed comes straight from the GDAL library. You haven't configured any GDAL error or log message handler and so the messages go directly to your terminal. Message handlers are configured if you run your code within a `with rasterio.Env()` block.

    import rasterio

    with rasterio.Env():
        # your code here

Also if you do

    with memfile.open(...) as dataset:

you won't see this message.


On Fri, Mar 22, 2019 at 9:25 AM Amine Aboufirass <amine.aboufirass@...> wrote:
Hi All, 

I just listed an issue on the main github log. https://github.com/mapbox/rasterio/issues/1659

If anyone could take a look I would be very grateful.

Kind Regards,

Amine



--
Sean Gillies

Re: Cloud-optimized GDAL VRT

mail@...
 

Huh, interesting. So only birds-eye overviews of VRTs should be expensive, which is to expected. Did you observe that in your tests? Did you ever go further and experiment with external overviews for the VRT?

Best,
Dion

Re: Cloud-optimized GDAL VRT

vincent.sarago@...
 

Hi Dion,
I encountered the same result couple month ago and couldn't get it to work. I asked Evan about it and see his answer here https://lists.osgeo.org/pipermail/gdal-dev/2018-October/049200.html

Vincent 

Cloud-optimized GDAL VRT

mail@...
 

We were discussing whether it would be possible to use `gdalbuildvrt` on several cloud-optimized GeoTIFF and obtain something that behaves like one big COG (which would be awesome to serve up large collections of GeoTIFFS, e.g. via Terracotta). `rasterio` seems to be able to read GDAL VRTs just fine, but it does not seem to have access to the overviews of the datasets. Is there anything obvious we could try? Is this missing functionality in GDAL or Rasterio?

See also https://github.com/DHI-GRAS/terracotta/issues/10

Re: issue with opening/closing datasets

Sean Gillies
 

Hi Amine,

I think you have made in error in pasting code into the GitHub issue. The code you've given will fail at dataset = memfile.open because you haven't assigned memfile yet.

The message you see printed comes straight from the GDAL library. You haven't configured any GDAL error or log message handler and so the messages go directly to your terminal. Message handlers are configured if you run your code within a `with rasterio.Env()` block.

    import rasterio

    with rasterio.Env():
        # your code here

Also if you do

    with memfile.open(...) as dataset:

you won't see this message.


On Fri, Mar 22, 2019 at 9:25 AM Amine Aboufirass <amine.aboufirass@...> wrote:
Hi All, 

I just listed an issue on the main github log. https://github.com/mapbox/rasterio/issues/1659

If anyone could take a look I would be very grateful.

Kind Regards,

Amine



--
Sean Gillies

Re: Output to grib2?

Shane Mill - NOAA Affiliate
 

Ah I see, I was confused then just because these come up when you "rio insp" and then "src.tags(1)", these creation options show up under the tags for band1.

That's very helpful so thanks for the response!

Shane

On Fri, Mar 22, 2019 at 12:09 PM Sean Gillies <sean.gillies@...> wrote:
Hi Shane,

I think I understand now. BAND_1_PDS_PDTN isn't a tag, it's a "creation option". You can apply these at the time you open a dataset for writing like

with rasterio.open("file.grb2", "w", driver="GRIB", ...., BAND_1_PDS_PDTN=something) as dataset:
    dataset.write(data)

On Fri, Mar 22, 2019 at 9:06 AM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
Okay so I think that I can update the GRIB specific options using rio warp

ie: rio warp 1000_500_thick.grb 1000_500_thick_new2.grb --format GRIB --profile BAND_1_PDS_PDTN=

with BAND_1_PDS_PDTN being a GRIB specific option. I'm going to explore this more today.

Thanks,
Shane

On Tue, Mar 19, 2019 at 3:50 PM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
Hey Sean,

I know that you were saying that you don't use grib2 often, but I wanted to follow up with the previous question.

I am able to output the numpy array to GRIB. Say that you have one band, you can read the tags of that one band with src.tags(1). These tags will be the defaults that come from:
https://gdal.org/frmt_grib.html (if you scroll down to the section "GRIB2 Write Support"). The problem with this is that for the grib2 to be readable, you need to be able to update these tags.

I'm wondering how to change the tags either when originally creating the GRIB or after the fact when using "update_tags" explained here: https://github.com/mapbox/rasterio/blob/master/docs/topics/tags.rst
It appears that update_tags works fine for GTiff but not for GRIB. 

Let me know if you have any thoughts on this. Thanks as always!

-Shane

On Fri, Mar 15, 2019 at 2:16 PM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
On Fri, Mar 15, 2019 at 11:02 AM, Sean Gillies wrote:
rio convert -f GRIB ~/code/rasterio/tests/data/RGB.byte.tif /tmp/rgb.grb2


Sean, thanks for the quick response! You are correct, that does work. I think my issue may have been that I was using gdalinfo to check the resulting grib2 file, and my gdal version was older. When I check with rio info it seems to be correct.

Thanks again!

-Shane



--
Shane Mill

Meteorological Application Developer, AceInfo Solutions

Meteorological Development Laboratory (NWS)

Phone: 301-427-9452



--
Shane Mill

Meteorological Application Developer, AceInfo Solutions

Meteorological Development Laboratory (NWS)

Phone: 301-427-9452



--
Sean Gillies



--
Shane Mill

Meteorological Application Developer, AceInfo Solutions

Meteorological Development Laboratory (NWS)

Phone: 301-427-9452

Re: Output to grib2?

Sean Gillies
 

Hi Shane,

I think I understand now. BAND_1_PDS_PDTN isn't a tag, it's a "creation option". You can apply these at the time you open a dataset for writing like

with rasterio.open("file.grb2", "w", driver="GRIB", ...., BAND_1_PDS_PDTN=something) as dataset:
    dataset.write(data)

On Fri, Mar 22, 2019 at 9:06 AM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
Okay so I think that I can update the GRIB specific options using rio warp

ie: rio warp 1000_500_thick.grb 1000_500_thick_new2.grb --format GRIB --profile BAND_1_PDS_PDTN=

with BAND_1_PDS_PDTN being a GRIB specific option. I'm going to explore this more today.

Thanks,
Shane

On Tue, Mar 19, 2019 at 3:50 PM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
Hey Sean,

I know that you were saying that you don't use grib2 often, but I wanted to follow up with the previous question.

I am able to output the numpy array to GRIB. Say that you have one band, you can read the tags of that one band with src.tags(1). These tags will be the defaults that come from:
https://gdal.org/frmt_grib.html (if you scroll down to the section "GRIB2 Write Support"). The problem with this is that for the grib2 to be readable, you need to be able to update these tags.

I'm wondering how to change the tags either when originally creating the GRIB or after the fact when using "update_tags" explained here: https://github.com/mapbox/rasterio/blob/master/docs/topics/tags.rst
It appears that update_tags works fine for GTiff but not for GRIB. 

Let me know if you have any thoughts on this. Thanks as always!

-Shane

On Fri, Mar 15, 2019 at 2:16 PM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
On Fri, Mar 15, 2019 at 11:02 AM, Sean Gillies wrote:
rio convert -f GRIB ~/code/rasterio/tests/data/RGB.byte.tif /tmp/rgb.grb2


Sean, thanks for the quick response! You are correct, that does work. I think my issue may have been that I was using gdalinfo to check the resulting grib2 file, and my gdal version was older. When I check with rio info it seems to be correct.

Thanks again!

-Shane



--
Shane Mill

Meteorological Application Developer, AceInfo Solutions

Meteorological Development Laboratory (NWS)

Phone: 301-427-9452



--
Shane Mill

Meteorological Application Developer, AceInfo Solutions

Meteorological Development Laboratory (NWS)

Phone: 301-427-9452



--
Sean Gillies

issue with opening/closing datasets

Amine Aboufirass
 

Hi All, 

I just listed an issue on the main github log. https://github.com/mapbox/rasterio/issues/1659

If anyone could take a look I would be very grateful.

Kind Regards,

Amine

Re: Output to grib2?

Shane Mill - NOAA Affiliate
 

Okay so I think that I can update the GRIB specific options using rio warp

ie: rio warp 1000_500_thick.grb 1000_500_thick_new2.grb --format GRIB --profile BAND_1_PDS_PDTN=

with BAND_1_PDS_PDTN being a GRIB specific option. I'm going to explore this more today.

Thanks,
Shane


On Tue, Mar 19, 2019 at 3:50 PM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
Hey Sean,

I know that you were saying that you don't use grib2 often, but I wanted to follow up with the previous question.

I am able to output the numpy array to GRIB. Say that you have one band, you can read the tags of that one band with src.tags(1). These tags will be the defaults that come from:
https://gdal.org/frmt_grib.html (if you scroll down to the section "GRIB2 Write Support"). The problem with this is that for the grib2 to be readable, you need to be able to update these tags.

I'm wondering how to change the tags either when originally creating the GRIB or after the fact when using "update_tags" explained here: https://github.com/mapbox/rasterio/blob/master/docs/topics/tags.rst
It appears that update_tags works fine for GTiff but not for GRIB. 

Let me know if you have any thoughts on this. Thanks as always!

-Shane

On Fri, Mar 15, 2019 at 2:16 PM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
On Fri, Mar 15, 2019 at 11:02 AM, Sean Gillies wrote:
rio convert -f GRIB ~/code/rasterio/tests/data/RGB.byte.tif /tmp/rgb.grb2


Sean, thanks for the quick response! You are correct, that does work. I think my issue may have been that I was using gdalinfo to check the resulting grib2 file, and my gdal version was older. When I check with rio info it seems to be correct.

Thanks again!

-Shane



--
Shane Mill

Meteorological Application Developer, AceInfo Solutions

Meteorological Development Laboratory (NWS)

Phone: 301-427-9452



--
Shane Mill

Meteorological Application Developer, AceInfo Solutions

Meteorological Development Laboratory (NWS)

Phone: 301-427-9452

Rasterio 1.0.22

Sean Gillies
 

Hi all,

Rasterio wheels and sdist for 1.0.22 are on PyPI today. This release fixes two bugs reported in 1.0.21.

Changes
=======

1.0.22 (2019-03-20)
-------------------

- Add JPEG2000 to enums.Compression so that the compression of JP2 files can be
  reported (#1654).
- Remove mock import from compat and move to test code (#1651).

--
Sean Gillies

Re: Python 2/3 compatibility management in rasterio

adrienoyono@...
 

Awesome ! Many thanks !! I will stay tuned to the release of 1.0.22 then

Re: Output to grib2?

Shane Mill - NOAA Affiliate
 

Hey Sean,

I know that you were saying that you don't use grib2 often, but I wanted to follow up with the previous question.

I am able to output the numpy array to GRIB. Say that you have one band, you can read the tags of that one band with src.tags(1). These tags will be the defaults that come from:
https://gdal.org/frmt_grib.html (if you scroll down to the section "GRIB2 Write Support"). The problem with this is that for the grib2 to be readable, you need to be able to update these tags.

I'm wondering how to change the tags either when originally creating the GRIB or after the fact when using "update_tags" explained here: https://github.com/mapbox/rasterio/blob/master/docs/topics/tags.rst
It appears that update_tags works fine for GTiff but not for GRIB. 

Let me know if you have any thoughts on this. Thanks as always!

-Shane

On Fri, Mar 15, 2019 at 2:16 PM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
On Fri, Mar 15, 2019 at 11:02 AM, Sean Gillies wrote:
rio convert -f GRIB ~/code/rasterio/tests/data/RGB.byte.tif /tmp/rgb.grb2


Sean, thanks for the quick response! You are correct, that does work. I think my issue may have been that I was using gdalinfo to check the resulting grib2 file, and my gdal version was older. When I check with rio info it seems to be correct.

Thanks again!

-Shane



--
Shane Mill

Meteorological Application Developer, AceInfo Solutions

Meteorological Development Laboratory (NWS)

Phone: 301-427-9452

Re: Python 2/3 compatibility management in rasterio

Sean Gillies
 

Hi,

We've always intended to keeping testing dependencies separate and the recent inclusion of mock in runtime dependencies was a mistake. It's been corrected, but it will be a few more days until we have a 1.0.22 release. I'm sorry about the inconvenience!


On Tue, Mar 19, 2019 at 8:21 AM <adrienoyono@...> wrote:
Hello everyone,

First of all, this is my occasion to say thanks for the rasterio project. It has been very helpful simplifying my life last year. I write a sofware that depends on rasterio, should ensure compatibility with Python 2 and Python 3, and should be deployed behind a web server (Apache with mod_wsgi). This last requirement triggers a problem caused by the way rasterio ensures Python 2/3 compatibility as I understand it after reading the "compat" module and the "setup.py".

Below is the stacktrace I have from Apache when trying to load my app:

Traceback (most recent call last):
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/rest/server.wsgi", line 3, in <module>
  from eodag.rest.server import app as application
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/__init__.py", line 18, in <module>
  from .api.core import EODataAccessGateway      # noqa
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/api/core.py", line 27, in <module>
  from eodag.api.search_result import SearchResult
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/api/search_result.py", line 26, in <module>
  from eodag.api.product import EOProduct
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/api/product/__init__.py", line 18, in <module>
  from ._product import EOProduct  # noqa
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/api/product/_product.py", line 25, in <module>
  import rasterio
File "/opt/eodag/venv/lib/python2.7/site-packages/rasterio/__init__.py", line 22, in <module>
  from rasterio._base import gdal_version
File "rasterio/_base.pyx", line 20, in init rasterio._base
File "/opt/eodag/venv/lib/python2.7/site-packages/rasterio/compat.py", line 26, in <module>
  import mock
ImportError: No module named mock
I guess "mock" is used in rasterio for testing purpose, I do use it for that myself. If that's the case, it is weird to have rasterio brake if we are in Python 2 without the "mock" module being installed, which is clearly the case here. To have mock installed at the same time as rasterio, one must do: "pip install rasterio[dev]". However, we always just do "pip install rasterio" (which fails to install "mock").

I think compatibility for testing purposes should be decoupled from compatibility for the core functionality. Is it possible to have it done in rasterio ? That would help avoid this kind of problem.

Thank you in advance for your time on this issue.



--
Sean Gillies

Python 2/3 compatibility management in rasterio

adrienoyono@...
 

Hello everyone,

First of all, this is my occasion to say thanks for the rasterio project. It has been very helpful simplifying my life last year. I write a sofware that depends on rasterio, should ensure compatibility with Python 2 and Python 3, and should be deployed behind a web server (Apache with mod_wsgi). This last requirement triggers a problem caused by the way rasterio ensures Python 2/3 compatibility as I understand it after reading the "compat" module and the "setup.py".

Below is the stacktrace I have from Apache when trying to load my app:

Traceback (most recent call last):
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/rest/server.wsgi", line 3, in <module>
  from eodag.rest.server import app as application
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/__init__.py", line 18, in <module>
  from .api.core import EODataAccessGateway      # noqa
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/api/core.py", line 27, in <module>
  from eodag.api.search_result import SearchResult
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/api/search_result.py", line 26, in <module>
  from eodag.api.product import EOProduct
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/api/product/__init__.py", line 18, in <module>
  from ._product import EOProduct  # noqa
File "/opt/eodag/venv/lib/python2.7/site-packages/eodag/api/product/_product.py", line 25, in <module>
  import rasterio
File "/opt/eodag/venv/lib/python2.7/site-packages/rasterio/__init__.py", line 22, in <module>
  from rasterio._base import gdal_version
File "rasterio/_base.pyx", line 20, in init rasterio._base
File "/opt/eodag/venv/lib/python2.7/site-packages/rasterio/compat.py", line 26, in <module>
  import mock
ImportError: No module named mock
I guess "mock" is used in rasterio for testing purpose, I do use it for that myself. If that's the case, it is weird to have rasterio brake if we are in Python 2 without the "mock" module being installed, which is clearly the case here. To have mock installed at the same time as rasterio, one must do: "pip install rasterio[dev]". However, we always just do "pip install rasterio" (which fails to install "mock").

I think compatibility for testing purposes should be decoupled from compatibility for the core functionality. Is it possible to have it done in rasterio ? That would help avoid this kind of problem.

Thank you in advance for your time on this issue.

Re: Output to grib2?

Shane Mill - NOAA Affiliate
 

On Fri, Mar 15, 2019 at 11:02 AM, Sean Gillies wrote:
rio convert -f GRIB ~/code/rasterio/tests/data/RGB.byte.tif /tmp/rgb.grb2


Sean, thanks for the quick response! You are correct, that does work. I think my issue may have been that I was using gdalinfo to check the resulting grib2 file, and my gdal version was older. When I check with rio info it seems to be correct.

Thanks again!

-Shane

Re: Output to grib2?

Sean Gillies
 

Hi Shane,

Rasterio's rio-convert command can convert from GeoTIFF to GRIB, so it should be possible

$ rio convert -f GRIB ~/code/rasterio/tests/data/RGB.byte.tif /tmp/rgb.grb2
$ rio info /tmp/rgb.grb2
{"bounds": [101985.0, 2611486.28, 339316.64, 2826915.0], "colorinterp": ["undefined", "undefined", "undefined"], "count": 3, "crs": "PROJCS[\"unnamed\",GEOGCS[\"Coordinate System imported from GRIB file\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-2072.483648],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"Metre\",1]]", "descriptions": ["0[-] RESERVED(0) (Reserved)", "0[-] RESERVED(0) (Reserved)", "0[-] RESERVED(0) (Reserved)"], "driver": "GRIB", "dtype": "float64", "height": 718, "indexes": [1, 2, 3], "lnglat": [84.75845370919801, 24.561589206961006], "mask_flags": [["nodata"], ["nodata"], ["nodata"]], "nodata": 0.0, "res": [300.04, 300.04], "shape": [718, 791], "tiled": false, "transform": [300.04, 0.0, 101985.0, 0.0, -300.04, 2826915.0, 0.0, 0.0, 1.0], "units": [null, null, null], "width": 791}

I'm not a regular GRIB user, but it looks to me that GDAL produces a GRIB2 file when you write using the GRIB format driver, and that the format driver reads either version without a need to specify. GRIB2 writing requires GDAL 2.3.0 or newer. The latest rasterio wheels on PyPI include GDAL 2.4.0, and should be sufficient.

On Fri, Mar 15, 2019 at 11:53 AM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
Is there anyway to output to grib2 after performing map algebra? For example, If you are performing a calculation between two grib2 bands, which rasterio does quite well, are you able to transform the numpy array into a grib2 file? I have been able to do this with geotiff but not grib2.



--
Sean Gillies

Output to grib2?

Shane Mill - NOAA Affiliate
 

Is there anyway to output to grib2 after performing map algebra? For example, If you are performing a calculation between two grib2 bands, which rasterio does quite well, are you able to transform the numpy array into a grib2 file? I have been able to do this with geotiff but not grib2.

Rasterio 1.0.20

Sean Gillies
 

Hi all,

Rasterio 1.0.20 is on PyPI today. I hope you'll like the changes: https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L4-L24.

I discovered that 1.0.19 was defective while building linux wheels to distribute. No distributions of 1.0.19 were uploaded.

--
Sean Gillies

Re: Implementing -scale functionality of gdal_translate

linden.kyle@...
 

That worked perfectly. Thanks for all your help.

Re: Implementing -scale functionality of gdal_translate

Sean Gillies
 

Hi,

On Wed, Feb 6, 2019 at 2:32 PM <linden.kyle@...> wrote:
So I'm seeing some odd behavior now. I ended up trying the masked_array solution, but I'm seeing breaks in my data at the tile boundaries.

Here's what the data looks like if I use reprojection with just the source band (code: https://pastebin.com/dZVRVKi2): https://imgur.com/WpFKfij

Here's what the data looks like I read into a masked array and use an ndarray source (code: https://pastebin.com/z1xd8mXu): https://imgur.com/2QkQQZf
Ironically enough I get the exact same result with my solution that doesn't use a masked array, but instead modifies the transform: https://pastebin.com/7ph1VYWX

Your code looks okay. I recommend padding your read window by a few pixels on each side so that tile_window, tile_data, and tile_transform are slight dilated. I think this is very likely to eliminate what looks like an off-by-one error in the reprojection, an error that you don't see when you use the entire source dataset as the reprojection source. Does that make sense?

--
Sean Gillies