rasterio src_transform ignored by WarpedVRT?

Expected behavior and actual behavior.

I use WarpedVRTs to reproject and crop Sentinel-1 GRD images on AWS. These images don't have a geotransform, but they have GCPs, and thanks to #1773 we can estimate a transform from the GCPs.

I use something like this (inspired by what @vincentsarago does in rio-tiler):

with rasterio.open("s1.tiff") as src:
    with WarpedVRT(
        src,
        src_crs=src.gcps[1],
        src_transform=from_gcps(src.gcps[0]),
        crs=32630,
    ) as vrt:
        ...

This works well, but I've noticed that when I change the src_transform that I pass to WarpedVRT, by doing for example:

src_transform=Affine.translation(0.1, 0) * from_gcps(src.gcps[0])

it doesn't seem to change the content of my cropped image. I can even remove the src_transform argument altogether, and the produced image is still the same!

This leads me to think that the src_transform keyword argument of WarpedVRT has no effect on the VRT dataset returned.

By reading the code of WarpedVRTReaderBase in rasterio/_warp.pyx, I see that src_transform is only used in three lines: https://github.com/mapbox/rasterio/blob/6c9816cb10876cd25dd5c802a52490628adc0ec1/rasterio/_warp.pyx#L824-L828

As far as I can tell, the src_transform is only used in the case where dst_height, dst_width and dst_transform are not None, through the src_gt variable here: https://github.com/mapbox/rasterio/blob/6c9816cb10876cd25dd5c802a52490628adc0ec1/rasterio/_warp.pyx#L909-L911

In all other cases, the variable is not used.

Steps to reproduce the problem.

import rasterio
from rasterio.session import AWSSession
from rasterio.transform import Affine, from_gcps
from rasterio.vrt import WarpedVRT

path = "s3://sentinel-s1-l1c/GRD/2020/3/17/IW/DV/S1B_IW_GRDH_1SDV_20200317T173959_20200317T174024_020735_0274FA_1C70/measurement/iw-vv.tiff"
session = AWSSession(requester_pays=True)
with rasterio.Env(session):
    with rasterio.open(path) as src:
        transform = from_gcps(src.gcps[0])
        print(f"Transform passed to WarpedVRT\n{transform}\n")
        with WarpedVRT(src, src_crs=4326, src_transform=transform, crs=4326) as warped:
            print(f"warped.src_transform\n{warped.src_transform}\n")
            print(f"warped.transform\n{warped.transform}\n")

        print()
        new_transform = Affine.translation(0.1, 0) * transform
        print(f"Transform passed to WarpedVRT\n{new_transform}\n")
        with WarpedVRT(src, src_crs=4326, src_transform=new_transform, crs=4326) as warped:
            print(f"warped.src_transform\n{warped.src_transform}\n")
            print(f"warped.transform\n{warped.transform}\n")

Transform passed to WarpedVRT
| 0.00,-0.00, 0.84|
| 0.00, 0.00, 48.20|
| 0.00, 0.00, 1.00|

warped.src_transform
| 0.00,-0.00, 0.84|
| 0.00, 0.00, 48.20|
| 0.00, 0.00, 1.00|

warped.transform
| 0.00, 0.00, 0.42|
| 0.00,-0.00, 50.09|
| 0.00, 0.00, 1.00|


Transform passed to WarpedVRT
| 0.00,-0.00, 0.94|
| 0.00, 0.00, 48.20|
| 0.00, 0.00, 1.00|

warped.src_transform
| 0.00,-0.00, 0.94|
| 0.00, 0.00, 48.20|
| 0.00, 0.00, 1.00|

warped.transform
| 0.00, 0.00, 0.42|
| 0.00,-0.00, 50.09|
| 0.00, 0.00, 1.00|

As you can see, regardless of the src_transform we pass, the .transform attribute of the warped VRT doesn't change, and I don't really know where it comes from (my gut feeling is that it is somehow estimated from the GPCs by GDAL, but in that case it's surprising that it's not equal to my from_gcps() transform).

Operating system

Ubuntu

Rasterio version and provenance

rasterio built from master (6c9816cb10876cd25dd5c802a52490628adc0ec1)

Mar 18 '20 17:03 glostis

@glostis indeed, this is interesting you can also totally remove any transform and you'll get the same result

