1. main@rasterio.groups.io
   Date   
Re: Open Santinel 2 archives with Rasterio

Luke
 

GDALs Sentinel2 driver exposes the data as subdatasets, with (sort of) one for each resolution - https://gdal.org/drivers/raster/sentinel2.html#level-2a

To access the subdataset that contains the 10m bands, you could use something like the following:

import rasterio

with rasterio.open('S2A_MSIL2A_20210214T022811_N0214_R046_T51TXL_20210214T044441.zip') as s2a:
    subdatasets = s2a.subdatasets

with rasterio.open(subdatasets[0]) as b10m:
    print(b10m.profile)
    
{'driver': 'SENTINEL2', 'dtype': 'uint16', 'nodata': None, 'width': 10980, 'height': 10980, 'count': 4, 'crs': CRS.from_epsg(32651), 'transform': Affine(10.0, 0.0, 600000.0,
       0.0, -10.0, 5100000.0), 'blockxsize': 128, 'blockysize': 128, 'tiled': True, 'compress': 'jpeg2000'}
Re: Open Santinel 2 archives with Rasterio

lucavlad50@...
 

Hello,
Thanks for the reply. I saw that but there is no example of a zip file being opened. I saw this: 
on: https://gdal.org/drivers/raster/sentinel2.html
and was wondering if it is possible to read directly from the zip file without unzipping it. 
Something like: 
with rio.open("satelite_datasets/S2A_MSIL2A_20201020T092031_N0214_R093_T34TGS_20201020T120739.zip") as f:
    print(f.read(4))
where 4 is band B04 in  R10mm.

Sorry if the question seems stupid, I am new to this. Thank you!

Re: Open Santinel 2 archives with Rasterio

Sean Gillies
 

Hi,

Please see https://rasterio.readthedocs.io/en/latest/quickstart.html#reading-raster-data for a quick tour of data access features.

On Fri, Feb 12, 2021 at 8:01 AM <lucavlad50@...> wrote:
Hello,
I saw that rasterio support reading Santinel 2 archives directly, but I can not find any example. How should I access the bands inside the archive after I read it. 
For example I have:

with rio.open("satelite_datasets/S2A_MSIL2A_20201020T092031_N0214_R093_T34TGS_20201020T120739.zip") as f:
    print(...)
What should I put instead of "..." to access for example R10 band 4? 

Thank you in advance!
_._,_._,

--
Sean Gillies
Open Santinel 2 archives with Rasterio

lucavlad50@...
 

Hello,
I saw that rasterio support reading Santinel 2 archives directly, but I can not find any example. How should I access the bands inside the archive after I read it. 
For example I have:

with rio.open("satelite_datasets/S2A_MSIL2A_20201020T092031_N0214_R093_T34TGS_20201020T120739.zip") as f:
    print(...)
What should I put instead of "..." to access for example R10 band 4? 

Thank you in advance!
Re: Reading and plotting rotated rasters

Sean Gillies
 

Hi Denis,

On Fri, Feb 5, 2021 at 4:16 PM Denis Rykov <rykovd@...> wrote:
Is it expected that during reading process rasterio doesn't take into account rotation angle for rotated rasters (transform has non-zero rotation) ? In order to make it works properly I had to do gdalwarp.

Same question for the "show" method. It just draws array as it is without considering the rotation angle.


--
Sean Gillies
reproject: partially covered raster cells reset fully to nodata

Fabian Hofmann
 

Dear community,


I'm trying to tackle the following problem: Raster A has a high resolution and a small extent. I should be projected to Raster B with small resolution and large extent. Cells from Raster A which partially cover cells in Raster B are getting ignored in the reprojection. Only cells in Raster B which are fully (!) covered by Raster A are getting non-nodata entries. I am not sure if this is intended.


Reproduce the behavior with


