SyntaxHighlighter

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.

Thursday, December 4, 2008

Grep and punctuation marks

I've often tried to use grep to filter search results for locate and several times I've had to search for a string with locate:
locate svn
This search would obviously include all the '.svn' folders too. Trying to get rid with them with
locate svn | grep -v '.svn'
does not work, since grep uses basic regular expressions, and the period mark will simply match any character.

To explicitly indicate the usage of a punctuation mark, one has to put it in brackets. The statement will now look something like:
locate svn | grep -v '[.]svn'

I just scanned through the grep man pages, and everything is explained pretty clearly there. I'd recommend reading the bit on regular expressions.

Wednesday, December 3, 2008

MySQL and the joys of security

I've been running my local MySQL server on Gentoo without hassles for a while and recently tried to connect to it from a remote machine. Obviously there was bound to be some issues and here's how I solved them.

When you first try to connect remotely to the machine, you'll be greeted by this lovely message:
# mysql --user=root --password=some_password --host=192.168.0.1 MYDB
ERROR 2003 (HY000): Can't connect to MySQL server on '
192.168.0.1' (111)
To get this sucker out of the way, edit your /etc/mysql/my.cnf file and comment the following line:
# bind-address = 127.0.0.1
Note to all the security girls: Your server will now be accessible from any other networked machine. Oooohh nooooo!!

Anyhow, now we try to connect again, and whooopa!! Another error:
mysql --user=root --password=some_password --host=192.168.0.1 MYDB
ERROR 1130 (00000): Host 'remote_pc_hostname' is not allowed to connect to this MySQL server
Ah bite me. Ok, so you have to grant his Lordship permissions to connect to thy server. In order to do that logon to sql on the sql server itself:
# mysql --user=root --password=some_password
Welcome to the MySQL monitor. Commands end with ; or \g.
Your MySQL connection id is 11
Server version: 5.0.70-log Gentoo Linux mysql-5.0.70-r1

Type 'help;' or '\h' for help. Type '\c' to clear the buffer.

mysql> GRANT ALL ON *.* TO 'root'@'remote_pc_hostname';
Query OK, 0 rows affected (0.00 sec)
Now all should be well in fairyland, but alas. The code securo freaks over at Sun have presented ye with another trial.
# mysql --user=root --password=some_password --host=192.168.0.1 MYDB
ERROR 1045 (28000): Access denied for user 'root'@'remote_pc_hostname' (using password: YES)
Ok I'm taking the fairytale rubbish too far. Here's what you should've done on the sql server in the first place:
# mysql --user=root --password=some_password
Welcome to the MySQL monitor. Commands end with ; or \g.
Your MySQL connection id is 11
Server version: 5.0.70-log Gentoo Linux mysql-5.0.70-r1

Type 'help;' or '\h' for help. Type '\c' to clear the buffer.

mysql> GRANT ALL ON *.* TO 'root'@'remote_pc_hostname' IDENTIFIED BY 'some_password';
Query OK, 0 rows affected (0.00 sec)

Sheesh! You little girls are wasting my time. Take a que from MS SQL Server.

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)