Raster data array output flipped on x-axis using python/gdal?

Question

I'm trying to create a raster using the python gdal libraries and I've gotten to the point where data is being output, but the output data is flipped on the origin point's x-axis. I know I must be overlooking something, but I can't figure out where I'm going wrong. Any ideas?

When creating the raster I set the upper-left x/y values, and the array appears to be indexed from the upper-left and continue down to lower-right. In the code below I'm filling the array with the value of the row.

When printing the array out it looks like this:

[[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1.]
 [ 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.
 2. 2. 2. 2. 2. 2.]
 [ 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3.
 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3.
 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3. 3.
 3. 3. 3. 3. 3. 3.]
 [ 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.
 4. 4. 4. 4. 4. 4.]
...

And this data writes successfully to the raster band. However when viewed in MapWindow GIS, the data appears to go in the opposite direction with the originally set upper-left origin point, appearing as the lower-left value.

In other words, the data is flipped on the origin point's x-axis.

import gdal
import osr
import numpy
OUTPUT_FORMAT = "GTiff"
def create_raster(filename="test.tif"):
 driver = gdal.GetDriverByName(OUTPUT_FORMAT)
 band_type = gdal.GDT_Byte
 number_of_bands = 1
 x_rotation = 0 # not supported
 y_rotation = 0 # not supported
 cell_width_meters = 50
 cell_height_meters = 50
 (min_lon, min_lat, max_lon, max_lat) = _get_point_bounds() # retrieve bounds for point data 
 srs = osr.SpatialReference()
 srs.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon 
 srs.SetUTM( 54, True) # Set projected coordinate system to handle meters 
 # create transforms for point conversion
 wgs84_coordinate_system = srs.CloneGeogCS() # clone only the geographic coordinate system
 wgs84_to_utm_transform = osr.CoordinateTransformation(wgs84_coordinate_system, srs)
 # convert to UTM
 top_left_x, top_left_y, z = wgs84_to_utm_transform.TransformPoint(min_lon, max_lat, 0) 
 lower_right_x, lower_right_y, z = wgs84_to_utm_transform.TransformPoint(max_lon, min_lat, 0) 
 cols, rows = _get_raster_size(top_left_x, lower_right_y, lower_right_x, top_left_y, cell_width_meters, cell_height_meters)
 dataset = driver.Create(filename, cols, rows, number_of_bands, band_type) #
 # GeoTransform parameters
 # --> need to know the area that will be covered to define the geo tranform
 # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
 geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, cell_height_meters ]
 dataset.SetGeoTransform(geo_transform)
 dataset.SetProjection(srs.ExportToWkt())
 dataset_band = dataset.GetRasterBand(1)
 data = dataset_band.ReadAsArray(0, 0, cols, rows).astype(numpy.float32) # returns empty array 
 for row in xrange(rows):
 for col in xrange(cols):
 data[row][ col] = row + 1
 dataset_band.WriteArray(data, 0, 0)
 dataset_band.SetNoDataValue(0.0)
 dataset_band.FlushCache()
 dataset = None # Close file

I've also noticed when I calculate the pixel position for a given lat/lon the y-value results in a negative index, which seems sort-of correct considering that the array is from upper-left to lower-right.

inverse_geo_transform = gdal.InvGeoTransform(self.geo_transform)[1] # for mapping lat/lon to pixel
pixel_x, pixel_y = gdal.ApplyGeoTransform(self.inverse_geo_transform, utm_x, utm_y)

Answer 1

I found the problem....

The problem is in defining the geo_transform. I had the following:

x_rotation = 0 
y_rotation = 0 
cell_width_meters = 50
cell_height_meters = 50
geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, cell_height_meters ]
dataset.SetGeoTransform(geo_transform)

The Gdal documentation isn't real clear about what these values are. (See SetGeoTransform) Searching around the internets I derived that the passed values should be (in order):

• top_left_x
• cell_width_meters
• x_rotation
• top_left_y
• y_rotation
• cell_height_meters

Which seems right, BUT re-reviewing the GDAL API Tutorial I noticed that the last value, cell_height_meters was shown being given in a negative value. It seems that this was all that was needed to properly output the data in the expected orientation.

So now I've changed the geo_transform definition line to:

(Notice the added "-")

geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, -cell_height_meters ]

• Ian Turton: this is the traditional way of dealing with the image worlds like of origins in the top left and geography's way of using the lower left.
• monkut: It makes sense once you know, but approaching the problem from the code examples it's difficult to pick up the reasoning. I've found that the ArcGIS documentation has some great documentation explaining about rasters: webhelp.esri.com/arcgisdesktop/9.3/…