```

import numpy as np
import rasterio as rio
from rasterio.warp import reproject, transform_bounds


src_crs = 3035
src_top = 3200000.0
src_left = 4000000.0
src_res = 100
src_shape = (100, 100)
src_transform = rio.Affine(src_res, 0.0, src_left, 0.0, -src_res, src_top)
src_right, src_bottom = src_transform * src_shape
src = np.ones(src_shape)


dst_crs = 4326
dst_top = 52
dst_left = 5
dst_res = 0.25
dst_shape = (10, 10)
dst_transform = rio.Affine(dst_res, 0, dst_left, 0, -dst_res, dst_top)
dst_right, dst_bottom = dst_transform * dst_shape
dst = np.zeros(dst_shape)


# check destination bounds comprise source bounds
covered = prj_left, prj_bottom, prj_right, prj_top = (
    transform_bounds(src_crs, dst_crs, src_left, src_bottom, src_right, src_top))
assert prj_left > dst_left and prj_bottom > dst_bottom
assert prj_right < dst_right and prj_top < dst_top


reproject(src, dst,
          src_transform=src_transform, src_crs=src_crs,
          dst_transform=dst_transform, dst_crs=dst_crs,
          resampling=rio.warp.Resampling.average,
          src_nodata=0, dst_nodata=0, init_dest_nodata=True)

assert dst.max() > 0

```

As soon as you change the target resolution to dst_res=0.1, the reprojected values contain non-zero values. However cells which should have some values, because Raster A partially covers it, get only nodata (in this case 0).


What I'm doing to avoid the issue is to pad around the src values. But this needs an extra calculation to get the needed pad_width in order to cover the whole cells in Raster B. (The operation is done many times). Is there a better way? Or is this even a bug?


Hope the explanation was clear, thanks in advance!


Best,

Fabian Hofmann


--
Fabian Hofmann
Frankfurt Institute for Adavanced Studies (FIAS)
FIAS Renewable Energy Systems & Network Analysis (FRESNA)

Ruth-Moufang-Straße 1, 60438 Frankfurt am Main
T +49 69 798 47605 | hofmann@fias.uni-frankfurt.de
Reading and plotting rotated rasters

Denis Rykov
 

Is it expected that during reading process rasterio doesn't take into account rotation angle for rotated rasters (transform has non-zero rotation) ? In order to make it works properly I had to do gdalwarp.

Same question for the "show" method. It just draws array as it is without considering the rotation angle. 
Re: HDF4 format support

Alan Snow
 

HDF4 files are currently not supported out of the box with rasterio wheels (https://github.com/rasterio/rasterio-wheels/issues/36).

To support HDF4 files, you will need to build GDAL yourself with HDF4 support enabled and install rasterio with"--no-binary rasterio" or conda with the conda-forge channel to install rasterio (I believe conda has HDF4 support in their distribution).

Your specific instance is something I cannot speak to as it seems to be a bit of an anomaly.
HDF4 format support

bsmith@...
 


```python
import rasterio as rio

from osgeo import gdal
```

I cannot open HDF4 files without including the osgeo gdal library.
Is there a reason that osgeo isn't a dependency in order to open HDF4 files?
Re: Clarification on usage with QGIS

Ari Meyer
 

Looking through the older messages, it seems like integration with QGIS, in particular, is a common enough use case that it really should be discussed separately in the docs.  I will be making a request to the QGIS team that Rasterio be included in the core distribution.  A couple questions for the group:

1)  Does anyone have experience writing QGIS plugins using Rasterio?  If so, did you run into any problems, particularly those related to issue I pointed to above?

2)  Is there a technical reason why Shapely was included in QGIS but Rasterio was not?

Thanks,
Ari
Re: Clarification on usage with QGIS

Ari Meyer
 

FYI, Sean Gillies responded on https://github.com/mapbox/rasterio/issues/2096 :

I'm certain that until rasterio is included with QGIS, shared library compatibility ("DLL Hell") will be a problem for your plugin. For a compiled program like QGIS, it's fairly simple to make sure that only one instance of a library is loaded, this is the job of the linker. In a Python program, there's no built in way to prevent conflicts when different extension modules cause different instances of the same shared library to be loaded. If I were writing QGIS plugins, I would try to stick to using https://qgis.org/pyqgis/master/.
Clarification on usage with QGIS

Ari Meyer
 

Note: I originally submitted this as a github issue, as I'd like to see an update to the documentation regarding using Rasterio with QGIS (and perhaps for integrating with any apps that use GDAL).

