On mercredi 29 avril 2020 08:40:07 CEST Guillaume Lostis wrote:

> Like Ciaran I'm not very familiar with GDAL's python bindings, but I've

> tried running a little snippet which I believe does the same thing as what

> was done with `rasterio`: it writes a JP2 to a `vsimem` in-memory file, and

> then writes the bytes to a file on disk.

>

> ```python

> from osgeo import gdal

>

> mem, out = "/vsimem/B01.jp2", "new_B01.jp2"

> gdal.Translate(mem, "B01.jp2")

>

> # Taken from

> https://lists.osgeo.org/pipermail/gdal-dev/2016-August/045030.html f =

> gdal.VSIFOpenL(mem, "rb")

> gdal.VSIFSeekL(f, 0, 2) # seek to end

> size = gdal.VSIFTellL(f)

> gdal.VSIFSeekL(f, 0, 0) # seek to beginning

> data = gdal.VSIFReadL(1, size, f)

> gdal.Unlink(mem)

>

> with open(out, "wb") as ff:

> ff.write(data)

> ```

>

> When I look at the `gdalinfo` of `new_B01.jp2`, it has a CRS and a

> transform, so I am not able to reproduce the behavior we have with

> `rasterio`.

>

> (Side note: Interestingly enough, the pixel values of `new_B01.jp2` have

> slightly changed with respect to the original `B01.jp2` file, so there is

> some data loss somewhere in the process. But maybe that is expected of the

> JP2 format and could be avoided by passing extra arguments to

> `gdal.Translate`?

By default, the JP2OpenJPEG driver uses a lossy compression (my personal take on JPEG2000 is that it has very limited interested when used in its lossless profile. Better use more conventional algorithms that are way faster).

See https://gdal.org/drivers/raster/jp2openjpeg.html#lossless-compression for options to pass to get a lossless JPEG2000 as output.

So all in all that should be something like

gdal.Translate(mem, "B01.jp2",

format = "JP2OpenJPEG", # in case you have several JP2K drivers available

options = ["REVERSIBLE=YES", "QUALITY=100"])

--

Spatialys - Geospatial professional services

http://www.spatialys.com

Like Ciaran I'm not very familiar with GDAL's python bindings, but I've tried running a little snippet which I believe does the same thing as what was done with rasterio: it writes a JP2 to a vsimem in-memory file, and then writes the bytes to a file on disk.

from osgeo import gdal

mem, out = "/vsimem/B01.jp2", "new_B01.jp2"
gdal.Translate(mem, "B01.jp2")

# Taken from https://lists.osgeo.org/pipermail/gdal-dev/2016-August/045030.html
f = gdal.VSIFOpenL(mem, "rb")
gdal.VSIFSeekL(f, 0, 2)  # seek to end
size = gdal.VSIFTellL(f)
gdal.VSIFSeekL(f, 0, 0)  # seek to beginning
data = gdal.VSIFReadL(1, size, f)
gdal.Unlink(mem)

with open(out, "wb") as ff:
    ff.write(data)

When I look at the gdalinfo of new_B01.jp2, it has a CRS and a transform, so I am not able to reproduce the behavior we have with rasterio.

(Side note: Interestingly enough, the pixel values of new_B01.jp2 have slightly changed with respect to the original B01.jp2 file, so there is some data loss somewhere in the process. But maybe that is expected of the JP2 format and could be avoided by passing extra arguments to gdal.Translate? I mainly have experience with writing GeoTIFFs, writing JP2s is new to me)

Guillaume

Would you say it's worth raising an issue on GDAL then?

I'm also not familiar with the Python GDAL bindings, Rasterio is my window in that haha.

A workaround could be to dismantle the geometry collections into separate objects. It would work but it's not ideal as the collections are there to group various styles of object together.

Does anyone have experience with this issue, or can suggest another solution? My code etc below.

[{'type': 'LineString',
  'coordinates': [(144.937247997, -26.3759745069999, 0.0),

{'type': 'GeometryCollection',
  'geometries': [{'type': 'LineString',
    'coordinates': [(145.619783254, -38.5831492689999, 0.0),


My code: 
with fiona.open('./Methane gas detection docs/NationalOnshoreGasPipelines.geojson') as fi:
    shapes = [each['geometry']for each in fi]

    # isolate a pipeline area
    pipes = rasterio.features.geometry_window(win, shapes, pad_x = 0.001, pad_y=0.001)



Error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-183-b9dff7a7db24> in <module>
     10 
     11     # isolate a pipeline area
---> 12     pipes = rasterio.features.geometry_window(win, shapes, pad_x = 0.001, pad_y=0.001)
     13 
     14     # print(shapes)

~\Anaconda3\envs\sat\lib\site-packages\rasterio\features.py in geometry_window(dataset, shapes, pad_x, pad_y, north_up, rotated, pixel_precision)
    419 
    420     if not rotated:
--> 421         all_bounds = [bounds(shape, north_up=north_up) for shape in shapes]
    422         lefts, bottoms, rights, tops = zip(*all_bounds)
    423 

~\Anaconda3\envs\sat\lib\site-packages\rasterio\features.py in <listcomp>(.0)
    419 
    420     if not rotated:
--> 421         all_bounds = [bounds(shape, north_up=north_up) for shape in shapes]
    422         lefts, bottoms, rights, tops = zip(*all_bounds)
    423 

~\Anaconda3\envs\sat\lib\site-packages\rasterio\features.py in bounds(geometry, north_up, transform)
    374 
    375     geom = geometry.get('geometry') or geometry
--> 376     return _bounds(geom, north_up=north_up, transform=transform)
    377 
    378 

rasterio\_features.pyx in rasterio._features._bounds()

KeyError: 'coordinates'

Hello, I have an application that calls rasterio to open rasters in its main process. This application allows me to parallelize the process, so it will call the function where rasterio open is many times in parallel.
If I open and close the image every time by doing

with rasterio.open(self.file) as src:
    raster = src.read(window=window)

The process goes fine, but it has to open and close at every step of the raster (losing a bit of time). Then I tried to let the image opened at the beginning of the process but it is not stable. If I don't have too many workers, the application might work well, but if I increase the number of workers it will crash with the following error:

RasterioIOError: Caught RasterioIOError in DataLoader worker process 0.
Original Traceback (most recent call last):
File "rasterio/_io.pyx", line 707, in rasterio._io.DatasetReaderBase._read
File "rasterio/shim_rasterioex.pxi", line 133, in rasterio._shim.io_multi_band
File "rasterio/_err.pyx", line 182, in rasterio._err.exc_wrap_int
rasterio._err.CPLE_AppDefinedError: ./sentinel2_tiled.tif, band 1: IReadBlock failed at X offset 58, Y offset 92: TIFFReadEncodedTile() failed.

What do you recommend me to do?

Thank you.

As with Guillaume's answer, it's a Sentinel 2 band and has no side-car file, when I write straight to a file, the georeferencing is maintained.

I tried this with GDAL 2.4.2 and GDAL 3.0 too

If anyone can point me to where I might look further to diagnose this, I can create an issue with some useful info, whether that's in Rasterio/GDAL :)

Revised: Don't require the seek call:

s3_r = boto3.resource('s3')
profile = source_image.profile.copy()

with rasterio.io.MemoryFile() as stack:
    with stack.open(**profile) as stack_out:
        stack_out.write(source_image.read(1), 1)
    s3_r.Object(<my-bucket>, <my-key>).upload_fileobj(stack)

For reference, doing:

s3_r = boto3.resource('s3')
profile = source_image.profile.copy()

with rasterio.io.MemoryFile() as stack:
    profile['driver'] = 'GTiff'
    with stack.open(**profile) as stack_out:
        stack_out.write(source_image.read(1), 1)
    s3_r.Object(<my-bucket>, <my-key>).upload_fileobj(stack)

results in a correctly geocoded GeoTiff. Appears it is an issue with .jp2s and Profiles/MemoryFiles.

On further inspection, when I switched to saving my outputs as GeoTIFF, the resultant images are both geo-coded fine (Disk vs. In Memory to S3).

I'm guessing something is going wrong when .jp2 is involved...

Any suggestions on next steps?

.jp2 is both preferable (because of file sizes) and a requirement for the work I'm currently doing...

s3_r = boto3.resource('s3')
profile = source_image.profile.copy()

with rasterio.io.MemoryFile() as stack:
    with stack.open(**profile) as stack_out:
        stack_out.write(source_image.read(1), 1)
    stack.seek(0)
    s3_r.Object(<my-bucket>, <my-key>).upload_fileobj(stack)

If I download that image and open it in QGIS, the image has no geo-coding at all.

If I process the same image, but write to disk like:

profile = source_image.profile.copy()

with rasterio.open('thing.jp2', 'w', **profile) as stack_out:
    stack_out.write(source_image.read(1), 1))

I get a geo-coded image and it appears where I'd expect.

If I compare file sizes:
MemoryFile(): 49,211KB
Disk: 49,214KB

I can only assume I'm losing any profile information when using a MemoryFile, can someone direct me to what I'm doing wrong?

Thanks!