[gdal-dev] Computing a geo-transform with gdalwarp from a set of GCPs also rectifies the input images

Thu Mar 25 11:00:45 PDT 2021

>Since GDAL's code base is quite large, it is not easy to find the appropriate location.

I have found that using Visual Studio Code is an unvaluable tool to navigate in large source codes.
To find a definition of a function one only have to right-click it and ask to go to its definition.

From: gdal-dev <gdal-dev-bounces at lists.osgeo.org> On Behalf Of Bullinger, Sebastian
Sent: Thursday, March 25, 2021 5:31 PM
To: Even Rouault <even.rouault at spatialys.com>; gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] Computing a geo-transform with gdalwarp from a set of GCPs also rectifies the input images

Hi,

@Andrew

> Have you noticed whether the output image has no-data pixels at the
> edges, or crops the original ? I would expect it to do a little of each.

The output image has no-data values - it does not crop the image - which is inline with the description by Even (see below).

@Even

> When creating an output file from scatch (either it doesn't exist, or -overwrite is used),
> gdalwarp will always produce a "north-up" image in the target CRS.
> ...
> How gdalwarp rectifies depends on content of the source dataset (GCP, geolocation arrays, RPCs),
> user switches and default behaviours. For a source dataset with GCPs, if you have < 6 GCPs,
> a linear fit will be used by default. For >= 6 GCPs, a 2nd order polynomial fit will be used.

Thank you for clarifying this - this explains the results I'm observing.

Is it possible access the transformation used to rectify the images during the gdalwarp call?
At the moment, I feel that gdalwarp acts (w.r.t. the rectification) like a black box.
Maybe it would be reasonable to (optionally) provide some information about the rectification process?

If it is not possible to get this information, could you point out the code snippet that performs the rectification step?
Since GDAL's code base is quite large, it is not easy to find the appropriate location.

Thank you for your help - and sorry for the noise.

Best regards,
Sebastian

--
Dr. Sebastian Bullinger
Department Object Recognition
Fraunhofer Institute of
Optronics, Sytem Technologies and Image Exploitation IOSB
Gutleuthausstr. 1, 76275 Ettlingen, Germany
Phone +49 7243 992-197
sebastian.bullinger at iosb.fraunhofer.de<mailto:sebastian.bullinger at iosb.fraunhofer.de>
www.iosb.fraunhofer.de<https://webmail.iosb.fraunhofer.de/owa/redir.aspx?C=GM65JhP8Gk6mqlokIktWtqcWz_hm-dIIiaYxp7LDE5w39h4r54rwzTEQR1dXSKQtdNkO601Flpk.&URL=http%3a%2f%2fwww.iosb.fraunhofer.de%2f>

________________________________
From: Even Rouault <even.rouault at spatialys.com<mailto:even.rouault at spatialys.com>>
Sent: Thursday, March 25, 2021 5:08 PM
To: Bullinger, Sebastian; gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
Subject: Re: [gdal-dev] Computing a geo-transform with gdalwarp from a set of GCPs also rectifies the input images

Hi,

When creating an output file from scatch (either it doesn't exist, or -overwrite is used), gdalwarp will always produce a "north-up" image in the target CRS. This is by design of the utility. If you just want to subset, preserving the original orientation, use gdal_translate. In update mode of an existing target dataset however, gdalwarp will use the georeferencing attached to the target (so potentially not north-up).

How gdalwarp rectifies depends on content of the source dataset (GCP, geolocation arrays, RPCs), user switches and default behaviours. For a source dataset with GCPs, if you have < 6 GCPs, a linear fit will be used by default. For >= 6 GCPs, a 2nd order polynomial fit will be used. You may also use the -tps switch to ask for thin plate splines. The -to switch can also be used for finer control: see link pointed by https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-to

Even
Le 23/03/2021 à 18:47, Bullinger, Sebastian a écrit :
Dear gdal community,

I'm currently working with a set of satellite images that are geo-registered with a set of GCPs. To compute a geo-transform from the GCPs, I've been using gdalwarp with:

gdalwarp -of GTiff path/to/original/file path/to/warped/file

While "gdalwarp" correctly computes a transformation from the GCPs, at the same it also performs some kind of rectification / north-up image transformation. I'm not sure about the applied operation.

Presumably, the reason for this is that the original GCPs reflect a transformation with skew factors.
See the following transformations of the original and the warped image.

>>> import gdal
>>> dataset = gdal.Open("path/to/original/file")
>>> print(gdal.GCPsToGeoTransform(dataset.GetGCPs()))
(-58.57294039342205, -3.991598267026281e-06, 0.0, -34.451175177186684, -3.6492767328575985e-07, 3.1987810722132996e-06)

>>> import gdal
>>> dataset = gdal.Open("path/to/warped/file")
>>> print(dataset.GetGeoTransform())
(-58.577917916461026, 3.435846284510554e-06, 0.0, -34.44693359348493, 0.0, -3.435846284510554e-06)

The "rectification" result is very convenient, since it allows to use the images with an existing satellite image segmentation pipeline (which tiles the satellite images using geo-tiles)
However, in order to overlay the segmentations (performed on the "rectified" images) with the original images, I need some information about the transformation applied during "gdalwarp".

I've searched the documentation for more information - unfortunately without success.
For example, parameters like "-novshiftgrid" do not affect this transformation.

It would be very helpful, if someone could shed some light on the used "rectification".
What kind of operation is applied? What would be the correct term to search for?
Is it "only" a decomposition of the transform? Or does "gdalwarp" something else too?
Is there a possibility to access this information, while using "gdalwarp"? (Or can this operation also be done by hand?)
Does GDAL contain any functions to invert this "rectification" step - given the corresponding transformation used by "gdalwarp"?

Best regards,
Sebastian

--
Dr. Sebastian Bullinger
Department Object Recognition
Fraunhofer Institute of
Optronics, Sytem Technologies and Image Exploitation IOSB
Gutleuthausstr. 1, 76275 Ettlingen, Germany
Phone +49 7243 992-197
sebastian.bullinger at iosb.fraunhofer.de<mailto:sebastian.bullinger at iosb.fraunhofer.de>
www.iosb.fraunhofer.de<https://webmail.iosb.fraunhofer.de/owa/redir.aspx?C=GM65JhP8Gk6mqlokIktWtqcWz_hm-dIIiaYxp7LDE5w39h4r54rwzTEQR1dXSKQtdNkO601Flpk.&URL=http%3a%2f%2fwww.iosb.fraunhofer.de%2f>

_______________________________________________

gdal-dev mailing list

gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/gdal-dev

--

http://www.spatialys.com

My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20210325/a116f86f/attachment-0001.html>