Date   

Re: Clarification on usage with QGIS

Ratcliff, Christina (A&F, Waite Campus)
 

Ari, I've been running the Geopandas (pandas, shapely fiona), RasterIO configuration since early 2019, and as long as you have access to those gdal backwards compatibility dll's you should be fine.

With libraries depending on GDAL, where possible use those available in osgeo4w. Osgeo4w will then properly install and get them talking to both GDAL and QGIS.

Good luck.


Re: How to specify the name of the input files with the geographical metadata?

juanjo.gomeznavarro@...
 

Thanks. It's a bit strange that this is not a native feature of GDAL, as it seems pretty natural to work on bunches of files that share the same geographical information. But anyway, I have finally opted for using programmatically generated VRT files.


Re: Georeference and warp a drone image based on sensor orientation and location using transform.from_gcps()

Patrick Gray
 

Hi Sean and thanks for the reply! I've changed my code to use rasterio.warp.reproject, but the final geotiff is still a rectangle rather than a trapezoid and not actually warping to the GCPs. It is roughly where it should be but not even in the middle of the GCPs which are at each corner of the image. Should I be using that rasterio.transform.from_gcps() function to get the transform for the reproject() function? I must be missing something... Here is my code:

import rasterio
from rasterio.warp import reproject
from rasterio.enums import Resampling
 
# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():
    
    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile
        
        print(profile)
        # create the transform from the GCPs
        tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
        print(np.array(src.read()))
        
        src_crs = "EPSG:4326"  # This is the crs of the GCPs
 
        # Destination: a 1024 x 1024 dataset in WGS84
        dst_raster = np.zeros((1024, 1024), dtype=np.uint16)
        dst_crs = "EPSG:4326"
 
        tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
 
        new_raster, dst_transform = reproject(
            source=np.array(src.read()),
            destination=dst_raster,
            dst_transform=tsfm,
            gcps=[gcp1,gcp2,gcp3,gcp4],
            src_crs=src_crs,
            dst_crs=dst_crs,
            resampling=Resampling.nearest,
            #**kwargs
        )
 
        profile.update(
            dtype=rasterio.uint16,
            transform = dst_transform,
            crs=dst_crs,
            width=dst_raster.shape[0], # TODO unsure if this is correct order
            height=dst_raster.shape[1]
        )
 
        print(profile)
 
        with rasterio.open('example_uint16_tsfm.tif', 'w', **profile) as dst:
            dst.write(src.read(1).astype(rasterio.uint16), 1)


And here is what the image looks like in WGS84:


Re: How to specify the name of the input files with the geographical metadata?

Sean Gillies
 

Hi,

On Tue, Feb 23, 2021 at 9:28 AM <juanjo.gomeznavarro@...> wrote:
Good morning,

I'm working with PNG files, and the metadata is stored in two separate files: .pgw and .png.aux.xml. Currently I'm forced to use the same name for the three files, i.e. the image and the two auxiliary files with the geographical information. This is cumbersome and inconvenient, because I have a set of images that share the same georeferentiation and currently I'm forced to clone the .pgw and .png.aux.xml files in a somewhat silly way for introducing each image. Is there a way to tell rasterio.open that the two files with the metadata are named different from the PNG image?

Thanks a lot!

The naming of the auxiliary files is a GDAL convention and cannot, to my knowledge, be changed. Neither does rasterio allow a user to open an image in read only mode and add or replace the spatial referencing and other metadata. I agree that would be handy in your case. 

Workarounds include:

Opening your PNGs in "r+" mode and patching them, like

with rasterio.open("example.png", "r+") as dataset:
    dataset.crs = rasterio.crs.CRS.from_epsg(3857)
    dataset.transform = rasterio.transform.Affine.translation(ulx, uly) * rasterio.transform.Affine.scale(dx, -dy)

Note that this would produce .pgw and .png.aux.xml files on your system.

Or you could make a generic VRT XML wrapper template and substitute PNG paths into it. See https://gdal.org/drivers/raster/vrt.html. The XML string can then be opened by rasterio.

with rasterio.open("<VRTDataset>...</VRTDataset>") as dataset:
    ...

That's not a well known feature of GDAL/rasterio but I use it quite a bit.

--
Sean Gillies


How to specify the name of the input files with the geographical metadata?

juanjo.gomeznavarro@...
 

Good morning,

I'm working with PNG files, and the metadata is stored in two separate files: .pgw and .png.aux.xml. Currently I'm forced to use the same name for the three files, i.e. the image and the two auxiliary files with the geographical information. This is cumbersome and inconvenient, because I have a set of images that share the same georeferentiation and currently I'm forced to clone the .pgw and .png.aux.xml files in a somewhat silly way for introducing each image. Is there a way to tell rasterio.open that the two files with the metadata are named different from the PNG image?

