Date   
How to extract DATA_ENCODING creation option type from an input grib2 file for use in the output grib2 file

Shane Mill - NOAA Affiliate
 

Hey everyone,

I know that grib2 isn't commonly used by this community so I apologize that I keep bringing it up... You can probably tell that I am super fond of it! Anyways, all kidding aside, I am very impressed with the abilities of Rasterio and how responsive this community has been.

In a previous topic, Sean and I discussed that creation options for a specific driver can be assigned when opening a grib file in writing mode.
ie. "with rasterio.open(file.grb2, 'w', driver='GRIB',dtype='float64'.....DATA_ENCODING='COMPLEX_PACKING',SPATIAL_DIFFERENCING_ORDER=2) as ds:"

The grib2 specific creation options can be found here:
https://www.gdal.org/frmt_grib.html

This works well in terms of creating a grib2 file encoded in a way that can be read by gdalinfo, toolsUI, and other grib2 readers. What I want/need to figure out now is how to be able to extract the DATA_ENCODING and SPATIAL_DIFFERENCING_ORDER from an input grib2 file so that these creation options don't need to be hardcoded when creating a new grib2 file. In other words, say that you have model data in an input grib2 file, and you create a new grib2 file with a new derived parameter that is derived from two bands in the input grib2 file. I would want the DATA_ENCODING to be the same for both the input and output grib2 file.

Many of the creation options are tracked in the metadata, but as far as I can tell, DATA_ENCODING is not. I'm wondering if it actually is somewhere and I'm just missing it. 