path = "s3://sentinel-s1-l1c/GRD/2020/3/17/IW/DV/S1B_IW_GRDH_1SDV_20200317T173959_20200317T174024_020735_0274FA_1C70/measurement/iw-vv.tiff"  
session = AWSSession(requester_pays=True)  
with rasterio.Env(session):  
    with rasterio.open(path) as src: 
        print("source transform") 
        print(src.transform) 
        print() 
        print("transform from gcps") 
        print(from_gcps(src.gcps[0])) 
        print() 
        with WarpedVRT(src, src_crs=4326) as warped:  
            print(f"warped.src_transform\n{warped.src_transform}\n")  
            print(f"warped.transform\n{warped.transform}\n")

source transform
| 1.00, 0.00, 0.00|
| 0.00, 1.00, 0.00|
| 0.00, 0.00, 1.00|

transform from gcps
| 0.00,-0.00, 0.84|
| 0.00, 0.00, 48.20|
| 0.00, 0.00, 1.00|

warped.src_transform
| 1.00, 0.00, 0.00|
| 0.00, 1.00, 0.00|
| 0.00, 0.00, 1.00|

warped.transform
| 0.00, 0.00, 0.42|
| 0.00,-0.00, 50.09|
| 0.00, 0.00, 1.00|

As you pointed because we don't pass width and height to WarpedVRT, rasterio will go directly to https://github.com/mapbox/rasterio/blob/6c9816cb10876cd25dd5c802a52490628adc0ec1/rasterio/_warp.pyx#L940-L943 and thus ignore input transform

Mar 18 '20 18:03 vincentsarago

@vincentsarago yes I've also seen that removing the argument has the same effect.

I've naively tried to force the src_dataset's transform with this diff:

diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx
index acc4c92c..19e8d79e 100644
--- a/rasterio/_warp.pyx
+++ b/rasterio/_warp.pyx
@@ -819,6 +819,9 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
         cdef int mask_block_xsize = 0
         cdef int mask_block_ysize = 0
 
+        if self.src_transform:
+            self.src_dataset.transform = self.src_transform
+
         hds = exc_wrap_pointer((<DatasetReaderBase?>self.src_dataset).handle())
 
         if not self.src_transform:

but that doesn't work because src_dataset.transform is a read-only attribute.

I could circumvent the problem by creating a VRT of my dataset where I override it's transform, but that defeats the whole purpose of using a WarpedVRT :smiley:

Mar 18 '20 18:03 glostis

other interesting notes, the transform resulting from gdalwarp is not equal to the one we get with from_gcps

gdalwarp -ot float32 -multi -t_srs EPSG:4326 -tps iw-vv.tiff iw-vv.vrt

with rasterio.open("iw-vv.vrt") as src:
    print(src.transform)
| 0.00, 0.00, 75.61|
| 0.00,-0.00, 11.19|
| 0.00, 0.00, 1.00|

when using rasterio WarpedVRT we get

with rasterio.open("iw-vv.tiff") as src:
    with WarpedVRT(src, src_crs="epsg:4326") as vrt:
        print(vrt transform)
        print(vrt.transform)

vrt transform
| 0.00, 0.00, 75.61|
| 0.00,-0.00, 11.19|
| 0.00, 0.00, 1.00|

with rasterio.open("iw-vv.tiff") as src:
    print("transform from gcps") 
    print(from_gcps(src.gcps[0]))

transform from gcps
|-0.00,-0.00, 78.20|
| 0.00,-0.00, 10.74|
| 0.00, 0.00, 1.00|

Not sure which one is the best one though.

When passing src_transform, transform, width and height it finally have the transfrom set to the one created from gcps but that's expected (because we specifically say we want dst_transform)

path = "iw-vv.tiff"  
with rasterio.open(path) as src:
    print("source transform") 
    print(src.transform) 
    print() 
    print("transform from gcps") 
    t = from_gcps(src.gcps[0])
    print(t)
    print() 
    with WarpedVRT(src, src_transform=t, transform=t, width=28359, height=21436, src_crs=4326) as warped:
        print(f"warped.src_transform\n{warped.src_transform}\n")  
        print(f"warped.transform\n{warped.transform}\n")

source transform
| 1.00, 0.00, 0.00|
| 0.00, 1.00, 0.00|
| 0.00, 0.00, 1.00|

transform from gcps
|-0.00,-0.00, 78.20|
| 0.00,-0.00, 10.74|
| 0.00, 0.00, 1.00|

warped.src_transform
|-0.00,-0.00, 78.20|
| 0.00,-0.00, 10.74|
| 0.00, 0.00, 1.00|

warped.transform
| 0.00, 0.00, 75.61|
| 0.00,-0.00, 11.19|
| 0.00, 0.00, 1.00|

if we just use transform=t then we are back on the one calculated by defaut

