23.10.2011 12:57 · GIS · qgis, gdal, howto
Quite often, I have come across the question: How can I use the georeferenced rasters (raster + map file) from OziExplorer in QGIS? In order to use such data in a GIS, it is usually necessary to convert them first, e.g., with GlobalMapper. With the release of GDAL 1.7.1, there is no need for proprietary software, as GDAL now supports OZI map files.
Let’s see how to do this in QGIS.
Read more ››GDAL is gradually moving towards full Unicode support: RFC 23 has already been implemented, and there are ongoing efforts to implement RFC 5.
Another step was the implementation of recoding of Shapefile attributes to UTF-8 when reading and from UTF-8 when writing. Except for that… the encoding is determined by reading the LDID (Language Driver ID) from the DBF header. In general, this is the right approach, but I can’t remember the last time I saw shapfiles with the encoding specified correctly. Mostly, there are files with LDID set to 87, which corresponds to the default value.
This is where the most interesting part begins. It is clear that this default is different for everyone. And in the current implementation, the LDID/87 value is interpreted as ISO8859_1 (Latin-1). The problem with this approach, I think, is clear to everyone. The proposed solution is either to edit the existing files and set the required DBF encoding or to override the interpretation of the LDID value by setting an environment variable. The first method requires more work (it is necessary to find out the encoding of each Shape file and write it in the DBF header), but it is also the most correct. The second one is actually an ugly workaround because, with this approach, only files in one, overridden, encoding will be read correctly. All others will still be displayed as unreadable garbage, which is unacceptable when using data in different encodings.
A long time ago, in a galaxy far, far away… Almost a year ago, I was involved in a project that is now almost dead. It was notable not only for its ideas, but also for the fact that its main development tool was Delphi. In my opinion, this was not the best choice: if you need ObjectPascal so much, you should use FreePascal, which supports 64-bit and is cross-platform, open, and free. But that’s not the point.
The GDAL library was supposed to be used to work with vector and raster data, but… As it turned out, there are no bindings for Pascal. There are bindings for Python, R, Perl, PHP and a few others, but not for Pascal. Well, there is a ticket in trac from 4 years ago, but the solution proposed there is basically a hack, because it requires changes in the GDAL code and a custom GDAL build. We had to do something… So I remembered my youth, remembered how I used to code in Pascal and tried to make a wrapper for the library.
I tried to keep the code compatible with both FreePascal and Delphi. Unfortunately, for various reasons, I did not manage to create a complete wrapper. Maybe when I have the time and resources, I will slowly add the missing bits. In the meantime, I have decided to publish the current implementation on GitHub. I hope someone will find it useful.
27.07.2011 15:41 · GIS · python, gdal, howto
It is no secret that one of the most common raster formats, TIFF, is able to store georeferencing information within the file, turning it into a GeoTIFF. Everything would be fine if it weren’t for one nuance: the vast majority of graphic editors suffer from “star fever”, believing that unrecognised tags have no place in the file and happily removing this valuable information when saving the file.
This feature of graphic editors is often used to “reset” georeferencing (for example, if the image is incorrectly georeferenced and needs to be georeferenced from scratch). A natural question arises: what to do if you need to edit the image in a graphics editor but do not want to lose the georeferencing? The answer is not original, it is necessary to save it in an external file (the so-called world file).
Let’s consider three possible solutions.
Read more ››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.
And this is the result of the same operation performed by gdal_calc
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
myval=BandReadAsArray(myFiles[i].GetRasterBand(myBands[i]), xoff=myX, yoff=myY, win_xsize=nXValid, win_ysize=nYValid)
change to this
myval=BandReadAsArray(myFiles[i].GetRasterBand(myBands[i]), xoff=myX, yoff=myY, win_xsize=nXValid, win_ysize=nYValid).astype(numpy.float32)
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 my 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 rasters
gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="A+B"
# average of two rasters
gdal_calc.py -A input1.tif -B input2.tif --outfile=result.tif --calc="(A+B)/2"
# difference of raster bands
gdal_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.
29.01.2011 09:21 · GIS · qgis, gdal, osgeo4w
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.
Follow up to the GDAL RFC 29 post.
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.
12.10.2010 18:01 · GIS · gdal, python, tips
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 their coordinates. The code will look like this:
def dump_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 polygon 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.
def dump_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:
So, this is a classical recursion and the final function will look like this:
def dump_geometry(geom, file):
geomType = geom.GetGeometryType() # geometry type
# singlepart polygon
if 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 polygon
elif 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.