More with the GDAL/OGR perl bindings

So my last post talked about my struggles to finally get something saved in the database using the native perl bindings into the GDAL/OGR library. Once I got that working and pushed out the post, I immediately started loading up multiple files and playing around with the data. One thing I noticed was that it was impossible to separate different “trips” within the data without playing around with space and time. What I wanted was an easy way to flag each batch of points with a field identifying the run.

The auto-generated schema for the GPX data looks like this:

d testogr.track_points
                                              Table "testogr.track_points"
       Column       |           Type           |                               Modifiers                                
--------------------+--------------------------+------------------------------------------------------------------------
 ogc_fid            | integer                  | not null default nextval('testogr.track_points_ogc_fid_seq'::regclass)
 wkb_geometry       | geometry(Point,4326)     | 
 track_fid          | integer                  | 
 track_seg_id       | integer                  | 
 track_seg_point_id | integer                  | 
 ele                | double precision         | 
 time               | timestamp with time zone | 
 magvar             | double precision         | 
 geoidheight        | double precision         | 
 name               | character varying        | 
 cmt                | character varying        | 
 desc               | character varying        | 
 src                | character varying        | 
 link1_href         | character varying        | 
 link1_text         | character varying        | 
 link1_type         | character varying        | 
 link2_href         | character varying        | 
 link2_text         | character varying        | 
 link2_type         | character varying        | 
 sym                | character varying        | 
 type               | character varying        | 
 fix                | character varying        | 
 sat                | integer                  | 
 hdop               | double precision         | 
 vdop               | double precision         | 
 pdop               | double precision         | 
 ageofdgpsdata      | double precision         | 
 dgpsid             | integer                  | 
 speed              | double precision         | 
Indexes:
    "track_points_pkey" PRIMARY KEY, btree (ogc_fid)
    "track_points_wkb_geometry_geom_idx" gist (wkb_geometry)

There are three fields that are completely blank: src, desc, and name. I decided to use src to identify the source of the data as the file name it came from.

First I modified my previous program to parse the command line options using Getopt::Long. I don’t use all of its power in this example, but in the past I’ve been well served by starting with that in case the script grows and mutates.

