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

Kevin B. McCarty kmccarty at gmail.com
Tue Jul 18 12:27:51 PDT 2017


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 B. McCarty
<kmccarty at gmail.com>


More information about the gdal-dev mailing list