Not long ago I wrote about the raster calculator in GDAL. It is a useful tool and has enough functionality to solve most tasks. But as it turns out, it is not without its bugs.
I was playing with this calculator and looked at the code. Almost immediately I noticed that it reads data “as is”, which was something to worry about, as I had already stepped on that rake some time ago (/me is immersed in memories for a few minutes. It was a fun time: two parallel projects, both related to GDAL; lots of new stuff to learn…). A little test confirmed my fears. Here is the result of an NDVI calculation for a Byte raster, performed with RasterCalc and saved as Float32.
NDVI calculated with RasterCalc for inputs with Byte data type
And this is the result of the same operation performed by gdal_calc
NDVI calculated with gdal_calc for inputs with Byte data type
You can’t get this across with an image, but the white areas aren’t really white at all — they contain a NODATA value.
What can I say, the results are very interesting. But the explanation is simple: since the data is read “as is”, when performing arithmetic operations, you can run into things like integer overflow and rounding of the result (in Python, the result of dividing two integers is also an integer).
Obviously, this is a bug. So let’s write to the author, report the issue and propose a patch.
If you are using this tool, you may not want to wait for the next version of GDAL and fix the bug yourself. To do this, open gdal_calc.py in your favourite text editor and lines 230–232
Having a raster calculator in QGIS is cool, but if you need to process a lot of rasters, it won’t be much help. There is no batch mode in the calculator (neither in the plugin nor in the QGIS core). Well, there is an entry about it in my TODO. But knowing that won’t help if you have hundreds of images for which you need to calculate, e.g., NDVI, or perform pixel replacement by a tricky condition.
This is where GDAL and Python come to the rescue. It literally takes 10 minutes to write a script, as there is nice and detailed documentation. But again, writing almost identical scripts for every task is not very practical.
Let me tell you a little secret. GDAL 1.8.0 has a wonderful tool, modestly named gdal_calc.py. This small (about 300 lines) Python script works with rasters of the same size (no check for mismatching CRS is performed) and supports basic arithmetic and logical operations. It is easy to use:
# sum of two rastersgdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"# average of two rastersgdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2"# difference of raster bandsgdal_calc.py -A input.tif --A_band=1 -B input.tif --B_band=2 --outfile=result.tif --calc="A-B"
It accepts up to 26 images as input, which should be sufficient for most use cases. As you can see from the examples, there is support for parentheses, you can access individual bands, and the console nature makes it easy to use the script for batch processing. It still lacks some features, but even in its current state, it is a great console calculator.
Alessandro Furieri has announced the SpatiaLite Cookbook, an excellent guide to SpatiaLite with many practical examples, optimisation tips and information on using SpatiaLite with different programming languages. You can check it out here.
Frank Warmerdam has announced that OSGeo4W is moving to the recently released GDAL 1.8.0. Instead of GDAL 1.5.4, the default version will be 1.8.0, additional packages like gdal-python, gdal-autotest, gdal-ecw etc. will also be updated or removed if no longer needed.
For compatibility with some packages that still require GDAL 1.5.4, a new package gdal15dll (with the necessary libraries) will be added. This package will be automatically installed if needed.
Since the transition to a new version of GDAL is a rather laborious and time-consuming process, the update of some packages included in OSGeo4W is temporarily suspended. In particular, this affects QGIS nightly builds.
The maintainers of the QGIS OSGeo4W package have already announced that they will rebuild the package with GDAL 1.8.0 as soon as possible, GRASS will be rebuilt after GRASS 6.4.1 is released.
Frank Warmerdam announced the release of GDAL 1.8.0.
GDAL is a free translation library for raster and vector geospatial formats. The library provides applications with a single generic data model for all supported formats. In addition to the library itself, GDAL includes a set of powerful command line utilities for data translation and processing.
The list of changes is impressive, and the adoption of RFC 29 and RFC 30 is particularly pleasing. For example, the changes described in RFC 29 will significantly improve performance and reduce memory consumption when processing vector data (about 1.5-2.5 times, according to tests). So far, SetIgnored is only implemented for shapefiles (by Martin) and SpatiaLite (my patch).
Once threading_branch is merged into trunk, QGIS will open large shapefiles much faster and many geoprocessing scripts will speed up. By the way, I will have to rewrite my data extraction script to take advantage of these changes.
I am looking forward to adding GDAL 1.8.0 to OSGeo4W, especially as Frank promised a global update of OSGeo4W after the release. GDAL 1.8.0 will be the main version instead of the outdated 1.5.4, and Python 2.5 will be replaced by 2.7.
Last night, I submitted a patch to add SetIgnored support to the SQLite driver. Literally an hour later, it was applied. This means that the two most important drivers have already been made to work faster.
Frank Warmerdam today applied changes proposed by Martin Dobias in RFC 29: OGR Set Ignored Fields. Martin’s patch was developed as part of the GSoC project and aims to improve the speed of reading data from OGR datasources. The new functionality is currently only available in the shapefile driver. I think we should expect this functionality to appear in other drivers soon.
For QGIS, this news marks the beginning of a new stage in its development, as the inclusion of the patch in the GDAL code, along with other improvements made as part of the QGIS on Steroids GSoC project, will significantly increase performance when rendering vector layers.
Imagine you have a vector layer containing polygons, both single- and multipart. Your task is to get vertices of each geometry, including all rings and parts. Here is how it can be done using GDAL and Python.
Let’s start with simplest case — singlepart polygons. In that case everything is straightforward — we need to find out the number of vertices, iterate over them and get coordinated. For demostration purposes the code will look like this:
defdump_geometry(geom):
print "Exterior ring" ring = geom.GetGeometryRef(0)
for i in range(ring.GetPointCount()):
print i, ring.GetX(i), ring.GetY(i)
A polygon with “holes” is considered as a set of “rings”: one outer ring defining the boundary and several inner rings describing the “holes”. The outer ring always comes first. Thus, the algorithm will transform as follows: loop over all rings; if it is the first iteration - this is an outer ring, otherwise - inner ring; extract vertices from the ring and output their coordinates.
defdump_geometry(geom):
nRings = geom.GetGeometryCount()
for i in range(nRings):
ring = geom.GetGeometryRef(i)
if i ==0:
print "Exterior ring"else:
print "Interior ring", i
for j in range(ring.GetPointCount()):
print j, ring.GetX(j), ring.GetY(j)
Handling multipart polygons is more complicated: they can be either singlepart polygons or contain “holes”, “islands”, or even “holes” and “islands” at the same time. Let’s try to figure out how to deal with them. A polygon consisting of several “islands” can be considered as a set of independent polygons. In turn, each of those polygons can have one or more rings. So what we have now:
if the polygon is a singlepart polygon and contains only one ring, we extract its vertices immediately
if the number of rings in polygon is greater than 1 — we loop through the rings and extract vertices from them
if the polygon is a multipart polygon — loop through the parts and for each for part perform checks starting from the step 1
So, this is a classical recursion and the final function will look like this:
defdump_geometry(geom, file):
geomType = geom.GetGeometryType() # geometry type# singlepart polygonif geomType == ogr.wkbPolygon or geomType == ogr.wkbPolygon25D:
nRings = geom.GetGeometryCount() # number of rings print "Rings", nRings
if nRings ==1: # there are no "holes" in the polygons print "Exterior ring" ring = geom.GetGeometryRef(0)
for i in range(ring.GetPointCount()):
print i, ring.GetX(i), ring.GetY(i)
else: # polygon with "holes"for i in range(nRings):
ring = geom.GetGeometryRef(i)
if i ==0:
print "Exterior ring"else:
print "Interior ring", i
for j in range(ring.GetPointCount()):
print j, ring.GetX(j), ring.GetY(j)
# multipart polygonelif geomType == ogr.wkbMultiPolygon or geomType == ogr.wkbMultiPolygon25D:
for i in range(geom.GetGeometryCount()):
subGeom = geom.GetGeometryRef(i) # for each polygon part dump_geometry(subGeom) # recursively call dump_geometry() function
The code is a bit bloated, but my aim was to demonstrate a possible way of approaching the task, not to come up with an ideal solution.
GDAL is a free library for working with raster and vector data, OGR is a part of the GDAL and is used to work with vector data. The command line utilities included in the library are widely used to perform a variety of tasks. Thanks to the developed API, you can work with OGR functions from many programming languages. This article is dedicated to using the OGR API in Python and is based on the GDAL Vector API Tutorial.
GDAL is a free library for working with raster and vector data. The command line utilities included in the library are widely used to perform a variety of tasks. Thanks to the developed API, you can work with GDAL functions from many programming languages. This article is dedicated to using the GDAL API in Python and is based on the GDAL Raster API Tutorial.
A serious bug has been found in the recently released GDAL 1.7.0. HFA files created with the new version of GDAL are not compatible with Erdas ViewFinder, ArcGIS, ArcExplorer and other programs. Such files can only be opened with GDAL itself.
The interim bugfix release 1.7.1 is scheduled for Monday. It is not recommended to use version 1.7.0, it will be declared deprecated after the bugfix release.
A fix was committed just two hours ago: trunk (r18728), branch 1.7 (r18729).