At this link (https://github.com/mapbox/rasterio/issues/405), Sean, you have indicated that the GDAL API doesn't track creation options, but you were able to develop this in Rasterio. Going back to https://www.gdal.org/frmt_grib.html, It looks like all of the Product Identification and Definition is either tracked in src.meta, src.profile, or src.tags(band#). It doesn't appear that Data Encoding is.

I just wanted to see if you (or anyone else) had any feedback on this or know if there is anyway of extracting the type of data encoding from an input file.

Thanks as always!
Shane



Re: How to extract DATA_ENCODING creation option type from an input grib2 file for use in the output grib2 file

Sean Gillies
 

Hi Shane,

Before Rasterio 1.0, the project did record creation options in a custom tag (metadata) namespace. We discovered problems when converting from GeoTIFF to other formats: creation options are very format specific and few of the GeoTIFF ones make any sense for other formats. It turned out that if you, for example, translated from GeoTIFF to something else and then back to GeoTIFF, you could end up with recorded creation options that didn't reflect what existed in the data. We decided to remove this feature in Rasterio 1.0 and leave it up to developers to do it for themselves if they need it.

I see what you mean about the GRIB driver's DATA_ENCODING option. I don't see a good way to get this from the GDAL API. You might have to use a native GRIB API or begin to store the value of this creation option in your application or in custom namespace in your dataset metadata.


On Tue, Apr 2, 2019 at 3:22 PM Shane Mill - NOAA Affiliate via Groups.Io <shane.mill=noaa.gov@groups.io> wrote:
Hey everyone,

I know that grib2 isn't commonly used by this community so I apologize that I keep bringing it up... You can probably tell that I am super fond of it! Anyways, all kidding aside, I am very impressed with the abilities of Rasterio and how responsive this community has been.

In a previous topic, Sean and I discussed that creation options for a specific driver can be assigned when opening a grib file in writing mode.
ie. "with rasterio.open(file.grb2, 'w', driver='GRIB',dtype='float64'.....DATA_ENCODING='COMPLEX_PACKING',SPATIAL_DIFFERENCING_ORDER=2) as ds:"

The grib2 specific creation options can be found here:
https://www.gdal.org/frmt_grib.html

This works well in terms of creating a grib2 file encoded in a way that can be read by gdalinfo, toolsUI, and other grib2 readers. What I want/need to figure out now is how to be able to extract the DATA_ENCODING and SPATIAL_DIFFERENCING_ORDER from an input grib2 file so that these creation options don't need to be hardcoded when creating a new grib2 file. In other words, say that you have model data in an input grib2 file, and you create a new grib2 file with a new derived parameter that is derived from two bands in the input grib2 file. I would want the DATA_ENCODING to be the same for both the input and output grib2 file.

Many of the creation options are tracked in the metadata, but as far as I can tell, DATA_ENCODING is not. I'm wondering if it actually is somewhere and I'm just missing it. 

At this link (https://github.com/mapbox/rasterio/issues/405), Sean, you have indicated that the GDAL API doesn't track creation options, but you were able to develop this in Rasterio. Going back to https://www.gdal.org/frmt_grib.html, It looks like all of the Product Identification and Definition is either tracked in src.meta, src.profile, or src.tags(band#). It doesn't appear that Data Encoding is.

I just wanted to see if you (or anyone else) had any feedback on this or know if there is anyway of extracting the type of data encoding from an input file.

Thanks as always!
Shane





--
Sean Gillies

Re: How to extract DATA_ENCODING creation option type from an input grib2 file for use in the output grib2 file

Shane Mill - NOAA Affiliate
 

Hey Sean,

Thanks for the response, I can definitely see why that would be problematic, especially for converting between formats. That makes perfect sense to me. It's unfortunate that it isn't originally provided in the GDAL API, but I'll look at other ways of getting around this with a native GRIB API or within the application.

Again, thanks as always for the context and background. It's very much appreciated.

Shane

Re: Rewriting uint16 headers with rasterio / applying rio color makes them unreadable by Preview, Photoshop

Edward Boyda
 


Thanks, Sean, that works for rio color.  I use the photometric creation option with gdal_translate but didn't connect the dots to rio color. 

I can see why, given all the possible creation options, this isn't mentioned explicitly in the docs. But maybe photometric interpretation makes a worthy special case - it's counterintuitive that an operation to adjust the color of an image doesn't by default preserve its header structure. 

Is there a similarly simple fix for setting colorinterp manually in update mode (my third example)? 

Much appreciated.
Ed


On Tue, Apr 2, 2019 at 1:55 PM Sean Gillies <sean.gillies@...> wrote:
Hi Ed,

Can you try the following variation on your first command?

$ rio color -j 1 uint16_image.tif uint16_brightened.tif gamma RGB 1.5 --co photometric=RGB

Note the addition of "--co photometric=RGB". GDAL automatically sets the photometric tag (which other apps need) to RGB when the created image data type is uint8 (see https://www.gdal.org/frmt_gtiff.html) but does not do the same for other data types including uint16.


On Thu, Mar 28, 2019 at 5:48 PM Edward Boyda <boyda@...> wrote:

Hi all, I'm new at this, more of a computer vision person than a developer, so please bear with me....

I see the behavior I'm about to describe running pip-installed rasterio (1.0.22) on my Mac (OSX Mojave; homebrewed python and gdal) and also running rasterio (1.0.22) in a dockerized Ubuntu platform, on images from a variety of sources (DigitalGlobe, Planet, Landsat). 

Example 1: 
$ rio color -j 1 uint16_image.tif uint16_brightened.tif gamma RGB 1.5

When I try to open the output file, uint16_brightened.tif, in Preview or Photoshop, I get a message like "Could not complete your request because of a problem parsing the TIFF file." (That's from Photoshop; Preview is equivalent.) 

Example 2:  
$ rio color -j -1 uint16_image.tif uint16_brightened_v2.tif gamma RGB 1.5

The number of cores has changed from the first example. Now the output, uint16_brightened_v2.tif, is readable by Photoshop but has had its color interpretation changed to (reading from rasterio):

(<ColorInterp.gray: 1>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>)

When I open the file with Preview or view the thumbnail with Mac Finder, there are dark vertical lines interspersed with the actual pixels, and about a third of the original pixels have been pushed out of the frame. See screenshot attached.

Example 3:
I take a file that has a (gray, undefined, undefined) color interpretation and try to change that to RGB, now in the interpreter:
>>> with rasterio.open('uint16_noCI.tif', 'r+') as f:
    f.colorinterp = (ColorInterp.red, ColorInterp.green, ColorInterp.blue)

Again the edited file is unreadable by Preview and Photoshop. 

A couple of caveats:

1) I can read the data from files output from any of the above examples with rasterio or skimage, and the resulting numpy array is uncorrputed. I can resave it with skimage, show it with matplolib, etc., and the image looks as expected.
2) With any of the above outputs, I can read and then rewrite a new_image.tif, using rasterio, and the resulting files open as expected with Photoshop and Preview. This is my current (obviously inefficient) workaround:

with rasterio.open('uint16_brightened.tif') as f:
    prof = f.profile
    img = f.read()

with rasterio.open('new_image.tif', 'w', photometric='rgb', **prof) as f:
    f.write(img)

As far as I know these failures happen only with uint16 images (at least not with uint8), and it would seem to have to do with the way color interpretation is written into the headers via the different write mechanisms.  Has anyone come across similar behavior? I've been reproducing and beating my head against this for months and would really appreciate a sanity check.

Since the Docker container is likely cleaner than what I've installed on my Mac, I've run the tests for the attached output there. Here is the Dockerfile and some possibly relevant release details: 

FROM ubuntu:latest

RUN apt-get update && apt-get install -y software-properties-common

RUN apt-get install -y python3-pip python3-dev build-essential

RUN pip3 install --upgrade pip

RUN apt-get install -y gdal-bin libgdal-dev python3-gdal

RUN apt-get install -y libssl-dev libffi-dev libcurl4-openssl-dev

ENV DEBIAN_FRONTEND noninteractive

RUN apt-get install -y python3-tk


ADD ./requirements.txt /tmp/requirements.txt

RUN pip install -r /tmp/requirements.txt


---

DISTRIB_ID=Ubuntu

DISTRIB_RELEASE=18.04

python 3.6.6

pip 18.1 from /usr/local/lib/python3.6/dist-packages/pip (python 3.6)


affine==2.2.2

gbdx-auth==0.4.0

gbdxtools==0.16.0

GDAL==2.2.2

rasterio==1.0.22

rio-color==1.0.0

rio-mucho==1.0.0

Shapely==1.6.4

tifffile==2018.11.28



Thanks everyone!


Ed



--
Sean Gillies

Re: Rewriting uint16 headers with rasterio / applying rio color makes them unreadable by Preview, Photoshop

Sean Gillies
 

For GeoTIFF (at least), the photometric creation option is different from the color interpretation that we can get/set using the GDAL API. The file's layout and compression strategy is influenced by the photometric creation option, so it is needed up front and has a permanent impact on the file. Setting the color interpretation through the rasterio API (and thereby GDAL API) will change file metadata but will not change the file's structure. It's confusing, to be sure. Even Rouault (the expert) provides a bit more about it in this email to gdal-dev:


On Wed, Apr 3, 2019 at 10:31 AM Edward Boyda <boyda@...> wrote:

Thanks, Sean, that works for rio color.  I use the photometric creation option with gdal_translate but didn't connect the dots to rio color. 

I can see why, given all the possible creation options, this isn't mentioned explicitly in the docs. But maybe photometric interpretation makes a worthy special case - it's counterintuitive that an operation to adjust the color of an image doesn't by default preserve its header structure. 

Is there a similarly simple fix for setting colorinterp manually in update mode (my third example)? 

Much appreciated.
Ed


On Tue, Apr 2, 2019 at 1:55 PM Sean Gillies <sean.gillies@...> wrote:
Hi Ed,

Can you try the following variation on your first command?

$ rio color -j 1 uint16_image.tif uint16_brightened.tif gamma RGB 1.5 --co photometric=RGB

Note the addition of "--co photometric=RGB". GDAL automatically sets the photometric tag (which other apps need) to RGB when the created image data type is uint8 (see https://www.gdal.org/frmt_gtiff.html) but does not do the same for other data types including uint16.


On Thu, Mar 28, 2019 at 5:48 PM Edward Boyda <boyda@...> wrote:

Hi all, I'm new at this, more of a computer vision person than a developer, so please bear with me....

I see the behavior I'm about to describe running pip-installed rasterio (1.0.22) on my Mac (OSX Mojave; homebrewed python and gdal) and also running rasterio (1.0.22) in a dockerized Ubuntu platform, on images from a variety of sources (DigitalGlobe, Planet, Landsat). 

Example 1: 
$ rio color -j 1 uint16_image.tif uint16_brightened.tif gamma RGB 1.5

When I try to open the output file, uint16_brightened.tif, in Preview or Photoshop, I get a message like "Could not complete your request because of a problem parsing the TIFF file." (That's from Photoshop; Preview is equivalent.) 

Example 2:  
$ rio color -j -1 uint16_image.tif uint16_brightened_v2.tif gamma RGB 1.5

The number of cores has changed from the first example. Now the output, uint16_brightened_v2.tif, is readable by Photoshop but has had its color interpretation changed to (reading from rasterio):

(<ColorInterp.gray: 1>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>)

When I open the file with Preview or view the thumbnail with Mac Finder, there are dark vertical lines interspersed with the actual pixels, and about a third of the original pixels have been pushed out of the frame. See screenshot attached.

Example 3:
I take a file that has a (gray, undefined, undefined) color interpretation and try to change that to RGB, now in the interpreter:
>>> with rasterio.open('uint16_noCI.tif', 'r+') as f:
    f.colorinterp = (ColorInterp.red, ColorInterp.green, ColorInterp.blue)

Again the edited file is unreadable by Preview and Photoshop. 

A couple of caveats:

1) I can read the data from files output from any of the above examples with rasterio or skimage, and the resulting numpy array is uncorrputed. I can resave it with skimage, show it with matplolib, etc., and the image looks as expected.
2) With any of the above outputs, I can read and then rewrite a new_image.tif, using rasterio, and the resulting files open as expected with Photoshop and Preview. This is my current (obviously inefficient) workaround:

with rasterio.open('uint16_brightened.tif') as f:
    prof = f.profile
    img = f.read()

with rasterio.open('new_image.tif', 'w', photometric='rgb', **prof) as f:
    f.write(img)

As far as I know these failures happen only with uint16 images (at least not with uint8), and it would seem to have to do with the way color interpretation is written into the headers via the different write mechanisms.  Has anyone come across similar behavior? I've been reproducing and beating my head against this for months and would really appreciate a sanity check.

Since the Docker container is likely cleaner than what I've installed on my Mac, I've run the tests for the attached output there. Here is the Dockerfile and some possibly relevant release details: 

FROM ubuntu:latest

RUN apt-get update && apt-get install -y software-properties-common

RUN apt-get install -y python3-pip python3-dev build-essential

RUN pip3 install --upgrade pip

RUN apt-get install -y gdal-bin libgdal-dev python3-gdal

RUN apt-get install -y libssl-dev libffi-dev libcurl4-openssl-dev

ENV DEBIAN_FRONTEND noninteractive

RUN apt-get install -y python3-tk


ADD ./requirements.txt /tmp/requirements.txt

RUN pip install -r /tmp/requirements.txt


---

DISTRIB_ID=Ubuntu

DISTRIB_RELEASE=18.04

python 3.6.6

pip 18.1 from /usr/local/lib/python3.6/dist-packages/pip (python 3.6)


affine==2.2.2

gbdx-auth==0.4.0

gbdxtools==0.16.0

GDAL==2.2.2

rasterio==1.0.22

rio-color==1.0.0

rio-mucho==1.0.0

Shapely==1.6.4

tifffile==2018.11.28



Thanks everyone!


Ed



--
Sean Gillies



--
Sean Gillies

Re: How to extract DATA_ENCODING creation option type from an input grib2 file for use in the output grib2 file

Shane Mill - NOAA Affiliate
 

In case anyone ends up needing to do this, I ended up using eccodes and you can retrieve and set 'packingType' which is the same as DATA_ENCODING with codes_get and codes_set

Thanks,
-Shane

Re: Rewriting uint16 headers with rasterio / applying rio color makes them unreadable by Preview, Photoshop

Edward Boyda
 


Got it, thanks!!

On Wed, Apr 3, 2019 at 1:12 PM Sean Gillies <sean.gillies@...> wrote:
For GeoTIFF (at least), the photometric creation option is different from the color interpretation that we can get/set using the GDAL API. The file's layout and compression strategy is influenced by the photometric creation option, so it is needed up front and has a permanent impact on the file. Setting the color interpretation through the rasterio API (and thereby GDAL API) will change file metadata but will not change the file's structure. It's confusing, to be sure. Even Rouault (the expert) provides a bit more about it in this email to gdal-dev:


On Wed, Apr 3, 2019 at 10:31 AM Edward Boyda <boyda@...> wrote:

Thanks, Sean, that works for rio color.  I use the photometric creation option with gdal_translate but didn't connect the dots to rio color. 

I can see why, given all the possible creation options, this isn't mentioned explicitly in the docs. But maybe photometric interpretation makes a worthy special case - it's counterintuitive that an operation to adjust the color of an image doesn't by default preserve its header structure. 

Is there a similarly simple fix for setting colorinterp manually in update mode (my third example)? 

Much appreciated.
Ed


On Tue, Apr 2, 2019 at 1:55 PM Sean Gillies <sean.gillies@...> wrote:
Hi Ed,

Can you try the following variation on your first command?

$ rio color -j 1 uint16_image.tif uint16_brightened.tif gamma RGB 1.5 --co photometric=RGB

Note the addition of "--co photometric=RGB". GDAL automatically sets the photometric tag (which other apps need) to RGB when the created image data type is uint8 (see https://www.gdal.org/frmt_gtiff.html) but does not do the same for other data types including uint16.


On Thu, Mar 28, 2019 at 5:48 PM Edward Boyda <boyda@...> wrote:

Hi all, I'm new at this, more of a computer vision person than a developer, so please bear with me....

I see the behavior I'm about to describe running pip-installed rasterio (1.0.22) on my Mac (OSX Mojave; homebrewed python and gdal) and also running rasterio (1.0.22) in a dockerized Ubuntu platform, on images from a variety of sources (DigitalGlobe, Planet, Landsat). 

Example 1: 
$ rio color -j 1 uint16_image.tif uint16_brightened.tif gamma RGB 1.5

When I try to open the output file, uint16_brightened.tif, in Preview or Photoshop, I get a message like "Could not complete your request because of a problem parsing the TIFF file." (That's from Photoshop; Preview is equivalent.) 

Example 2:  
$ rio color -j -1 uint16_image.tif uint16_brightened_v2.tif gamma RGB 1.5

The number of cores has changed from the first example. Now the output, uint16_brightened_v2.tif, is readable by Photoshop but has had its color interpretation changed to (reading from rasterio):

(<ColorInterp.gray: 1>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>)

When I open the file with Preview or view the thumbnail with Mac Finder, there are dark vertical lines interspersed with the actual pixels, and about a third of the original pixels have been pushed out of the frame. See screenshot attached.

Example 3:
I take a file that has a (gray, undefined, undefined) color interpretation and try to change that to RGB, now in the interpreter:
>>> with rasterio.open('uint16_noCI.tif', 'r+') as f:
    f.colorinterp = (ColorInterp.red, ColorInterp.green, ColorInterp.blue)

Again the edited file is unreadable by Preview and Photoshop. 

A couple of caveats:

1) I can read the data from files output from any of the above examples with rasterio or skimage, and the resulting numpy array is uncorrputed. I can resave it with skimage, show it with matplolib, etc., and the image looks as expected.
2) With any of the above outputs, I can read and then rewrite a new_image.tif, using rasterio, and the resulting files open as expected with Photoshop and Preview. This is my current (obviously inefficient) workaround:

with rasterio.open('uint16_brightened.tif') as f:
    prof = f.profile
    img = f.read()

with rasterio.open('new_image.tif', 'w', photometric='rgb', **prof) as f:
    f.write(img)

As far as I know these failures happen only with uint16 images (at least not with uint8), and it would seem to have to do with the way color interpretation is written into the headers via the different write mechanisms.  Has anyone come across similar behavior? I've been reproducing and beating my head against this for months and would really appreciate a sanity check.

Since the Docker container is likely cleaner than what I've installed on my Mac, I've run the tests for the attached output there. Here is the Dockerfile and some possibly relevant release details: 

FROM ubuntu:latest

RUN apt-get update && apt-get install -y software-properties-common

RUN apt-get install -y python3-pip python3-dev build-essential

RUN pip3 install --upgrade pip

RUN apt-get install -y gdal-bin libgdal-dev python3-gdal

RUN apt-get install -y libssl-dev libffi-dev libcurl4-openssl-dev

ENV DEBIAN_FRONTEND noninteractive

RUN apt-get install -y python3-tk


ADD ./requirements.txt /tmp/requirements.txt

RUN pip install -r /tmp/requirements.txt


---

DISTRIB_ID=Ubuntu

DISTRIB_RELEASE=18.04

python 3.6.6

pip 18.1 from /usr/local/lib/python3.6/dist-packages/pip (python 3.6)


affine==2.2.2

gbdx-auth==0.4.0

gbdxtools==0.16.0

GDAL==2.2.2

rasterio==1.0.22

rio-color==1.0.0

rio-mucho==1.0.0

Shapely==1.6.4

tifffile==2018.11.28



Thanks everyone!


Ed



--
Sean Gillies



--
Sean Gillies

Memory error in rio calc that equivalent gdal_calc.py performs quickly without issue

lawr@...
 

Today I tried to replace a `gdal_calc.py` command with the equivalent `rio calc`, because I wanted a UInt8 output, and `gdal_calc.py` restricts me to UInt16.

 

The GDAL command is relatively quick on my 15m² raster of New Zealand, at about 2-3 minutes processing time on my machine, but the equivalent `rio calc` that I came up with seemed to be attempting to load everything into memory, and the command eventually crashed with a Memory Error after struggling for about ten minutes.

Here's the command:

```

gdal_calc.py -A a.tif -B b.tif< \
--outfile=out.tif --overwrite --calc="where(A != 1,0,1) * B" \
--format=GTiff --NoDataValue=0 --type=UInt16\
--co TILED=YES --co COMPRESS=LZW --co TFW=YES \
--co BLOCKXSIZE=128 --co BLOCKYSIZE=128
```

gdalinfo for a.tif:

```
Size is 66814, 96417
Coordinate System is:
PROJCS["NZGD_2000_New_Zealand_Transverse_Mercator",
    GEOGCS["GCS_NZGD_2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (1089805.593673222931102,6194224.566799569875002)
Pixel Size = (15.000000000000000,-15.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 1089805.594, 6194224.567) (167d27'39.65"E, 34d16' 4.53"S)
Lower Left  ( 1089805.594, 4747969.567) (166d15'35.49"E, 47d13'23.08"S)
Upper Right ( 2092015.594, 6194224.567) (178d20'32.62"E, 34d16'36.05"S)
Lower Right ( 2092015.594, 4747969.567) (179d30' 5.76"E, 47d14'12.93"S)
Center      ( 1590910.594, 5471097.067) (172d53'31.45"E, 40d54'40.25"S)
Band 1 Block=128x128 Type=Byte, ColorInterp=Gray
  NoData Value=255
  Image Structure Metadata:
    NBITS=2
```

a.tif is a classification (0/null, 1, 2) representing two land use classes of interest.

b.tif contains small integer values (0-21) representing intensity of some phenomenon. It's on the same grid as a.tif, and was made with the same GeoTiff driver creation options.

The command essentially masks a.tif to pick class 1, then gets the value of b where the class is 1, nodata otherwise. I believe the equivalent `rio calc` command is the following:

```

rio calc "(* (where (!= (take A 1) 1) 0 1) (take B 1))" --dtype uint8 --name "A=a.tif" --name "B=b.tif" --masked --overwrite --co TILED=YES --co COMPRESS=LZW --co TFW=YES --co BLOCKXSIZE=128 --co BLOCKYSIZE=128 a.tif b.tif out.tif

```

(Incidentally, lisp syntax is a lot harder to read and write, IMO. Getting the brackets right was a challenge.)

I understand that both are using numpy masked arrays internally, so I'm not sure why rio calc fails but gdal_calc.py does not. Can anyone offer any insight?

Re: Memory error in rio calc that equivalent gdal_calc.py performs quickly without issue

Sean Gillies
 

rio-calc does indeed read the entire file. The requirement to use it on files too large to fit in memory hasn't come up before. Modifying it to work on only one (or several, using a thread or process pool) block at a time would not be complicated. I'll add an issue about this to the project tracker.

On Thu, Apr 4, 2019 at 7:40 PM <lawr@...> wrote:

Today I tried to replace a `gdal_calc.py` command with the equivalent `rio calc`, because I wanted a UInt8 output, and `gdal_calc.py` restricts me to UInt16.

 

The GDAL command is relatively quick on my 15m² raster of New Zealand, at about 2-3 minutes processing time on my machine, but the equivalent `rio calc` that I came up with seemed to be attempting to load everything into memory, and the command eventually crashed with a Memory Error after struggling for about ten minutes.

Here's the command:

```

gdal_calc.py -A a.tif -B b.tif< \
--outfile=out.tif --overwrite --calc="where(A != 1,0,1) * B" \
--format=GTiff --NoDataValue=0 --type=UInt16\
--co TILED=YES --co COMPRESS=LZW --co TFW=YES \
--co BLOCKXSIZE=128 --co BLOCKYSIZE=128
```

gdalinfo for a.tif:

```
Size is 66814, 96417
Coordinate System is:
PROJCS["NZGD_2000_New_Zealand_Transverse_Mercator",
    GEOGCS["GCS_NZGD_2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (1089805.593673222931102,6194224.566799569875002)
Pixel Size = (15.000000000000000,-15.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 1089805.594, 6194224.567) (167d27'39.65"E, 34d16' 4.53"S)
Lower Left  ( 1089805.594, 4747969.567) (166d15'35.49"E, 47d13'23.08"S)
Upper Right ( 2092015.594, 6194224.567) (178d20'32.62"E, 34d16'36.05"S)
Lower Right ( 2092015.594, 4747969.567) (179d30' 5.76"E, 47d14'12.93"S)
Center      ( 1590910.594, 5471097.067) (172d53'31.45"E, 40d54'40.25"S)
Band 1 Block=128x128 Type=Byte, ColorInterp=Gray
  NoData Value=255
  Image Structure Metadata:
    NBITS=2
```

a.tif is a classification (0/null, 1, 2) representing two land use classes of interest.

b.tif contains small integer values (0-21) representing intensity of some phenomenon. It's on the same grid as a.tif, and was made with the same GeoTiff driver creation options.

The command essentially masks a.tif to pick class 1, then gets the value of b where the class is 1, nodata otherwise. I believe the equivalent `rio calc` command is the following:

```

rio calc "(* (where (!= (take A 1) 1) 0 1) (take B 1))" --dtype uint8 --name "A=a.tif" --name "B=b.tif" --masked --overwrite --co TILED=YES --co COMPRESS=LZW --co TFW=YES --co BLOCKXSIZE=128 --co BLOCKYSIZE=128 a.tif b.tif out.tif

```

(Incidentally, lisp syntax is a lot harder to read and write, IMO. Getting the brackets right was a challenge.)

I understand that both are using numpy masked arrays internally, so I'm not sure why rio calc fails but gdal_calc.py does not. Can anyone offer any insight?



--
Sean Gillies

Is it possible to open a binary file with rasterio.open?

htaidirt@...
 

Hi all,

I've been using rasterio for few days, and now experiencing a situation where I can't find any help.

I read in the doc that when using `rasterio.open(fp,...)` it is possible to pass a binary content as `fp` (doc in the source file -> https://github.com/mapbox/rasterio/blob/c2df12979a5e07f96f108b0be8329e79fe950532/rasterio/__init__.py#L74

So instead of passing a file path or URL, I wanted to pass the content of the jp2 file to rasterio. Here is the full code for reproductability:


import boto3
import rasterio
s3 = boto3.resource('s3')

bucket = 'sentinel-s2-l1c'
key = 'tiles/31/U/DQ/2019/3/2/0/B01.jp2'

content = s3.Object(bucket, key).get(RequestPayer='requester')['Body'].read()
rio = rasterio.open(content, driver="JP2OpenJPEG").read(1)

But I'm getting an invalid path or file: b'\x00\x00\x00\x0cjP... error.

Any help on how to do so?

Thanks

Re: Is it possible to open a binary file with rasterio.open?

Sean Gillies
 

Hi,

If you pass a Python file object opened in "r" mode to rasterio.open, the contents will be read and stored in a MemoryFile. See https://github.com/mapbox/rasterio/blob/master/rasterio/__init__.py#L176. Rasterio's open function doesn't accept bytes as an argument. Note that if the JP2 file is large, on the order of your computer's RAM, your program will be at risk of running out of memory.


On Sun, Apr 7, 2019 at 7:51 AM <htaidirt@...> wrote:
Hi all,

I've been using rasterio for few days, and now experiencing a situation where I can't find any help.

I read in the doc that when using `rasterio.open(fp,...)` it is possible to pass a binary content as `fp` (doc in the source file -> https://github.com/mapbox/rasterio/blob/c2df12979a5e07f96f108b0be8329e79fe950532/rasterio/__init__.py#L74

So instead of passing a file path or URL, I wanted to pass the content of the jp2 file to rasterio. Here is the full code for reproductability:


import boto3
import rasterio
s3 = boto3.resource('s3')

bucket = 'sentinel-s2-l1c'
key = 'tiles/31/U/DQ/2019/3/2/0/B01.jp2'

content = s3.Object(bucket, key).get(RequestPayer='requester')['Body'].read()
rio = rasterio.open(content, driver="JP2OpenJPEG").read(1)

But I'm getting an invalid path or file: b'\x00\x00\x00\x0cjP... error.

Any help on how to do so?

Thanks



--
Sean Gillies

Is it possible to create a rasterio object without exporting a file?

Leah Wasser
 

Hi All,
I am running into a consistent workflow issue. i'd like to crop a numpy array. I have metadata dict and a numpy array but then want to crop it. is there a way to crop a numpy array to a new extent (using rio.mask()) using just the meta file rather than needing to export the array and reimport it to instantiate a rasterio object and create a new numpy array and associated meta?

Perhaps i'm missing something very basic!!
thank you
Leah

Re: Is it possible to create a rasterio object without exporting a file?

Alan Snow
 

Hi Leah,

Here is an example of what you want to do:
https://github.com/corteva/geocube/blob/582ea0c2f3f0ba91326de098b75412e5033bee21/geocube/xarray_extensions/rioxarray.py#L617

It uses xarray, but the concept is the same.

Hopefully this helps,
Alan

Re: Is it possible to create a rasterio object without exporting a file?

Leah Wasser
 

Hey Alan,
this is super helpful.

digging through your code. how does this line work:

cropped_ds = self._obj.where(clip_mask_xray, drop=True).astype(self._obj.dtype)

i'm familiar with numpy.where however the drop=True is not an argument i'm familiar with.
thank you so much!
Leah

Re: Is it possible to create a rasterio object without exporting a file?

Alan Snow
 

Hi Leah,

So, this project uses xarray which has many similar operations to numpy, but with slight differences.
In the case of the example, the `drop=True` is an option they provide that is described here:
http://xarray.pydata.org/en/stable/generated/xarray.DataArray.where.html

For numpy, it is not as simple if I remember correctly for a 2D array. An example of doing this is here:
https://github.com/CI-WATER/gsshapy/blob/00fd4af0fd65f1614d75a52fe950a04fb0867f4c/gsshapy/grid/grid_to_gssha.py#L616-L631

Hopefully this is helpful and gets you in the right direction.

Best,
Alan

tempfile.NamedTemporaryFile behaving as /vsimem and eating all the machine memory

vincent.sarago@...
 

While working on https://github.com/cogeotiff/rio-cogeo/pull/75 we noticed strange behaviors with `vsimem` driver (this could be a GDAL but TBH). 

1. When using `
tempfile.NamedTemporaryFile()` Rasterio uses `vsimem` driver

with tempfile.NamedTemporaryFile() as tmpfile:

   print(tmpfile)

   with rasterio.open(tmpfile, "w", **meta) as tmp_dst:

       print(tmp_dst)

<tempfile._TemporaryFileWrapper object at 0x1151bb2e8>

<open DatasetWriter name='/vsimem/d26d4650-3010-4d39-8b0c-3d947b94f1d5.' mode='w+'>

Here I was expecting Rasterio/GDAL to behave as `tempfile` was a regular file.

2. When closing
 a `vsimem`  (`MemoryFile` or `tempfile`) we observe a huge memory surge when working with big images.

code: 
https://github.com/cogeotiff/rio-cogeo/pull/75#issuecomment-482745580



Tested on Mac OS and linux with python 3.7 (gdal 2.4 and 2.3) 

Thanks 

Re: Is it possible to create a rasterio object without exporting a file?

Alan Snow
 

Hi Leah,

If you would like to try the rioxarray as an xarray extension, it is now available on pypi.

An example of what you would like to do is here:
https://corteva.github.io/rioxarray/html/examples/clip_geom.html

Best,
Alan

Re: Is it possible to create a rasterio object without exporting a file?

Leah Wasser
 

Thank you, Alan!
that example is very helpful. i'll dig into this a bit deeper. xarray seems like a powerful tool for raster processing.
i've wanted to plan with it a bit more regardless.
Leah

Re: tempfile.NamedTemporaryFile behaving as /vsimem and eating all the machine memory

Sean Gillies
 

Hi Vincent,

This is expected (if not well-documented) behavior. tempfile.NamedTemporaryFile() returns an open Python file object, not a filename. GDAL can't use a Python file object, so in that case rasterio.open reads all the bytes from the file object, copies them to the vsimem filesystem, and opens that vsimem file.

I think what you want do do is pass the name of the temp file object to GDAL. Like this:

with tempfile.NamedTemporaryFile() as temp:
    with rasterio.open(temp.name) as dataset:
        print(dataset)

No copy in the vsimem filesystem will be made.

On Tue, Apr 16, 2019 at 6:55 AM <vincent.sarago@...> wrote:
While working on https://github.com/cogeotiff/rio-cogeo/pull/75 we noticed strange behaviors with `vsimem` driver (this could be a GDAL but TBH). 

1. When using `
tempfile.NamedTemporaryFile()` Rasterio uses `vsimem` driver

with tempfile.NamedTemporaryFile() as tmpfile:

   print(tmpfile)

   with rasterio.open(tmpfile, "w", **meta) as tmp_dst:

       print(tmp_dst)

<tempfile._TemporaryFileWrapper object at 0x1151bb2e8>

<open DatasetWriter name='/vsimem/d26d4650-3010-4d39-8b0c-3d947b94f1d5.' mode='w+'>

Here I was expecting Rasterio/GDAL to behave as `tempfile` was a regular file.

2. When closing
 a `vsimem`  (`MemoryFile` or `tempfile`) we observe a huge memory surge when working with big images.

code: 
https://github.com/cogeotiff/rio-cogeo/pull/75#issuecomment-482745580



Tested on Mac OS and linux with python 3.7 (gdal 2.4 and 2.3) 

Thanks 



--
Sean Gillies

Re: tempfile.NamedTemporaryFile behaving as /vsimem and eating all the machine memory

vincent.sarago@...
 

Thanks Sean this is really helpful and love the `temp.name` solution. 

About the second point, do you have any idea why `/vsimem` driver need so much memory when exiting/closing ? Should I raise this to the gdal list?