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 have an array of data, and for each datapoint I know the latitude and longitude. I'd like to save it as a GTiff with the same projection as other rasters I have. This is what I've tried so far, but no luck.

import numpy as np
import gdal
from gdalconst import *
from osgeo import osr

def GetGeoInfo(FileName):
    SourceDS = gdal.Open(FileName, GA_ReadOnly)
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())    
return GeoT, Projection

def CreateGeoTiff(Name, Array, driver, 
                  xsize, ysize, GeoT, Projection):
    DataType = gdal.GDT_Float32
    NewFileName = Name+'.tif'
    # Set up the dataset
    DataSet = driver.Create( NewFileName, xsize, ysize, 1, DataType )
            # the '1' is for band 1.
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    # Write the array
    DataSet.GetRasterBand(1).WriteArray( Array )
    return NewFileName

def ReprojectCoords(x, y,src_srs,tgt_srs):
    trans_coords=[]
    transform = osr.CoordinateTransformation( src_srs, tgt_srs)
    x,y,z = transform.TransformPoint(x, y)
    return x, y

# Some Data
Data = np.random.rand(5,6)
Lats = np.array([-5.5, -5.0, -4.5, -4.0, -3.5])
Lons = np.array([135.0, 135.5, 136.0, 136.5, 137.0, 137.5])

# A raster file that exists in the same approximate aregion.
RASTER_FN = 'some_raster.tif'

# Open the raster file and get the projection, that's the
# projection I'd like my new raster to have, it's 'projected',
# i.e. x, y values are numbers of pixels.
GeoT, TargetProjection, DataType = GetGeoInfo(RASTER_FN)
# Meanwhile my raster is currently in geographic coordinates.
SourceProjection = TargetProjection.CloneGeogCS()

# Get the corner coordinates of my array
LatSize, LonSize = len(Lats), len(Lons)
LatLow, LatHigh = Lats[0], Lats[-1]
LonLow, LonHigh = Lons[0], Lons[-1]
# Reproject the corner coordinates from geographic
# to projected...
TopLeft = ReprojectCoords(LonLow, LatHigh, SourceProjection, TargetProjection)
BottomLeft = ReprojectCoords(LonLow, LatLow, SourceProjection, TargetProjection)
TopRight = ReprojectCoords(LonHigh, LatHigh, SourceProjection, TargetProjection)
# And define my Geotransform
GeoTNew = [TopLeft[0],  (TopLeft[0]-TopRight[0])/(LonSize-1), 0,
           TopLeft[1], 0, (TopLeft[1]-BottomLeft[1])/(LatSize-1)]

# I want a GTiff
driver = gdal.GetDriverByName('GTiff')
# Create the new file...
NewFileName = CreateGeoTiff('Output', Data, driver, LatSize, LonSize, GeoTNew, TargetProjection)

But this results in the following error message:

File "TES2GtifBBB.py", line 25, in CreateGeoTiff
  DataSet.GetRasterBand(1).WriteArray( Array )
File "/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site packages/osgeo/gdal.py", line 1082, in WriteArray
  return gdalnumeric.BandWriteArray( self, array, xoff, yoff )
File "/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site packages/osgeo/gdal_array.py", line 256, in BandWriteArray
  raise ValueError("array larger than output file, or offset off edge")
ValueError: array larger than output file, or offset off edge
share|improve this question

1 Answer 1

up vote 2 down vote accepted

Eddy, try a different approach to write the array as a raster, the following code works for me and for your case you would only need to change the corner coordinate, prixel size and number of pixels for the reprojected values as you did in your code, so the first part would be replaced by some of your procedures.

from osgeo import gdal

def array_to_raster(array):
    """Array > Raster
    Save a raster from a C order array.

    :param array: ndarray
    """
    dst_filename = '/a_file/name.tiff'


    # You need to get those values like you did.
    x_pixels = 16  # number of pixels in x
    y_pixels = 16  # number of pixels in y
    PIXEL_SIZE = 3  # size of the pixel...        
    x_min = 553648  
    y_max = 7784555  # x_min & y_max are like the "top left" corner.
    wkt_projection = 'a projection in wkt that you got from other file'

    driver = gdal.GetDriverByName('GTiff')

    dataset = driver.Create(
        dst_filename,
        x_pixels,
        y_pixels,
        1,
        gdal.GDT_Float32, )

    dataset.SetGeoTransform((
        x_min,    # 0
        PIXEL_SIZE,  # 1
        0,                      # 2
        y_max,    # 3
        0,                      # 4
        -PIXEL_SIZE))  

    dataset.SetProjection(wkt_projection)
    dataset.GetRasterBand(1).WriteArray(array)
    dataset.FlushCache()  # Write to disk.
    return dataset, dataset.GetRasterBand(1)  #If you need to return, remenber to return  also the dataset because the band don`t live without dataset.
share|improve this answer
    
Great, thanks. I've adapted it and it works a charm! –  EddyThe B Apr 19 '13 at 18:49

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.