OS OpenData Contours on Garmin

I've added contour lines from OS OpenData Landform Panorama to my Garmin OpenStreetMap mapset. Here's a quick summary of how I did it, based on a Debian system:

DXF to Shapefile

aptitude install libdime-dev libgdal-dev

Compile the following using 'g++ -Wall -g $(gdal-config --cflags) -o dxf2shp dxf2shp.cc $(gdal-config --libs) -ldime -lm -lstdc++':

/**************************************************************************\
 *
 *  Copyright (C) 2010 Toby Speight.  All rights reserved.
 *
 *  Author: Toby Speight <T.M.Speight.90@cantab.net>
 *
 *  This program is free software; you can redistribute it and/or modify it
 *  under the terms of the GNU General Public License, version 2, as
 *  published by the Free Software Foundation.
 *
 *  This library is distributed in the hope that it will be useful, but
 *  WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License (the accompanying file named COPYING) for more
 *  details.
\**************************************************************************/

#include <stdio.h>
#include <stdlib.h>

#include <dime/Input.h>
#include <dime/Model.h>
#include <dime/entities/Entity.h>
#include <dime/entities/Polyline.h>
#include <dime/entities/Vertex.h>

#include <gdal/ogr_core.h>
#include <gdal/ogr_feature.h>
#include <gdal/ogrsf_frmts.h>

static bool parseEntity(const class dimeState * const, class dimeEntity *, void *);

static OGRFieldDefn *idField = new OGRFieldDefn("ID", OFTInteger);
static OGRFieldDefn *typeField = new OGRFieldDefn("TYPE", OFTString);
static OGRFieldDefn *eleField = new OGRFieldDefn("NAME", OFTReal);

int main(int argc, char **argv)
{
  if (argc != 3) {
      fprintf(stderr, "Usage: %s source dest\n", argv[0]);
      return EXIT_FAILURE;
  }

  char *src = argv[1];
  char *dest = argv[2];

  dimeInput in;
  if (!in.setFile(src)) {
      fprintf(stderr, "stdin: open failed\n");
      return EXIT_FAILURE;
  }

  dimeModel model;
  if (!model.read(&in)) {
    fprintf(stderr,"DXF read error in line: %d\n", in.getFilePosition());
    return EXIT_FAILURE;
  }

  const char *driverName = "ESRI Shapefile";
  OGRRegisterAll();
  OGRSFDriverRegistrar *poRegistrar = OGRSFDriverRegistrar::GetRegistrar();
  OGRSFDriver *poDriver = poRegistrar->GetDriverByName(driverName);
  if (!poDriver) {
      fprintf(stderr, "Driver not found for %s\n", driverName);
      fprintf(stderr, "The following drivers are available:\n");
      for (int i = 0;  i < poRegistrar->GetDriverCount();  i++)
          fprintf(stderr, "%d: `%s'\n", i+1, poRegistrar->GetDriver(i)->GetName());
      return EXIT_FAILURE;
  }
  if(!poDriver->TestCapability(ODrCCreateDataSource)) {
      fprintf(stderr,  "%s driver does not support data source creation.\n", driverName);
      return EXIT_FAILURE;
  }

  if (unlink(dest) && errno != ENOENT) {
      perror(dest);
      return EXIT_FAILURE;
  }
  OGRDataSource *source = poDriver->CreateDataSource(dest);
  if (!source) {
      fprintf(stderr, "%s: open failed: %s\n", dest, CPLGetLastErrorMsg());
      return EXIT_FAILURE;
  }
  OGRSpatialReference epsg27700;
  int ogrError = epsg27700.importFromEPSG(27700);
  if (ogrError != OGRERR_NONE) {
      fprintf(stderr, "couldn't initialise EPSG:27700: %d: %s\n", ogrError, CPLGetLastErrorMsg());
      return EXIT_FAILURE;
  }

  OGRLayer *layer = source->CreateLayer("Landform", &epsg27700, wkbLineString);

  idField->SetWidth(8);
  typeField->SetWidth(8);
  eleField->SetWidth(7);
  eleField->SetPrecision(1);

  layer->CreateField(idField);
  layer->CreateField(typeField);
  layer->CreateField(eleField);

  if (!model.traverseEntities(parseEntity, layer, false, false, false))
    fprintf(stderr,"Conversion error\n");

  OGRDataSource::DestroyDataSource(source);


}

