GDAL – Rotate ENVI Hyperspectral Imagery Using GDAL

coordinate systemenvigdalrasterremote sensing

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.

Best Answer

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
Related Question