6

I am developing an application with the GDAL library to combine georeferenced raster hyperspectral imagery in ENVI format from NASA's Airborne Visual/Infrared Imaging Spectrometer (AVIRIS) with other geographic data sets. When I import the imagery into a GIS application such as QGIS, I find that it is rotated incorrectly:

Incorrectly rotated flightline in QGIS

The flightline should go from southwest to northeast, but it is displayed from south to north. The map info field of the ENVI header file contains a rotation parameter of -66° which I assume is relative to the top left pixel:

$ grep "map info" ang20150422t163638_corr_v1e_img_zero.hdr
map info = { UTM , 1.000 , 1.000 , 736600.089 , 4078126.750 , 2.7000000000e+00 , 2.7000000000e+00 , 12 , North , WGS-84 , units=Meters , rotation=-66.00000000 }

However, this field is not documented and is apparently ignored by GDAL. Here is the current projection of the image:

$ gdalsrsinfo ang20150422t163638_corr_v1e_img_zero.img

PROJ.4 : '+proj=utm +zone=12 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

OGC WKT :
PROJCS["UTM Zone 12, Northern Hemisphere",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-111],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

What I would like to do is rotate the image using either the GDAL utilities or preferably the GDAL Python API so that it is oriented north-up without changing anything else about the image.

I have read that rotations can be achieved with some combination of gdalwarp and gdal_translate, either by applying a geotransform or changing the +towgs84 parameters.

However, I am a GIS novice and seek specific guidance on what parameters to set to accomplish this.

You may find a single-band version of the image and its header here.

First Try

Based on lpdudley's suggestion, I made a copy of the image:

$ cp ang20150422t163638_corr_v1e_img_zero.hdr ang20150422t163638_corr_v1e_img_zero_rot.hdr
$ cp ang20150422t163638_corr_v1e_img_zero.img ang20150422t163638_corr_v1e_img_zero_rot.img

And then applied the geotransform in the Python shell:

>>> image_name = "ang20150422t163638_corr_v1e_img_zero_rot.img"
>>> dataset = gdal.Open(image_name, gdal.GA_Update)
>>> gt = dataset.GetGeoTransform()
>>> rotation = (math.pi/180) * -66
>>> new_geotransform = (gt[0],
                        math.cos(rotation) * gt[1],
                        -math.sin(rotation) * gt[1],
                        gt[3],
                        math.sin(rotation) * gt[5],
                        math.cos(rotation) * gt[5])
... ... ... ... ... >>> dataset.SetGeoTransform(new_geotransform)
0
>>> dataset = None

However, the image is scaled rather than rotated: Incorrectly scaled flightline in QGIS

Bug

Rotating ENVI seems to be a bug in GDAL:

The bug appears to come from the ENVI driver. If I find a solution I'll share it here.

Fix

I have posted a pull request and shared it on the GDAL mailing list using the values suggested by lpdudley. Now when I apply gdalwarp to the image:

$ gdalwarp -of envi ang20150422t163638_corr_v1e_img_zero.img ang20150422t163638_corr_v1e_img_zero_rot.img
Creating output file that is 4048P x 2461L.
Processing input file ang20150422t163638_corr_v1e_img_zero.img.
0...10...20...30...40...50...60...70...80...90...100 - done.

The rotation will be detected correctly:

Correctly rotated flightline in QGIS

Rotated ENVI files will be supported in GDAL version 2.2.0 or greater.

5

The rotation can be changed from the python API.

The raster grid size, position and rotation parameters can be accessed with the GetGeoTransform() method, and they can be changed with the SetGeoTransform method. The geotransform is a six element tuple:

(x_origin,
cos(rotation) * x_pixel_size,
-sin(rotation) * x_pixel_size,
y_origin,
sin(rotation) * y_pixel_size,
cos(rotation) * y_pixel_size)

rotation is in radians, not degrees.

So for your example, the code could be:

import math
from osgeo import gdal

# make sure to open in update mode
dataset = gdal.Open(image_name, gdal.GA_Update)

gt = dataset.GetGeoTransform()

# convert -66 degrees to radians
rotation = (math.pi/180) * -66 

new_geotransform = (gt[0],
                    math.cos(rotation) * gt[1],
                    -math.sin(rotation) * gt[1],
                    gt[3],
                    math.sin(rotation) * gt[5],
                    math.cos(rotation) * gt[5])

# should return 0
dataset.SetGeoTransform(new_geotransform)

dataset = None
  • I appreciate your comment. I had neglected to open the image in update mode. However, when I load the new image into QGIS it is scaled rather than rotated: drive.google.com/open?id=0BxysdOuBmaIGRk1ZQ3QzeUdtZnc . – Taylor Alexander Brown Feb 26 '17 at 11:09
  • I updated my question to show the steps I used to get that result. – Taylor Alexander Brown Feb 26 '17 at 11:17
  • I am not sure, but I think for some reason the HFA driver does not like floats with many decimal places in the geotransform array. For me it works fine with a gtiff, and fine with an .img if I round the geotransform elements to 10 d.p. or so before I set it, but not with a .img if I do not round these numbers. So I would recommend you give this a go! – lpd2 Feb 26 '17 at 12:09
  • Thanks for the suggestion. There was no effect when I rounded the numbers. Your comment about the driver made me realize I was using an old version of the library, so I installed the latest stable release from source. However, the image is still scaled and not rotated when I load it into QGIS. Recall that the image is in ENVI format. – Taylor Alexander Brown Feb 26 '17 at 23:01
  • It turned out to be a bug in GDAL. I have updated my post with the solution. Thanks for the help! – Taylor Alexander Brown Feb 27 '17 at 12:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.