static bool parseEntity(const class dimeState * const state, class dimeEntity *entity, void *userData)
{
    static int id = 0;
    OGRLayer *layer = (OGRLayer *)userData;

    switch (entity->typeId()) {
    case dimeBase::dimePolylineType: {
        class dimePolyline *poly = (class dimePolyline *)entity;
        OGRLineString *line = new OGRLineString();
        int vertexCount = poly->getNumCoordVertices();
        for (int i = 0;  i < vertexCount;  i++) {
            dimeVertex *v = poly->getCoordVertex(i);
            dimeVec3f coords = v->getCoords();
            line->addPoint(coords.x, coords.y, coords.z);
        }
        OGRFeature *feature = new OGRFeature(layer->GetLayerDefn());
        feature->SetField(0, id++);
        feature->SetField(1, poly->getLayerName());
        feature->SetField(2, poly->getElevation().z);
        feature->SetGeometryDirectly(line);
        layer->CreateFeature(feature);
        return true;
    }
    default:
        return true;
    }
}

You should now have an executable that you can run on an input DXF file:

<code>./dxf2shp data/sd/sd00.dxf data/sd/sd00.shp</code>

This shapefile can be loaded into Merkaartor etc. to check for correctness.

Shapefile to OSM

A little bit of perl:

aptitude install libgdal-perl

#!/usr/bin/perl

use strict;
use warnings;

use Geo::GDAL;
#use Geo::Proj4;
use XML::Writer;

# Mark the contours with "contour" and "contour_ext" tags
my $med = 50;                   # modulus for medium contours - set to 1 to not have minors
my $large = 250;                # modulus for major contours - set to 100000 to not have majors

my $contour_file = shift;

BEGIN { $SIG{'__WARN__'} = sub {  print STDERR "Warning: @_"; } }

my $wgs84 = new Geo::OSR::SpatialReference();
eval { $wgs84->SetWellKnownGeogCS('WGS84') }; die "GDAL error: $@\n" if $@;
my $epsg27700 = new Geo::OSR::SpatialReference();
eval { $epsg27700->ImportFromEPSG(27700) }; die "GDAL error: $@\n" if $@;

my $transform = new Geo::OSR::CoordinateTransformation($epsg27700, $wgs84);

my $writer = new XML::Writer(DATA_MODE=>1, DATA_INDENT=>1);
$writer->xmlDecl();
$writer->startTag('osm', version=>"0.6");

my $id = 0;
my %types;


my $source = eval { Geo::OGR::DataSource::Open($contour_file) } or die "GDAL error: $@\n";
foreach my $layer_name($source->Layers()) {
    #print STDERR "Reading layer '$layer_name'\n";
    my $layer = $source->GetLayerByName($layer_name);
while (my $feature = $layer->GetNextFeature()) {
    my $name = int $feature->GetField('NAME');
    my $type = $feature->GetField('TYPE');
    my $geometry = $feature->GetGeometry();
    $geometry->Transform($transform);
    my $points;
    foreach (0..$geometry->GetPointCount()-1) {
        my $point = $geometry->GetPoint_2D($_);
        $writer->emptyTag('node', id=>--$id, lon=>$point->[0], lat=>$point->[1]);
        push @$points, $id;
    }
    $types{$type} = [] unless ($types{$type});
    push @{$types{$type}}, $name, $points;
}
}

sub write_ways(&$) {
    my $write_tags = shift;
    my $type = shift;
    my $ways = $types{$type};
    return unless $ways;
    while (@$ways) {
        my $name = shift @$ways;
        my $points = shift @$ways;

        $writer->startTag('way', id=>--$id);
        $write_tags->($writer, $name);
        foreach (reverse @$points) {
            $writer->emptyTag('nd', ref=>$_);
        }
        $writer->endTag('way');
    }

    delete $types{$type};
}

write_ways {
    # Contours
    my $writer = shift;
    my $name = shift;
    my $ext = ($name % $med) ? "minor" : ($name % $large) ? "medium" : "major";
    $writer->emptyTag('tag', k=>'contour', v=>"elevation");
    $writer->emptyTag('tag', k=>'contour_ext', v=>"elevation_$ext");
    $writer->emptyTag('tag', k=>'ele', v=>"$name");
} "G8040201";

write_ways {
    # Lakes
    my $writer = shift;
    $writer->emptyTag('tag', k=>'line', v=>"lake");
} "G8040202";

write_ways {
    # Lakes
    my $writer = shift;
    $writer->emptyTag('tag', k=>'line', v=>"break_line");
} "G8040203";

write_ways {
    # Coastlines
    my $writer = shift;
    $writer->emptyTag('tag', k=>'line', v=>"coast");
} "G8040204";

write_ways {
    # Lakes
    my $writer = shift;
    $writer->emptyTag('tag', k=>'line', v=>"ridge_line");
} "G8040205";

write_ways {
    # Lakes
    my $writer = shift;
    $writer->emptyTag('tag', k=>'line', v=>"form_line");
} "G8040207";

$writer->endTag('osm');
$writer->end();

print STDERR 'No output written for ', join(', ', keys %types), "\n" if (%types) ;

OSM to Garmin

TODO: mkgmap invocation

A Makefile to drive it all

TODO: insert the makefile