SyntaxHighlighter

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.

4 comments:

heavenlyfodder said...

Just wanted to let you know that your post helped me solve an issue I was having. In my case, I needed a way to set the GeoTIFF fields within an existing TIFF file. I'd tried using geotifcp, but it kept crashing without explanation.

Fortunately, I was able to use a variation on the technique you described. Specifically, I used GDAL to read the GeoTIFF tags from the source GeoTIFF file, then write those tags to an existing GeoTIFF file.

Jerry said...

this script worked really nicely. I am currently writing some Python code to interpolate and clip NASA TRMM precipitation data based on arbitrary shapefile input (decimal degrees or linear meter units), and output each numpy array of the TRMM data to GeoTIFFs.

Thank you so much for this little bit of code ... very useful.

most respectfully,

Gerasimos Michalitsianos (geri)

Anonymous said...

Andre a very helpful introduction to using GDAL to work with the GeoTiff file format.
I have one query though, if the lower left y coordinate is -34.54 should not the value for top left y in the SetGeoTransform parameter be -21.89
and the n-s pixel resolution -0.11?
All the best, Pete

Anonymous said...

how do I update a tag value in a geotiff? thanls