SyntaxHighlighter

Showing posts with label Python. Show all posts
Showing posts with label Python. Show all posts

Friday, December 12, 2008

Python, GDAL: Set JPEG quality values and other att

Here I shall show how to set the quality values of the JPEG format in GDAL. This post will build upon this previous article which describes how to create a basic JPEG.

It is pretty trivial to set the quality (and other) JPEG values. Here is a list of parameters one can set for the JPEG format. These parameters are simply passed to the CreateCopy() method as a list of options. Here's an example:
# ... Some other code ... #

if (cmp(output_file_ext.lower(),".jpg") == 0 or cmp(output_file_ext.lower(),".jpeg") == 0):
# Add options
options = []
options.append("QUALITY=90")

jpg_driver = gdal.GetDriverByName("JPEG")

# Pass options as parameter
jpg_driver.CreateCopy( output_file, dst_ds, options=options )
os.remove(output_file_tmp)

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.

Wednesday, December 3, 2008

Creating JPG's and other formats with Python and GDAL

Ok for some daft reason the driver.Create() method is not supported by all image drivers, but the driver.CreateCopy() is. The main difference being that Create() actually creates the format from scratch whereas CreateCopy() makes a copy of another existing image. I'm sure they have their reasons, but I don't understand why everything isn't created in a generic image type and then converted to the destintation image format through CreateCopy().

Anyhow, this is how I created a JPEG with GDAL. The same approach can be followed for multiple image formats. See this link for a complete list of supported formats.
# Import libs
import numpy, os
from osgeo import osr, gdal

# Set file vars
output_file = "out.jpg"
output_file_root = os.path.splitext(output_file)[0]
output_file_ext = os.path.splitext(output_file)[1]
output_file_tmp = output_file_root + ".tmp"

# Create tmp gtif
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file_tmp, 512, 512, 1, gdal.GDT_Byte )
raster = numpy.zeros( (512, 512) )
dst_ds.GetRasterBand(1).WriteArray(raster)

# Create jpeg or rename tmp file
if (cmp(output_file_ext.lower(),"jpg" ) == 0 or cmp(output_file_ext.lower(),"jpeg") == 0):
jpg_driver = gdal.GetDriverByName("JPEG")
jpg_driver.CreateCopy( output_file, dst_ds, 0 )
os.remove(output_file_tmp)
else:
os.rename(output_file_tmp, output_file)

Create a GeoTiff with Python and GDAL

It looks like the GDAL tutorial on creating GeoTiff's is slightly outdated (2008/05). I received this error with Python 2.5.2, GDAL 1.53 and NumPy 1.0.4:
Traceback (most recent call last):
File "test.py", line 9, in <module>
dst_ds.GetRasterBand(1).WriteArray(raster)
File "/usr/lib/python2.5/site-packages/gdal.py", line 796, in WriteArray

File "usr/lib/python2.5/site-packages/osgeo/gdal_array.py", line 154, in BandWriteArray
AttributeError: dtype

As far as I can tell this occurs because the matrix raster doesn't have the attribute (property) dtype. Something must have changed, but I managed to get it going by using the numpy lib rather than Numeric.

The code that worked for me looks like this:
import numpy
from osgeo import osr, gdal
format = "GTiff"
driver = gdal.GetDriverByName( format )
dst_ds = driver.Create("out.tif", 512, 512, 1, gdal.GDT_Byte )
raster = numpy.zeros( (512, 512) )
dst_ds.GetRasterBand(1).WriteArray(raster)

One can also explicitly state the matrix's data type by adding the dtype attribute:
raster = numpy.zeros( (512, 512), dtype = numpy.uint8)

Wednesday, July 30, 2008

Python: Multi-dimensional dynamic arrays

Another very poorly documented topic is how to create dynamic arrays of type array. Every second website will show you how to create static 2D arrays, which is pretty much useless for 90% of the practical applications out there.

I 'solved' it to a degree by creating a list of array items, but this is still not ideal. I would like to use a pure multi-deminsional array of type array.

Here's how I did it:
multi_arr = []

for i in range(some_range):
arr = array('f') # array of type float

# ... populate arr ...

multi_arr.append(arr)

Julian Day in Python

I found very little documentation on obtaining the Juilian day with the time module in Python. So this is how it's done:
import time

print time.strftime("%j")

As easy as that!