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

[Even Rouault]

On mardi 18 juillet 2017 14:45:31 CEST Kevin B. McCarty wrote:
> Hi Even,
> 
> On Tue, Jul 18, 2017 at 1:28 PM, Even Rouault
> 
> <even.rouault at spatialys.com> wrote:
> > 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 created a pretty artificial .hdr file that exhibits rotation
> together with xpixelsize != ypixelsize.  Here it is: (cut-and-paste it
> somewhere and name it to dum.hdr)
> 
> ENVI
> samples = 10
> lines   = 20
> bands   = 1
> header offset = 0
> file type  = ENVI Standard
> data type  = 1
> interleave = bsq
> byte order = 0
> map info = { UTM, 1, 1, 200, 200, 25, 12.5,
>              30, North, WGS-84, units=Meters, rotation=45 }
> 
> Note that the local X size (sample count) is 10 pixels and local Y
> size (line count) is 20 pixels; but since the local X pixel width is
> twice the local Y pixel width, the image in world coordinates should
> be a square of 250 meters on each side, rotated 45 degrees from due
> north about its upper left (in pixel-space) corner.
> 
> You can create an accompanying .bsq file for this .hdr file with the
> following Unix command, or via any other mechanism that creates a file
> having at least 200 bytes (exact contents don't matter):
> 
>   dd if=/dev/zero of=./dum.bsq bs=200 count=1
> 
> 
> Here is what gdal 2.2.1 gdalinfo reports for the corners of dum.bsq
> [that I claim is erroneous]:
> 
> Upper Left  (     200.000,     200.000) (  7d29'13.03"W,  0d 0' 6.49"N)
> Lower Left  (     553.553,      23.223) (  7d29' 1.62"W,  0d 0' 0.75"N)
> Upper Right (     376.777,     288.388) (  7d29' 7.33"W,  0d 0' 9.36"N)
> Lower Right (     730.330,     111.612) (  7d28'55.92"W,  0d 0' 3.62"N)
> Center      (     465.165,     155.806) (  7d29' 4.48"W,  0d 0' 5.06"N)
> 
> 
> Here is what it reports for the same input after my suggested fix for
> the transform matrix calculation:
> 
> Upper Left  (     200.000,     200.000) (  7d29'13.03"W,  0d 0' 6.49"N)
> Lower Left  (     376.777,      23.223) (  7d29' 7.33"W,  0d 0' 0.75"N)
> Upper Right (     376.777,     376.777) (  7d29' 7.33"W,  0d 0'12.23"N)
> Lower Right (     553.553,     200.000) (  7d29' 1.62"W,  0d 0' 6.49"N)
> Center      (     376.777,     200.000) (  7d29' 7.33"W,  0d 0' 6.49"N)
> 
> 
> If you view these two sets of five world coordinates each in your
> favorite map program, you can see that the former coordinates make the
> corners and center of a long thin parallelogram, whereas the latter
> coordinates form a 45-degree rotated square as expected.  And for the
> latter set of coordinates, the difference between min and max X
> coordinates, and between min and max Y coordinates, is sqrt(2) * 250
> as expected.
> 
> I've also tried some tests with our own (proprietary) software using
> each of the vanilla GDAL 2.2.1 .so file, and the version with my fix,
> in turn.
> 
> I've not tried as hard to verify my suggested fixes to computation of
> adfGeoTransform elements [0] and [3], since they are smaller in effect
> and only apply to files that are both rotated and in which the
> reference X and/or Y pixels are other than one, which are probably
> rare.  But it seems clear to me that some change to elements [0] and
> [3] of this sort is needed in that (obscure) case.
> 
> > 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
> 
> OK, thank you!

If we don't get feedback from him, please create a ticket about the issue with your changes, 
or a pull request at your convenience. Fixes in the writing side would be appreciated too. I 
guess the envi_15() test in autotest/gdrivers/envi.py will probably need some adaptations.

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/20170719/82cd7702/attachment.html>