Thanks a lot!


Re: Georeference and warp a drone image based on sensor orientation and location using transform.from_gcps()

Sean Gillies
 

Hi,

On Mon, Jan 11, 2021 at 11:10 AM <pgrayobx@...> wrote:
I have been flying some drone surveys over the ocean and need to properly project and georeference and warp my images. I have all the information I think I need: lat, lon, altitude, yaw, pitch, and roll along with sensor optics params. I have to imagine this can be done with some existing package and I had thought some rasterio functions were suitable, but I can't find any. So I've been using the cameratransform package to get the GPS positions of the corners of my image to feed into rasterio:
 
import cameratransform as ct
 
# camera parameters
 
cam = ct.Camera(ct.RectilinearProjection(focallength_mm=f,
                                         sensor=sensor_size,
                                         image=image_size),
               ct.SpatialOrientation(elevation_m=img.altitude,
                                     tilt_deg=pitch+sensor_offset,
                                     roll_deg=roll,
                                    heading_deg=yaw))
 
# gps pts are lat lon
cam.setGPSpos(img.latitude, img.longitude, img.altitude)
 
# these are the coordinates of the image corners
coords = np.array([cam.gpsFromImage([0 , 0]), \
    cam.gpsFromImage([image_size[0]-1 , 0]), \
    cam.gpsFromImage([image_size[0]-1, image_size[1]-1]), \
    cam.gpsFromImage([0 , image_size[1]-1])])

From there I'm not certain what to do. I've tried basically saying these corners are GCPs and trying to warp the image in rasterio:
 
import rasterio
 
gcp1 = rasterio.control.GroundControlPoint(row=0, col=0, x=coords[0,1], y=coords[0,0], z=coords[0,2], id=None, info=None)
gcp2 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=0, x=coords[1,1], y=coords[1,0], z=coords[1,2], id=None, info=None)
gcp3 = rasterio.control.GroundControlPoint(row=image_size[0]-1, col=image_size[1]-1, x=coords[2,1], y=coords[2,0], z=coords[2,2], id=None, info=None)
gcp4 = rasterio.control.GroundControlPoint(row=0, col=image_size[1]-1, x=coords[3,1], y=coords[3,0], z=coords[3,2], id=None, info=None)
 
# Register GDAL format drivers and configuration options with a
# context manager.
with rasterio.Env():
    
    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile
                    
    # create rasterio transform
    tsfm = rasterio.transform.from_gcps([gcp1,gcp2,gcp3,gcp4])
    
    # I also tried this function but to no avail
    #tsfm = rasterio.warp.calculate_default_transform(rasterio.crs.CRS({"init": "epsg:4326"}), rasterio.crs.CRS({"init": "epsg:4326"}), img.size()[0], img.size()[1], gcps=[gcp1,gcp2,gcp3,gcp4])
 
    crs = rasterio.crs.CRS({"init": "epsg:4326"})
 
    profile.update(
        dtype=rasterio.uint8,
        transform = tsfm,
        crs=crs)
        
    with rasterio.open('example.tif', 'w', **profile) as dst:
        dst.write(src.read().astype(rasterio.uint8), 1)

This does produce an image but it is not properly warped. It is just in the correct location. Since my sensor is a 40 deg off nadir the warping should be reasonably significant, when I warp it with cameratransform's cam.getTopViewOfImage() function I get this:

example image

So I'm curious if I'm missing some rasterio function or ability, isn't what the GCPs should be doing in the transform.from_gcps() function?

First of all, I apologize for not responding sooner.

You must use the rasterio.warp.reproject function to warp. Warping is never implicit in the dataset read/write methods, these are very low-level and always operate strictly in the dataset's row-column space.


--
Sean Gillies


Re: Clarification on usage with QGIS

Ari Meyer
 

That's exactly what I hoped to hear, Christina -- thank you so much!

One question: are there any other libraries commonly used in the QGIS space that depend upon a particular version of GDAL?  I'm hoping I won't be facing any version conflicts with this setup.
Ari


Re: Clarification on usage with QGIS

Ratcliff, Christina (A&F, Waite Campus)
 

I've successfully developed the PAT plugin for QGIS using rasterio,  geopandas (uses fiona, shapely). All of these libraries are available in the osgeo4w installer through the advanced install. You will also have to ensure the correct GDAL dll is also installed for the version of rasterio you are using.



My method was to create a windows bat file containing the relevant installation which then can be run as administrator to install the libraries. This will use a command line install of rasterio, geopandas and GDAL dll via OSGeo4W.

Have a look at the following osgeo4w and QGIS  tickets for answers to some of the issues I came across.

I hope this helps.
Christina


Re: Open Santinel 2 archives with Rasterio

lucavlad50@...
 

That's it, thank you very much!


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/.

101 - 120 of 828