path = "iw-vv.tiff"  
with rasterio.open(path) as src:
    print("transform from gcps") 
    t = from_gcps(src.gcps[0])
    print(t)
    print() 
    with WarpedVRT(src, transform=t, src_crs=4326) as warped:
        print(f"warped.src_transform\n{warped.src_transform}\n")  
        print(f"warped.transform\n{warped.transform}\n")

transform from gcps
|-0.00,-0.00, 78.20|
| 0.00,-0.00, 10.74|
| 0.00, 0.00, 1.00|

warped.src_transform
| 1.00, 0.00, 0.00|
| 0.00, 1.00, 0.00|
| 0.00, 0.00, 1.00|

warped.transform
| 0.00, 0.00, 75.61|
| 0.00,-0.00, 11.19|
| 0.00, 0.00, 1.00|

Summary of the BUG: WarpedVRT doesn't take the input src_transform so it has no impact. ~even when I pass width and height + the src_transform the output transform is similar to the one when I only pass~

Mar 18 '20 19:03 vincentsarago

The fact that from_gcps doesn't give the same result as gdalwarp is interesting indeed. Maybe the underlying GDAL function being used is not the same, or not called with the same parameters?

Concerning width and height, you need to also pass transform to WarpedVRT (which is then internally assigned to self.dst_transform) in order for the src_transform to be effectively taken into account, because of this if: https://github.com/mapbox/rasterio/blob/6c9816cb10876cd25dd5c802a52490628adc0ec1/rasterio/_warp.pyx#L904

Mar 18 '20 19:03 glostis

thanks @glostis I've updated my previous comment.

so with the file I'm using the transform created from gcps has 78.20 value while with gdalwarp it's 75.61

when checking the gcps, the value 78.20 comes from the first gcps point, which IMO make sense.

path = "iw-vv.tiff"
with rasterio.Env(GDAL_GCPS_TO_GEOTRANSFORM_APPROX_OK="FALSE"):
    with rasterio.open(path) as src:
        t = src.gcps[0][0]
        print(t)
GroundControlPoint(row=0.0, col=0.0, x=78.20390305293508, y=10.737548333634779, z=255.30090383719653, id='1', info='')
In [ ]:

edits: gdal and rasterio return the same transform when using GCPsToGeoTransform

path = "iw-vv.tiff"
with rasterio.Env(GDAL_GCPS_TO_GEOTRANSFORM_APPROX_OK="FALSE"):
    with rasterio.open(path) as src:
        t = from_gcps(src.gcps[0])
        print("rasterio")
        print(t.to_gdal())

import gdal
ds = gdal.Open("iw-vv.tiff")
print("gdal")
print(gdal.GCPsToGeoTransform(ds.GetGCPs()))

rasterio
(78.19888393140982, -8.947144038655639e-05, -1.8313004509073926e-05, 10.738373173025243, 1.7764971108527815e-05, -9.010977604836015e-05)
gdal
(78.19888393140982, -8.947144038655639e-05, -1.8313004509073926e-05, 10.738373173025243, 1.7764971108527815e-05, -9.010977604836015e-05)

this means that gdalwarp and WarpedVRT are doing something else to get the GeoTransfrom rom the gcps.

Mar 18 '20 19:03 vincentsarago

@vincentsarago GCPsToGeoTransform is defined here: https://github.com/OSGeo/gdal/blob/master/gdal/gcore/gdal_misc.cpp#L2371-L2690, but this file also contains what looks like a different implementation of more or less the same problem: https://github.com/OSGeo/gdal/blob/master/gdal/alg/gdal_crs.cpp

This could be why the two methods give different results?

But maybe we could take the GCPsToGeoTransform vs gdalwarp transform difference to a separate discussion because it is not directly relevant to the src_transform issue at hand?

Mar 19 '20 18:03 glostis

I'm labeling this as a documentation bug. WarpedVRT supports only 4 cases: https://github.com/mapbox/rasterio/pull/1936/files#diff-7519f03c7f3cca2b36a087ca08b1ffd1R879. The src_transform parameter is part of the solution for the 2nd case. If I understand yours, the following is the right usage.

with rasterio.open("s1.tiff") as src:
    with WarpedVRT(
        src,
        crs=32630,
    ) as vrt:
        ...

Sep 09 '20 21:09 sgillies

rasterio rasterio copied to clipboard

src_transform ignored by WarpedVRT?

Expected behavior and actual behavior.

Steps to reproduce the problem.

Operating system

Rasterio version and provenance

rasterio
rasterio copied to clipboard