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 <[email protected]> _______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
