Month: April 2012

  • geocoder.ca’s excellent database

    Update, late 2017: Looks like geocoder.ca’s database hasn’t been updated since late 2015. The download link is still active, though.

    geocoder.ca provides a crowdsourced postal code geocoder under the ODbL. You can download the database as CSV directly. Here’s a bash script to convert that text file into a (very large) point shapefile:

    #!/bin/bash
    # geocoder2shp.sh - convert geocoder.ca CSV to a shape file
    # NB: input CSV is UTF-8; it is passed through unchanged
    # Needs >= v1.7 of GDAL
    # scruss - 2012/04/15
    
    if [ ! -f Canada.csv.gz ]
    then
        echo ""
        echo " " Download \"Canada.csv.gz\" into the current directory from
        echo "  " http://geocoder.ca/onetimedownload/Canada.csv.gz
        echo " " and try again.
        echo ""
        exit 1
    fi
    
    # make input file with header
    echo PostalCode,Latitude,Longitude,City,Province > Canada2.csv
    gunzip -c Canada.csv.gz >> Canada2.csv
    
    # create GDAL VRT file
    cat > Canada2.vrt <<EOF
    <OGRVRTDataSource>
      <!-- note that OGRVRTLayer name must be basename of source file -->
      <OGRVRTLayer name="Canada2">
        <SrcDataSource>Canada2.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <LayerSRS>EPSG:4326</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="Longitude" y="Latitude"/>
        <Field name="PostalCode" type="String" width="6" />
        <Field name="Latitude" type="Real" />
        <Field name="Longitude" type="Real" />
        <Field name="City" type="String" width="60" />
        <Field name="Province" type="String" width="2" />
      </OGRVRTLayer>
    </OGRVRTDataSource>
    EOF
    
    # create shapefile
    ogr2ogr PostalCodes.shp Canada2.vrt
    
    # clean up
    rm -f Canada2.csv	Canada2.vrt
    

    Though the script is a bit Unix-centric, it’s just a simple list of instructions which could be run on any command line. What it does is add some headers to the geocoder.ca file, then sets up an OGR Virtual Format to convert the text into a fairly well-defined shapefile. When you use this shapefile, you should credit geocoder.ca as the ODbL requires.

    Eek! geocoder.ca has been sued by Canada Post! (News responses: Michael Geist, Boing Boing, CBC) I’ve donated to defend this useful service.

  • 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 =&amp;gt; 25.4;
    
    my $qgis = $ARGV[0];
    die &amp;quot;$qgis must exist\n&amp;quot; unless ( -f $qgis );
    
    my $q = XMLin($qgis) or die &amp;quot;$!\n&amp;quot;;
    my %composer = %{ $q-&amp;gt;{Composer} };
    my $image_width =
      int( $composer{Composition}-&amp;gt;{paperWidth} *
        $composer{Composition}-&amp;gt;{printResolution} /
        MM_PER_INCH );
    my $image_height =
      int( $composer{Composition}-&amp;gt;{paperHeight} *
        $composer{Composition}-&amp;gt;{printResolution} /
        MM_PER_INCH );
    
    # we need xpixelsize, ypixelsize, ulx and uly
    my $xpixelsize =
      ( $composer{ComposerMap}-&amp;gt;{Extent}-&amp;gt;{xmax} -
        $composer{ComposerMap}-&amp;gt;{Extent}-&amp;gt;{xmin} ) /
      int( $composer{ComposerMap}-&amp;gt;{ComposerItem}-&amp;gt;{width} *
        $composer{Composition}-&amp;gt;{printResolution} /
        MM_PER_INCH );
    my $ypixelsize =
      -1.0 *
      ( $composer{ComposerMap}-&amp;gt;{Extent}-&amp;gt;{ymax} -
        $composer{ComposerMap}-&amp;gt;{Extent}-&amp;gt;{ymin} ) /
      int( $composer{ComposerMap}-&amp;gt;{ComposerItem}-&amp;gt;{height} *
        $composer{Composition}-&amp;gt;{printResolution} /
        MM_PER_INCH );
    my $ulx =
      $composer{ComposerMap}-&amp;gt;{Extent}-&amp;gt;{xmin} -
      $xpixelsize *
      int( $composer{ComposerMap}-&amp;gt;{ComposerItem}-&amp;gt;{x} *
        $composer{Composition}-&amp;gt;{printResolution} /
        MM_PER_INCH ) - $xpixelsize;
    my $uly =
      $composer{ComposerMap}-&amp;gt;{Extent}-&amp;gt;{ymax} -
      $ypixelsize *
      int( $composer{ComposerMap}-&amp;gt;{ComposerItem}-&amp;gt;{y} *
        $composer{Composition}-&amp;gt;{printResolution} /
        MM_PER_INCH ) - $ypixelsize;
    
    printf( &amp;quot;%.12f\n%.12f\n%.12f\n%.12f\n%.12f\n%.12f\n&amp;quot;,
      $xpixelsize, 0.0, 0.0, $ypixelsize, $ulx, $uly );
    
    # FIXME? pixel scale seems a tiny bit off - allow for border?
    exit;