We're looking to simplify our code for a QGIS plugin to use Rasterio/Shapely in place of GDAL/OGR. Shapely is included in the standard QGIS distribution, but Rasterio is not (not sure why). In any case, I found on
https://rasterio.readthedocs.io/en/latest/topics/switch.html :
GDAL’s bindings (gdal for the rest of this document) and Rasterio are not entirely compatible and should not, without a great deal of care, be imported and used in a single Python program. The reason is that the dynamic library they each load (these are C extension modules, remember), libgdal.so on Linux, gdal.dll on Windows, has a number of global objects and the two modules take different approaches to managing these objects.

Does this imply that we may have problems if we try to port the GDAL Python code of our QGIS plugin to Rasterio that we may run into problems, as QGIS is built on GDAL? Can you add a section to rasterio.readthedocs.io clarifying interoperability with QGIS?
GDAL Error 4 triggers when working with geo-referenced jpeg

ashnair0007@...
 
Edited

I used gdal_translate to convert a geotiff to a jpg with jpg.aux.xml file.  The following is part of the script I used to get and write the tiles to disk

import rasterio as rio

im = rio.open('./rgb.tif')
meta = im.meta.copy()

for window, transform in get_tiles(w,h,im):
   tile_name = f"{window.col_off}_{window.row_off}"
   meta["transform"] = transform
   meta["width"], meta["height"] = window.width, window.height
   outpath = "./tiles/" + tile_name
   with rio.open(outpath, "w", **meta) as outds:
       outds.write(self.im.read(window=window))

Now when the input is the original geotiff, writing the tiles to disk works perfectly. However, if my input is the jpeg that I created via gdal_translate, writing to tiles still works however it keeps printing out 
ERROR 4: tiles/0_0.jpg: No such file or directory for each tile. So it'll go something like this:

ERROR 4: tiles/0_0.jpg: No such file or directory
ERROR 4: tiles/0_1000.jpg: No such file or directory
ERROR 4: tiles/0_2000.jpg: No such file or directory
.........................

and so on despite successfully writing the jpeg tiles. I'm confused as to why this error gets printed out and would appreciate any clarification regarding this behaviour.
Re: TypeError: open() takes at most 8 arguments (14 given) when writing

ashnair0007@...
 

Yes. That was exactly it. My bad. Thanks.
Re: TypeError: open() takes at most 8 arguments (14 given) when writing

Alan Snow
 

I think the issue is that you need to use `rasterio.open`.
TypeError: open() takes at most 8 arguments (14 given) when writing

ashnair0007@...
 

Hi guys,

I'm trying to write to an array to a tiff as follows:

profile.update(driver="GTiff", count=1)
with open('res.tiff', 'w', **profile) as dst:
dst.write(thresh.astype(rio.uint8), 1)

However, I'm encountering the following error:
TypeError: open() takes at most 8 arguments (14 given)

This is my profile:

{'driver': 'JPEG', 'dtype': 'uint8', 'nodata': None, 'width': 1000, 'height': 1000, 'count': 3, 'crs': CRS.from_wkt('PROJCS["unknown",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",-90],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1],AXIS["Easting",NORTH],AXIS["Northing",NORTH]]'), 'transform': Affine(25.610955512797485, 0.0, -2357528.4841969036,
       0.0, -25.610955512797485, 1593925.6448274788), 'tiled': False, 'compress': 'jpeg', 'interleave': 'pixel', 'photometric': 'ycbcr'}

I'm unsure what fields must be removed for open() to work.
16 bit to 8 bit conversion using rasterio

ashnair0007@...
 

How can I convert a 16 bit tiff to an 8 bit tiff using rasterio. I can do this via gdal_translate so I was wondering if there's a rasterio equivalent?
Re: Merge bug?

Marco
 


Yes, all rasters have nodata set to match the correct nodata value. QGIS correctly mask the nodata zone(transparent).
Re: Merge bug?

Alan Snow
 

One thing to check would be if the nodata values are properly set on the rasters. If they are not, then merging could cause the nodata values to fill in those areas instead of the correct values.
Merge bug?

Marco
 

I want to make a mosaic of satellite images.
I have 7 images clipped each one to fill my BBOX exactly. Loading that images in QGIS i can see that is totally fulfilled and without gaps.

Now, reading these images with rasterio and than merging(files_to_join it's a list of opened images with rasterio.open)

mosaic, out_trans = merge(files_to_join, precision=50)

and saving it with out_trans and sizes updated, I get:

that has gaps and missing pieces.

Why this is happening? (Same CRS and same ZONE)
41 - 60 of 759