[gdal-dev] Errors in geotransform computation from ENVI rotated header?

[Even Rouault]

On mardi 18 juillet 2017 12:27:51 CEST Kevin B. McCarty wrote:
> Hi,
> 
> We are upgrading to GDAL 2.2.1 where I work.  While we greatly
> appreciate the effort to start supporting the "rotation=xxx" syntax in
> the "map_info" field of an ENVI header, it appears to me that there
> are some inconsistencies with the computation of the GDAL GeoTransform
> matrix from rotated ENVI headers.
> 
> Here is the current implementation in frmts/raw/envidataset.cpp in
> ENVIDataset::ProcessMapinfo() (copied from github, commit
> ddbf6d39aa4b005a77ca4f27c2d61a3214f336f8) which I paste for reference:
> 
>   // Capture geotransform.
>   const double xReference = CPLAtof(papszFields[1]);
>   const double yReference = CPLAtof(papszFields[2]);
>   const double pixelEasting = CPLAtof(papszFields[3]);
>   const double pixelNorthing = CPLAtof(papszFields[4]);
>   const double xPixelSize = CPLAtof(papszFields[5]);
>   const double yPixelSize = CPLAtof(papszFields[6]);
> 
>   adfGeoTransform[0] = pixelEasting - (xReference - 1) * xPixelSize;
>   adfGeoTransform[1] = cos(dfRotation) * xPixelSize;
>   adfGeoTransform[2] = -sin(dfRotation) * xPixelSize;
>   adfGeoTransform[3] = pixelNorthing + (yReference - 1) * yPixelSize;
>   adfGeoTransform[4] = -sin(dfRotation) * yPixelSize;
>   adfGeoTransform[5] = -cos(dfRotation) * yPixelSize;
> 
> [My apologies in advance for any obnoxious breakage that gmail applies
> to the formatting of this snippet.]
> 
> The element [4] of adfGeoTransform is equivalent in semantics to the
> second line [that is, zero-based line 1] of a world file (to
> illustrate this I point for instance at the implementation of
> GDALLoadWorldFile() in gcore/gdal_misc.cpp).  Per Wikipedia:
> 
>   https://en.wikipedia.org/wiki/World_file
> 
> ... "A better description of the A, D, B and E parameters are: [...]
> Line 2: D: y-component of the pixel *width* (y-skew)"  [emphasis
> mine].
> 
> So it seems to me that the line setting adfGeoTransform[4] above
> should instead read:
> 
>   adfGeoTransform[4] = -sin(dfRotation) * xPixelSize; // NOTE
> yPixelSize --> xPixelSize
> 
> Similarly, in the computation of element [2], which is semantically
> equivalent to the third line [that is, zero-based line 2] of a world
> file: as Wikipedia puts it, the "x-component of the pixel *height*"
> [emphasis mine]; xPixelSize should be replaced by yPixelSize.
> 
> The error is apparent in practice only when xPixelSize != yPixelSize,
> of course; which is presumably why this wasn't noticed earlier.
> 
> 
> An additional issue is that I think elements [0] and [3] (the
> world-coordinate X,Y position of the pixel-space origin) of the
> transform are not being computed correctly when rotation != 0.  For
> instance, the world X origin (element [0]) should first of all take
> into account that the xPixelSize is partly in the world Y direction
> (so its extent in the world X direction must be reduced by the
> appropriate angular factor); second, it should take into account that
> the yPixelSize contributes to some extent in the world X direction.
> 
> To summarize all in one place, I think that the correct computations
> would look like this:
> 
>   adfGeoTransform[1] = cos(dfRotation) * xPixelSize;
>   adfGeoTransform[2] = -sin(dfRotation) * yPixelSize;
>   adfGeoTransform[4] = -sin(dfRotation) * xPixelSize;
>   adfGeoTransform[5] = -cos(dfRotation) * yPixelSize;
> 
>   adfGeoTransform[0] = pixelEasting - (xReference - 1) * adfGeoTransform[1]
>                                                            -
> (yReference - 1) * adfGeoTransform[2];
>   adfGeoTransform[3] = pixelNorthing - (xReference - 1) * adfGeoTransform[4]
> -
> (yReference - 1) * adfGeoTransform[5];
> 
> [again, apologies in advance for however Gmail may have messed this
> up; fortunately C++ is mostly whitespace-insensitive]
> 
> Note that for element [3], the + sign on (yReference - 1) * yPixelSize
> changes to a - sign on (yReference - 1) * adfGeoTransform[5] ... this
> follows since in the zero-rotation case, yPixelSize being positive
> implies adfGeoTransform[5] is negative.
> 
> 
> It is possible that there are similar reciprocal errors in the code in
> ENVIDataset::WriteProjectionInfo() to write out an ENVI dataset header
> given a rotated GeoTransform ... I have not looked into that code
> since we don't use it here.
> 
> Best regards, and many thanks for your work,

Kevin,

Guillaume Pasero also sent an email a few days ago about issues regarding this:
https://lists.osgeo.org/pipermail/gdal-dev/2017-July/046860.html
so it looks like we have a confirmed issue here

Did you check that the georeferencing in a GDAL based software (QGIS for example, or 
maybe by just running gdalwarp to get a north-up oriented output image) with your above 
proposed fix is consistant with the one you get in ENVI software ?

I've just left a message to the original author of ENVI rotation read support in:
https://github.com/OSGeo/gdal/pull/197#issuecomment-316185610

The writing side is very much likely affected since the formulas were deduced by reversing 
the ones of the read side.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170718/ec647654/attachment-0001.html>