Categories
GIS

Georeferencing QGIS output rasters

Some of the design packages I rely on use very crude GIS facilities. In fact, all they can support is a georeferenced raster as a background image, so it’s more of a rough map than GIS. It helps if these rasters are at a decent resolution, typically 1m/pixel or better.

A while back, I asked on the QGIS forum if the package could output high resolution georeferenced rasters. I received a rather terse response that no, it couldn’t (and I inferred from the tone that the poster thought that it shouldn’t, and I was wrong to want such a thing). I shelved the idea at the time.

After having to fix a lot of paths in a QGIS project file I’d moved to a new system, I noticed that all the map composer attributes are rather neatly defined in the XML file structure. Some messing around with Perl, XML::Simple and Data::Dumper::Simple, and I had a little script that would spit out an ESRI World File for the map composer raster defined in the project.

To run this, you have to create a project with just one Print Composer page, save the composed map as an image, save the project, then run the script like this:

./geoprint.pl project.qgs > image.pngw

There are some caveats:

  • This probably won’t work for projects with multiple print composers
  • It doesn’t quite get the scale right, but it’s within a pixel or so. I may not have corrected for image borders.

Though there’s some fairly hideous XML-mungeing in the code, what the script does is entirely trivial. If you feel you can use it, good; if you feel you can improve it, be my guest.

#!/usr/bin/perl -w
# geoprint - georef a QGIS output image by creating a world file
# one arg: qgis project file (xml)
# $Id: geoprint.pl,v 1.3 2012/04/06 03:32:01 scruss Exp $

use strict;
use XML::Simple;
use constant MM_PER_INCH => 25.4;

my $qgis = $ARGV[0];
die "$qgis must exist\n" unless ( -f $qgis );

my $q = XMLin($qgis) or die "$!\n";
my %composer = %{ $q->{Composer} };
my $image_width =
  int( $composer{Composition}->{paperWidth} *
    $composer{Composition}->{printResolution} /
    MM_PER_INCH );
my $image_height =
  int( $composer{Composition}->{paperHeight} *
    $composer{Composition}->{printResolution} /
    MM_PER_INCH );

# we need xpixelsize, ypixelsize, ulx and uly
my $xpixelsize =
  ( $composer{ComposerMap}->{Extent}->{xmax} -
    $composer{ComposerMap}->{Extent}->{xmin} ) /
  int( $composer{ComposerMap}->{ComposerItem}->{width} *
    $composer{Composition}->{printResolution} /
    MM_PER_INCH );
my $ypixelsize =
  -1.0 *
  ( $composer{ComposerMap}->{Extent}->{ymax} -
    $composer{ComposerMap}->{Extent}->{ymin} ) /
  int( $composer{ComposerMap}->{ComposerItem}->{height} *
    $composer{Composition}->{printResolution} /
    MM_PER_INCH );
my $ulx =
  $composer{ComposerMap}->{Extent}->{xmin} -
  $xpixelsize *
  int( $composer{ComposerMap}->{ComposerItem}->{x} *
    $composer{Composition}->{printResolution} /
    MM_PER_INCH ) - $xpixelsize;
my $uly =
  $composer{ComposerMap}->{Extent}->{ymax} -
  $ypixelsize *
  int( $composer{ComposerMap}->{ComposerItem}->{y} *
    $composer{Composition}->{printResolution} /
    MM_PER_INCH ) - $ypixelsize;

printf( "%.12f\n%.12f\n%.12f\n%.12f\n%.12f\n%.12f\n",
  $xpixelsize, 0.0, 0.0, $ypixelsize, $ulx, $uly );

# FIXME? pixel scale seems a tiny bit off - allow for border?
exit;
Categories
GIS

image georeferencing with QGIS

Quantum GIS (QGIS) has a very powerful image georeferencing module. What that allows you to do is convert screen pixels to a map of an area. The pixels could be from a scanned map or from a screen image. Scanned maps can sometimes be a little distorted, but QGIS’s georeferencer can handle that, within limits.

For an example, say I live in Redickville, ON (I don’t, but some folks do). I’ve heard from North Dufferin Agricultural & Community Taskforce (NDACT) that a large quarry is planned in my area. How close am I going to be from it?

NDACT has a helpful map (2.4 MB PNG image) had a helpful map (which I’ve kept here so you can work through this: MelancthonAerial), but it has no scale. It’s clearly derived from something like Google Maps, so it’s not exactly a free image I can throw about. One thing about georeferencing is that both the source image and the map from which you take you control points affect the licensing of the final map. You’re ending up with a derived work.

There are lots of sources of coordinates for your control points. You could always use Google Maps, but then you’re well down the derived work sinkhole. GeoGratis has a bunch of good data sources, and they are free to use. I’m going to use Toporama‘s digital image maps, as they’re clear and fairly accurate.

Redickville is on Toporama sheet 041A01. I’ve downloaded it as UTM, as it’s easier to measure distances that way. You want to set up your project so it uses the coordinate reference system (CRS) that the control point map uses. Toporama uses EPSG 26917 in the area (easily checked with gdalinfo; it’ll come up with something like AUTHORITY[“EPSG”,”26917″]]), so you should set the project to use that CRS:

