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

Kevin B. McCarty kmccarty at gmail.com
Tue Jul 18 14:45:31 PDT 2017


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!


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

I was afraid of that.  Thanks again.

-- 
Kevin B. McCarty
<kmccarty at gmail.com>


More information about the gdal-dev mailing list