With Getopt::Long, I understand there are ways to input a list of things into the arguments. You can have multiple invocations of the same option, for example, --file mydata.gpx --file moredata.gpx, or you can input them as a comma separated list and follow the recipe in the perldoc for the module. However, I wanted to use a glob, like –file data/*.gpx, so I instead decided to just stick all the files after a double dash on the command line. So really, in the following code, I’m only using Getopt::Long to parse out a –help command! However, it’s there if I need to expand functionality in the future.

use strict;
use warnings;
use Carp;

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

use Getopt::Long;
use Pod::Usage;

my $man = 0;
my $help = 0;

my @files;

my $result = GetOptions(
    'help|?' => $help,
    ) or pod2usage(2);

pod2usage(-exitval => 0, -verbose => 2) if $help;

@files = @ARGV;
...

With that, I have all of my input files in an array, and I can loop over them and store the filename in the source field in the db by using $new_feature->SetField('src',$_);, as follows:

foreach (@files){

    my $ds = Geo::OGR::Open($_);

    my $layer         = $ds->Layer($layer_name);
    my $feature_count = $layer->GetFeatureCount();
    carp "$layer_name, $feature_count";
    if ( $feature_count < 10 ) {
        next;
    }

    carp "saving $_ to pg";

    # 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);

        # write the filename as the src field, for making lines later
        $new_feature->SetField('src',$_);

        my $pgf = $pg_layer->CreateFeature($new_feature);

        $x += 1;
        if ( $x % 128 == 0 ) {
            carp $x;
            # uncomment the following to crash your program
            # $pg_layer->CommitTransaction();
            # StartTransaction() seems to auto commit prior transaction?
            $pg_layer->StartTransaction(); 
            $x = 0;
        }

    }
    if ($x) {
        carp "all done, $x remaining";
        $pg_layer->CommitTransaction(); # this one doesn't crash for some reason
        carp "last transaction committed";
    }
}

That does its magic, and the database now has distinct groups of points. Now if you want to make “lines” out of those points, you can do this in PostGIS:

SELECT ST_MakeLine(wkb_geometry ORDER BY track_seg_point_id ASC) AS linegeom, src
INTO table testogr.lines
FROM testogr.track_points
GROUP BY src;

Et voila

QGIS rendering the new lines table, on top of OSM lines data

QGIS rendering the new lines table, on top of OSM lines data

Of course, that isn’t at all helpful, as I want to see speeds, not just the lines. Next step is to try to figure out how to add a measure to each point, and then collect those (X,Y,M) type points into a line with a measure dimension. I guess that will be my next post.

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;

namespace::clean compile problem solved

I couldn’t get namespace::clean to compile and install properly on one of my servers. It was a clean Perl 5.10.1 install, so I couldn’t really see what the problem was. But just to be safe I really really cleaned out perl, reinstalled, and then ran Gentoo’s perl cleaner utility, ran a revdep-rebuild, etc etc. Then I found this bug, installed Term::ReadLine::Gnu, and all was well. So simple, and yet not something I’d ever figure out on my own. It would be pretty cool if the CPAN testers who failed on this test and those who passed could compare notes automatically and generate a list of possible dependencies that are missing.

But anyway, I really hate it when CPAN installs fail and I can’t figure out why right away. Maybe Google will pick up this blog and index it. Here Google, index this: Compilation failed in require at /usr/lib64/perl5/site_perl/5.10.1/Term/ReadLine/Perl.pm line 63. at /usr/lib64/perl5/site_perl/5.10.1/Term/ReadLine/Perl.pm line 63 Term::ReadLine::Perl::new(‘Term::ReadLine’, ‘perldb’, ‘GLOB(0x1caace8)’, ‘GLOB(0x1b9a150)

A simple “show me” program to understand what lazy_build, etc does

I’m beginning to use MooseX::Declare more and more, but this morning I realized I didn’t quite understand when the builder was getting called, and under what circumstances, so I wrote the following program to testing things out. Not so much a test of the code for some notion of correctness as a literal test of what is going on so I can be more informed.
Continue reading

Yikes!

Okay, two lessons. No, three. First iron man perl planet rocks. Second, from now on push up code examples with my code problems. Third, I need to turn off moderation on this blog, since the spam filter seems to work pretty well.

I’m working up an example of the problem I had with DateTime and my attempts to solve it, and I’ll fold in the suggestions I got on my previous post (especially if they solve the problem!)

I used perl today, and I can’t figure out how to get my paper man icon.

I used perl today. Moose. I ran into a problem. It was annoying. I am way too stressed and tired to blog more. But I will anyway, secure in the knowledge that no one reads this blog but google’s spiders.

Okay anyway I used MooseX::Declare, and couldn’t get the method signature stuff to work. I did something like

method weekend_or_vacation (DateTime $dt){
  # check if weekend or vacation
  # with vacation being the tricky bit
  if($vacation || $weekend){
   return 1;
  }else{
    return 0;
  }
}

But MooseX::Declare kept complaining that it didn’t know what DateTime was. I scanned the tests in t and sure enough, they all test simple things like Str and ArrayRef and so on, but none of the more magical parts of type checking.

I eventually solved it the old fashioned way by puking if the argument wasn’t a DateTime, but I’d rather do it the method signature way.

PUT problem solved

I had a problem linking up dojo/xhrPut and Catalyst::Controller::REST. As always, the answer was in the documentation, but I didn’t see it.

Catalyst::Controller::REST docs say that:

The HTTP POST, PUT, and OPTIONS methods will all automatically deserialize the contents of $c->request->body based on the requests content-type header. A list of understood serialization formats is below.

And the docs for dojo/xhrPut point to those for dojo/xhrGet for parameters, which include:

headers

A JavaScript object of name/string value pairs. These are the headers to send as part of the request. For example, you can use the headers option to set the Content-Type, X-Method-Override, or Content-Encoding headers of the HTTP request.

This parameter is optional

So all I had to do in my javascript code is

	    var xhrArgs = {
		url: ajaxurls.sort + '/' + editing.id,
		putData: dojo.toJson(data),
		handleAs: "json",
		headers: {'Content-Type':'application/json'},
		load: function(data){
		    // don't really need to do anything here
		    // uncomment for testing
                    // console.log("new sort order put to server");
		},
		error: function(error){
		    alert('warning, edits were not saved properly.  Proceed with caution');
		}
	    }
	    //Call the asynchronous xhrPost
	    var deferred = dojo.xhrPut(xhrArgs);
	

And the controller magically started to work as expected. Hooray, git commit and all that, but it is time to go to sleep and actually get things working some other day.

A Flash-related browser crash ate my bug report to MooseX::Declare

ggghhhhaaaa. I hate flash. I really like http://proquest.safaribooksonline.com/, or rather, I used to love it, but now they’ve switched to Flash and it is hateful hateful hateful. But I still can’t stop using it because the information is so awesome and handy and because UCI has an account and it is right there waiting for me whenever I have a question. But then *bang* one page too many and firefox just blinks off my desktop.

Continue reading