Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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)
share|improve this question

1 Answer 1

up vote 4 down vote accepted

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 ]
share|improve this answer
    
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. –  iant Apr 24 '12 at 10:16
    
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/… –  monkut Apr 25 '12 at 0:41

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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