And here’s the raster map loaded into QGIS;

Opening up the georeferencer plugin (which is now in the Raster menu, as of QGIS 1.8) gives you a whole lot of blank:

If we open the raster, you can zoom into the area into which you want to put control points. You want to have the highest zoom that the map’s still clear, as the accuracy of your final map depends on how well you placed control points. I’ve chosen road intersections as my control points. Here are the ones I’m going to use:

Point   E-W Road        N-S Road        Note
1       County Rd 21    2nd Line W      Honeywood
2       County Rd 21    4th Line       
3       20th Side Rd    5th Line       
4       15th Side Rd    County Rd 124   just S of where 124 straightens
5       15th Side Rd    Mulmur Townline
6       4th Line        5th Line        4th Line really runs SE-NW

These intersections are fairly well spaced apart, and are clear on both maps. So I choose a pixel on the map in the georeferencer, and a dialogue comes up:

(If you’ve downloaded the MelancthonAerial archive from here, the georeferencing points should load automatically from the “Melancthon Aerial July 22 09.png.points” file. If you want to go through the exercise of manually adding reference points, delete the points file.)

Here I’ve already clicked on the corresponding point on the Toporama map, and the coordinates have been filled in. If you know the coordinates, you can enter them in the boxes. Once you click OK, you have your first control point:

We can add the rest one by one. QGIS’s “Zoom Previous” is really useful for flipping between scales on both the plugin window and the main map. It’s probably a good idea to “Save GCP Points As …” every now and again so you don’t lose your work. Here are all of the points in the georeferencer plugin:

Now you want to modify the Transformation Settings; it’s the little wrench/spanner in the toolbar:

There are lots of options here:

  • Transformation Type: Linear is useful if you’re just adding georeference data to an already computer generated map. Helmert is a simple shear/scale/rotate transformation. The various Polynomial types will correct more gross distortions, but need many control points and can be processor intensive. Thin Plate Spline will distort your map locally to match GCPs; this can work well if your map’s a bit “vernacular”, but if your control points are wrong, your map will end up hilariously melty.
  • Resampling Method: This controls how the output pixels are calculated. Nearest Neighbour is quick and blocky, but useful if you’re mapping an image that has sharp transitions. Cubic maintains more detail, while applying some smoothing at the cost of some detail loss and a fair bit of processing power. The other options can look nice, but eat CPU. This is worth experimenting with many options, as there’s no one solution for all maps.
  • Compression: This controls the file size of the output GeoTiff file. JPEG and Deflate can result in small files, but there’s a chance that other GIS systems can’t read the data. Note that JPEG is lossy, and will lose some detail.

Don’t forget to set the output raster file, and make sure that the target SRS matches the CRS (EPSG 26917) you chose earlier.

Helmert is a useful transformation type to check your work. The plugin plots the residual difference between the two sets of map control points as red lines. Here, I’ve clearly made an error in point ID 2:

If you unselect the bad point, the plugin quickly calculates the residuals again. When this one bad point is removed, the residuals drop down to less than 1.0 for each point. I’ve got enough points for a Helmert transformation with 5 GCPs, but I’d probably want some more points for more complex transformations.

Once you hit “Start Transformation”, QGIS will create your referenced image. If you’ve chosen a simple linear transformation, it won’t create a GeoTiff file, but just a world file for your image.

So here’s the georeferenced image overlaid on the Toporama map, with a bit of transparency. It’s quite a good match:

(the above’s an image saved straight from QGIS. It helpfully creates a world file too, so here it is: redickville.jpgw)

So to answer the hypothetical question, Redickville’s pretty well surrounded by this proposed quarry.

Categories
GIS

MNR DEM tiles

So I got the DEM tiles for my zone, couriered from Peterborough. Even zipped, they are several gigabytes. Here’s a visual guide to the tiles I got:

It seems that every project I work on is split across several tiles. My fair city is no exception:

Most of the city is in tile 87, but the northern and eastern corners poke into tiles 91 & 92. Merging the files and pulling out the area isn’t that hard.

First off, we need to know the extents we want. Conveniently, I’d already transformed the Open Toronto data to the same CRS, so I could just go from the extents of one of the city shape files:

$ ogrinfo -al ../open-toronto-shp/dest/2958/city-wards/TCL3_ICITW.shp | grep '^Extent'
Extent: (609538.195527, 4826360.800158) - (651613.577420, 4857439.389359)

You can then use gdal_merge.py to merge the tiles. Annoyingly, rather than using the familiar ‘lower left – upper right’ coordinate convention for bounding box, it uses ‘upper left – lower right’, so working out a buffered bounding box is a bit of a pain:

gdal_merge.py -v -of EHdr -o Toronto.flt -n -9999 -ul_lr 609000 4858000 652000 4826000 `find . -iname 'd*.flt'`

The MNR DEMs are shipped as ESRI .hdr Labelled files, 32-bit floating point. The CRS used is NAD83 Zone 17N. For all that they’re large, they’re quite quick to load and move around.