Topics

Use of GCPS in Rasterio

vincent.sarago@...
 

Right now Rasterio has partial support of GCPS (meaning it can just list them). 

I've started a PR to add support of GCPS to retrieve the geotransform and crs from it when it's not present in the file (mostly to support Sentinel-1 data COG): https://github.com/mapbox/rasterio/pull/1749
While this works, the implementation and logic might not be the one expected from users. Because it's unlikely that a file will have a georeferenced coordinates system + gcps, I'd love if we could have something at least semi-automatic/seamless for users when opening the file (e.g with an environment variable ?) to handle gcps information.

GDAL docs says "The GDAL data model does not imply a transformation mechanism that must be generated from the GCPs … this is left to the application"
Ref: https://gdal.org/user/raster_data_model.html?highlight=gcps#gcps

GDAL behaviour:
- gdalinfo: list gcps, does not get transform, bounds and crs from GCPS
- gdal_translate: create a pure copy, does not get transform, bounds and crs from GCPS
- gdalwarp: get transform, bounds and crs from GCPS and copy those to the output file

Rasterio behaviour:
- rio info: list gcps, does not get transform, bounds and crs from GCPS
- rio convert: copy the file, loose GCPS info, does not get transform, bounds and crs from GCPS
- rio warp: copy the file, loose GCPS info, does not get transform, bounds and crs from GCPS

Note: 
when using `WarpedVRT`rasterio seems to be able to get the geotransform from the gcps 

```
with rasterio.open(src_path) as src_dst: 
    with WarpedVRT(src_dst, src_crs=src_dst.gcps[1], scrs=src_dst.gcps[1]) as vrt_dst: 
        print(vrt_dst.bounds)
        print(vrt_dst.meta)
        print(vrt_dst.transform)
```


Sean Gillies
 

Hi Vincent,

Doesn't GDAL get georeferencing for GeoTIFFs from GCPs already? https://gdal.org/drivers/raster/gtiff.html#georeferencing leads me to believe that it does. GDAL doesn't get georeferencing for a VRT from its GCPs, I have just checked this. Maybe the place to address this is in GDAL?

On Thu, Aug 22, 2019 at 8:16 AM <vincent.sarago@...> wrote:
Right now Rasterio has partial support of GCPS (meaning it can just list them). 

I've started a PR to add support of GCPS to retrieve the geotransform and crs from it when it's not present in the file (mostly to support Sentinel-1 data COG): https://github.com/mapbox/rasterio/pull/1749
While this works, the implementation and logic might not be the one expected from users. Because it's unlikely that a file will have a georeferenced coordinates system + gcps, I'd love if we could have something at least semi-automatic/seamless for users when opening the file (e.g with an environment variable ?) to handle gcps information.

GDAL docs says "The GDAL data model does not imply a transformation mechanism that must be generated from the GCPs … this is left to the application"
Ref: https://gdal.org/user/raster_data_model.html?highlight=gcps#gcps

GDAL behaviour:
- gdalinfo: list gcps, does not get transform, bounds and crs from GCPS
- gdal_translate: create a pure copy, does not get transform, bounds and crs from GCPS
- gdalwarp: get transform, bounds and crs from GCPS and copy those to the output file

Rasterio behaviour:
- rio info: list gcps, does not get transform, bounds and crs from GCPS
- rio convert: copy the file, loose GCPS info, does not get transform, bounds and crs from GCPS
- rio warp: copy the file, loose GCPS info, does not get transform, bounds and crs from GCPS

Note: 
when using `WarpedVRT`rasterio seems to be able to get the geotransform from the gcps 

```
with rasterio.open(src_path) as src_dst: 
    with WarpedVRT(src_dst, src_crs=src_dst.gcps[1], scrs=src_dst.gcps[1]) as vrt_dst: 
        print(vrt_dst.bounds)
        print(vrt_dst.meta)
        print(vrt_dst.transform)
```
,_._,_

--
Sean Gillies

vincent.sarago@...
 

We added the `from_gcps` function to create a transformation matrix from gcps (see: https://github.com/mapbox/rasterio/pull/1773
This is a good first step but now when we open a file we still need to do something like 

```
with rasterio.open(src_path) as src_dst: 
    with WarpedVRT(src_dst, src_crs=src_dst.gcps[1], src_transform=from_gcps(src_dst.gcps[0])) as vrt_dst: 
        ...
```

Question: I wonder if it will be possible to add an option (or env variable) for the `open` method to force the use of crs and transform created from gcps. 

Sean Gillies
 

Hi Vincent,

I think we can do this without a flag or environment variable if we agree on a convention. How about: when creating the GDAL transformer in the WarpedVRT class, we use the source dataset's transform *or* GCPs transform, unless the user provides the value?

I'll do this right now if you agree.


On Tue, Sep 10, 2019 at 7:29 AM <vincent.sarago@...> wrote:
We added the `from_gcps` function to create a transformation matrix from gcps (see: https://github.com/mapbox/rasterio/pull/1773
This is a good first step but now when we open a file we still need to do something like 

```
with rasterio.open(src_path) as src_dst: 
    with WarpedVRT(src_dst, src_crs=src_dst.gcps[1], src_transform=from_gcps(src_dst.gcps[0])) as vrt_dst: 
        ...
```

Question: I wonder if it will be possible to add an option (or env variable) for the `open` method to force the use of crs and transform created from gcps.

--
Sean Gillies

vincent.sarago@...
 

Sean, my main worry is that we need to use `WarpedVRT` so it adds one more layer, and thus not optimal. But it will be better than nothing :-) 

For files with gcps, not only the transform but also the CRS won't be set by default so basically, the WarpedVRT reader it should do something like: 

```
gcps, gcps_crs = src.gcps
in_transform = from_gcps(gcps) if gcps else src.transform
in_crs = gcps_crs if gcps else src.crs

# in https://github.com/mapbox/rasterio/blob/master/rasterio/_warp.pyx#L738
self.src_transform = src_transform if src_transform else in_transform

# in https://github.com/mapbox/rasterio/blob/master/rasterio/_warp.pyx#L736
self.src_crs = src_crs if src_crs else in_crs
```

Sean Gillies
 

Vincent.

After a little more thought, I'm thinking that working around GDAL with respect to geotransform and GCPs is not a good idea. GDALGetGeoTransform() does not use GCPs if they exist and the second to last paragraph in https://gdal.org/user/raster_data_model.html#gcps cautions that the meaning of the GCPs is up to the application. And by application, it means the application of the data defined by the provider, not our software project as an application. It could be the case that GCPs in one file are meant to be used differently than GCPs in another file, and if we made Rasterio overly clever we might get in the way of programmers or cause problems with the wrong default behavior.

On Tue, Sep 10, 2019 at 10:52 AM <vincent.sarago@...> wrote:
Sean, my main worry is that we need to use `WarpedVRT` so it adds one more layer, and thus not optimal. But it will be better than nothing :-) 

For files with gcps, not only the transform but also the CRS won't be set by default so basically, the WarpedVRT reader it should do something like: 

```
gcps, gcps_crs = src.gcps
in_transform = from_gcps(gcps) if gcps else src.transform
in_crs = gcps_crs if gcps else src.crs

# in https://github.com/mapbox/rasterio/blob/master/rasterio/_warp.pyx#L738
self.src_transform = src_transform if src_transform else in_transform

# in https://github.com/mapbox/rasterio/blob/master/rasterio/_warp.pyx#L736
self.src_crs = src_crs if src_crs else in_crs

--
Sean Gillies