Thursday, December 11, 2008

Python, GDAL: Adding GeoTiff Meta Data

I'm pretty much a noob when it comes coordinate systems and projections. I only know the standard GPS coordinates which are in the latitude/longitude (lat/long) format. It turns out that this system is known as the WGS84 (World Geodetic System) format.

The project I'm busy with at the moment, generates GeoTiff images and I only had a raster with the bottom lower left corner coordinates and the pixel sizes in km. It also happened that this is the only information I needed to generate the GeoTiff.

So this is the info I have:
ncols         174
nrows 115
xllcorner 14.97
yllcorner -34.54
cellsize 0.11
It's just the number of x and y pixels, the lower left corner and the pixel size.

Now we can go ahead and generate the GeoTiff like so:
if __name__ == '__main__':
# Import libs
import numpy, os
from osgeo import osr, gdal

# Set file vars
output_file = "out.tif"

# Create gtif
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, 174, 115, 1, gdal.GDT_Byte )
raster = numpy.zeros( (174, 115) )

# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform( [ 14.97, 0.11, 0, -34.54, 0, 0.11 ] )

# set the reference info
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
dst_ds.SetProjection( srs.ExportToWkt() )

# write the band
dst_ds.GetRasterBand(1).WriteArray(raster)
Easy as pie with GDAL.

0 comments: