Using GDAL/OGR perl bindings to load GPX files into PostgreSQL/PostGIS

Today I wrote a short perl program to import GPX files into PostgreSQL using the OGR library’s native perl bindings. This was a super pain to figure out because the naive way doesn’t work, and it appears all the documentation pushed out to mailing lists and on various wikis talks about Python.

OGR has an excellent tool called ogr2ogr that allows you to append data. However, I didn’t want to use that because I wanted to fiddle with the data first, the pipe it to SQL. Specifically, I wanted to delete long pauses at stop lights, etc., and I wanted to use some logic to make sure I didn’t blindly reload old GPX files.

My initial solution was to simply copy the GPX layer in, and then hunt around for a way to flip on an “append” option. My initial program looked like:

use strict;
use warnings;
use Carp;

use Geo::GDAL;
use Data::Dumper;

# Establish a connection to a PostGIS database
my $pg = Geo::OGR::GetDriverByName('PostgreSQL');
if ( !$pg ) {
    croak 'PostgreSQL driver not available';
}

my $conn = $pg->Open( "PG:dbname='osm' user='james' schemas=testogr", 1 );

if ( !$conn ) {
    croak 'choked making connection';
}

my $ds = Geo::OGR::Open('../test/2014-07-10_07-29-12.gpx');

my $pg_layer;
my $defn;

## I'm only interested in the track_points layer
my $layer_name = 'track_points';
my $layer      = $ds->Layer($layer_name);

# use copy
$pg_layer = $conn->CopyLayer( $layer, $layer_name, { 'overwrite' => 1 } );
if ( !$pg_layer ) {
    carp 'failed to copy';
}

1;

That works, but curiously the automatic FID doesn’t automatically increment when using CopyLayer. No matter, I don’t actually use that, because I like creating my own table definitions.

And even if that did work properly, it would only work once. Every other time, that “overwrite” option on the CopyLayer command is going to wipe the table.

Poring over the docs, I didn’t see any option for “append” as was used in the ogr2ogr utility. So I combed through the ogr2ogr source code, and discovered that the “-append” option actually causes the code to create each feature and add it to the existing layer inside of a loop by iterating over the each of the fields in the layer:

    if (papszFieldMap && bAppend)
    {
        int bIdentity = FALSE;

        if (EQUAL(papszFieldMap[0], "identity"))
            bIdentity = TRUE;
        else if (CSLCount(papszFieldMap) != nSrcFieldCount)
        {
            fprintf( stderr, "Field map should contain the value 'identity' or "
                    "the same number of integer values as the source field count.n");
            VSIFree(panMap);
            return NULL;
        }

        for( iField=0; iField < nSrcFieldCount; iField++)
        {
            panMap[iField] = bIdentity? iField : atoi(papszFieldMap[iField]);
            if (panMap[iField] >= poDstFDefn->GetFieldCount())
            {
                fprintf( stderr, "Invalid destination field index %d.n", panMap[iField]);
                VSIFree(panMap);
                return NULL;
            }
        }
    }

So I tried something like that, but for some reason I kept failing to be able to add the new feature to the existing PostgreSQL layer. My broken code looked like:

if ( !$append ) {
    $pg_layer = $conn->CopyLayer( $layer, $layer_name );
    if ( !$pg_layer ) {
        carp 'failed to copy';
    }
}
else {
    if ( !$pg_layer ) {

        # try to get the layer from db
        $pg_layer = $conn->GetLayerByName($layer_name);
        $defn     = $pg_layer->GetLayerDefn();
    }

    # now append each feature
    while ( my $feature = $layer->GetNextFeature() ) {

        my $newFeature = Geo::OGR::Feature->new($defn);

        # Add field values from input Layer
        for my $fi ( 0 .. $defn->GetFieldCount() - 1 ) {
            $newFeature->SetField( $defn->GetFieldDefn($fi)->GetNameRef(),
                $feature->GetField($fi) );

            # Set geometry
            $newFeature->SetGeometry( $feature->GetGeometryRef() );
        }

        # THIS BREAKS 
        my $pgf = $pg_layer->InsertFeature($newFeature);

    }
}

And many variations on that theme, including just trying to directly copy in the feature with $pg_layer->InsertFeature($feature).

The unhelpful error read:

RuntimeError Illegal field type value at /usr/local/lib64/perl5/Geo/OGR.pm line 1473.

I hacked out a little instrumentation around Geo/OGR.pm line 1473, but then I found out that the problem “field type value” changed every time, which made me think I was doing something wrong.

Finally, after giving up twice, I stumbled on an old mailing list posting here. Again, it was in Python, but I read Python well enough to translate into perl without problems. With a little bit of hacking around a buggy call to CommitTransaction(), it worked! My final code looks like:

use strict;
use warnings;
use Carp;

use Geo::GDAL;
use Data::Dumper;

# Establish a connection to a PostGIS database
my $pg = Geo::OGR::GetDriverByName('PostgreSQL');
if ( !$pg ) {
    croak 'PostgreSQL driver not available';
}

my $conn = $pg->Open( "PG:dbname='osm' user='james' schemas=testogr", 1 );

if ( !$conn ) {
    croak 'choked making connection';
}

my $ds = Geo::OGR::Open('../test/2014-07-14_17-56-45.gpx');

my $pg_layer;
my $defn;
my $layer_name = 'track_points';

my $layer         = $ds->Layer($layer_name);
my $feature_count = $layer->GetFeatureCount();
carp "$layer_name, $feature_count";
if ( $feature_count < 10 ) {
    croak;
}
carp "saving to pg";
if ( !$pg_layer ) {

    # try to get the layer from db
    $pg_layer = $conn->GetLayerByName( $layer_name, 1 );
    $defn = $pg_layer->GetLayerDefn();
    carp $pg_layer->GetFeatureCount();
}

# now append each feature
my $x = 0;
$pg_layer->StartTransaction();
while ( my $feature = $layer->GetNextFeature() ) {

    my $new_feature = Geo::OGR::Feature->new($defn);
    $new_feature->SetFrom($feature);
    my $pgf = $pg_layer->CreateFeature($new_feature);

    $x += 1;
    if ( $x % 128 == 0 ) {
        carp $x;
        # leaving this uncommented causes a crash.  Bug?
        # $pg_layer->CommitTransaction();
        $pg_layer->StartTransaction();
        $x = 0;
    }

}
if ($x) {
    carp "all done, $x remaining";
    # curiously, this call to CommitTransaction works okay
    $pg_layer->CommitTransaction();
    carp "last transaction committed";
}
1;

Advertisements