[Openstreetmap] london free maps...

Mon Nov 15 00:34:41 GMT 2004

[Matt Amos]
> any set of algorithms we put together is going to have difficiencies, 
> so i think it might be good to have Hough as stage 1 out of N of the 
> automatic processes?

Yes.  It might work fairly well for cities. :)

> i know what you mean, but i think lots of innacurate info is probably 
> as good as lots of accurate info. its just that little accurate info 
> is much better than little innaccurate info.

Well, my current collection script include the number of satellites
and the horizontal dilution info, so I am starting to collect more
points.  I include the script at the end of this email.

With this script, I'm able to collect points, export them as GPX and
import the GPX file into qgis.  This give me a nice view of the
collected data. :)

> from what i've heard GPSX doesn't allow for extra fields, such as DOP, 
> PR, ephemeris to be included and i really think that we want these 
> fields, if at all possible, to increase our accuracy.

It includes dilution of precision and a number of other fields.  Look
for hdop and vdop.  I have no idea what you mean with PR and
ephemeris, so I suspect those are missing from GPX.

> what i'm proposing is that we extend GPSX to have optional fields
> for these and (maybe) give it a full XML DTD, then contribute it
> back to GPSX so that it is "standard" in some way. disclaimer: i
> have done very little work with XML before and might be proposing
> something very difficult.

This might be a good idea, if the info is actually missing.

#!/usr/bin/perl
#
# Author:  Petter Reinholdtsen
# Date:    2004-10-07
# License: GPL
#
# This script serves two purposes:
#  - Connect to gpsd, and store coordinate info into an SQL database.
#  - Fetch data from the SQL database, export it as GPX (-g)

use strict;
use warnings;
use DBI;
use Socket;
use IO::Handle;
use Date::Manip;
use Getopt::Std;

use vars qw($debug $gpsdhost $gpsdport $gps_stream %opts);

getopt("dh:sgt:", \%opts);

$debug = $opts{'d'} || 0;

$gpsdhost = $opts{'h'} || "localhost";
$gpsdport = "2947";

$gps_stream = $opts{s} || 1;

my $period = 1; # How many seconds between each position check?

# How many seconds delay between points result in a new track? (gpx dump)
my $timelimit = $opts{'t'} || (2 * $period);

# How long a jump is too long?  When should we make a new track?
my $jumplimit = 0.005;

# Connect to database

my $dbname  |= "pere";
my $dbuser  |= "" ; # "pere";
my $dbpass  |= "" ; # "pere";
my $dbtable |= "tempPoints";

my $dbh = DBI->connect("dbi:Pg:dbname=$dbname", $dbuser, $dbpass,
		       {'AutoCommit' => 1}) ||
    die "Unable to connect to database";

if ($opts{'g'}) { # Generate GPX output
    generate_gpx();
} else {

    # Connect to gpsd
    my $socket = tcp_connect($gpsdhost, $gpsdport);

    # Turn off buffering on stdout
    $| = 1;

    print "Starting $0\n";

    print $socket "RMB\n" if $gps_stream;

    while (1) {
	print $socket "APS\n" unless $gps_stream;

	# fetch next line from socket
	my $line = <$socket>;

	if ( ! defined $line) { # Reconnect to gpsd.  It probably restarted
	    $socket->close();
	    sleep 1;
	    $socket = tcp_connect($gpsdhost, $gpsdport);
	    print $socket "RMB\n" if $gps_stream;
	    next;
	}
	$line =~ s/\r\n//;
	$line =~ s/\r//;
	$line =~ s/\n//;
	print "\nL: '$line'\n" if $debug > 2;

	# Each sentence starts with a "$", a two letter "talker ID", a
	# three letter "sentence ID", followed by a number of data
	# fields separated by commas, and terminated by an optional
	# checksum, and a carriage return/line feed.  A sentence may
	# contain up to 82 characters including the "$" and CR/LF.

	my ($tid, $data) = $line =~ m/^[\$](..)(.+)$/;
	$line = $data if ($tid);

	print "L1: '$line'\n" if $debug > 2;
	if ($line =~ m/GPSD,A=.+,P=.+,S=.+/) {
	    print "APS match\n" if $debug;
	    my ($lat, $lon, $altitude, $quality);
	    my @token = split(/,/, $line);
	    for (@token) {
		if (m/^P=(\S+) (\S+)/) {
		    $lat = $1;
		    $lon = $2;
		}
		$altitude = $1 if (m/^A=(\S+)/);
		$quality = $1 if (m/^S=(\d+)/);
	    }
	    add_point(
		      'lat' => $lat,
		      'lon' => $lon,
		      'altitude' => $altitude,
		      'quality' => $quality
		      );
	    sleep $period;
	} elsif ($line =~ /^GGA,/) {
# GGA,105434.016,5957.0233,N,01047.0231,E,1,04,05.1,00143.6,M,38.4,M,,*6D'

# GA - Global Positioning System Fix Data
#         GGA,123519,4807.038,N,01131.324,E,1,08,0.9,545.4,M,46.9,M, , *42
#            123519       Fix taken at 12:35:19 UTC
#            4807.038,N   Latitude 48 deg 07.038' N
#            01131.324,E  Longitude 11 deg 31.324' E
#            1            Fix quality: 0 = invalid
#                                      1 = GPS fix
#                                      2 = DGPS fix
#            08           Number of satellites being tracked
#            0.9          Horizontal dilution of position
#            545.4,M      Altitude, Metres, above mean sea level
#            46.9,M       Height of geoid (mean sea level) above WGS84
#                         ellipsoid
#            (empty field) time in seconds since last DGPS update
#            (empty field) DGPS station ID number
	my (undef,$when, $lat, $latdir, $lon, $londir,
	    $quality, $satellites, $hdilution,
	    $altitude, undef, $height, undef,
	    $dgpswhen, $dgpsid) = split(/,/, $line);
	$lat = convert_degree(2, $lat);
	$lon = convert_degree(3, $lon);
	add_point(
		  'lat' => $lat,
		  'lon' => $lon,
		  'altitude' => $altitude,
		  'hdilution' => $hdilution,
		  'satellites' => $satellites,
		  'quality' => $quality
		  );
    } else {
	print "Unrecognized format: '$line'\n" if $debug > 2;
    }
}
}
$dbh->disconnect;

sub convert_degree {
    my ($leading, $value) = @_;
    my ($deg, $min) = $value =~ m/^(.{$leading})(.+)$/;
    return $deg + $min/60.0;
}

sub add_point {
    my %args = @_;
    my $lat = $args{lat};
    my $lon = $args{lon};
    my $g = "'($lat, $lon)'::point";
    my %data;
    $data{g} = $g;
    for my $field (qw(altitude quality hdilution satellites)) {
	$data{$field} = $dbh->quote($args{$field}) if (exists $args{$field});
    }
    if (0 != $args{quality}) { # Satellites in view
	# This SQL is PostgreSQL specific
	my @fields;
	my @values;
	for my $field (sort keys %data) {
	    push @fields, $field;
	    push @values, $data{$field};
	}
	my $sql = "INSERT INTO $dbtable (". join(",", @fields) .
	    ") VALUES (". join(",", @values) . ");";
	print "SQL: $sql\n" if $debug > 3;
	if (! $debug ) {
	    $dbh->do($sql) || print "S: '$sql'\n";
	    print ".";
	} else {
	    print ",((lat,$lon),$args{altitude})\n";
	}
    } else {
	print "_";
    }
}

sub tcp_connect {
    my ($iaddr, $paddr, $proto);

    if ($gpsdport =~ /\D/) { $gpsdport = getservbyname($gpsdport, 'tcp') }
    die "No port" unless $gpsdport;
    $iaddr   = inet_aton($gpsdhost)             || die "no host: $gpsdhost";
    $paddr   = sockaddr_in($gpsdport, $iaddr);
    $proto   = getprotobyname('tcp');
    socket(SOCK, PF_INET, SOCK_STREAM, $proto)  || die "socket: $!";
    if (connect(SOCK, $paddr)) {
	# Turn off buffering on SOCK
	my $old_fh = select(SOCK);
	$| = 1;
	select($old_fh);
    } else {
	warn "connect: $!";
    }
    return *SOCK;
}

=head1 valid_point ($latitude, $longitude)

Check if the given point is valid.  Latitude is -90 -> 90, and
longitude is -180 -> 180.

=cut

sub valid_point {
    my ($lat, $lon) = @_;
    return "" if ($lat > 90 || $lat < -90);
    return "" if ($lon > 180 || $lon < -180);
    return 1;
}

sub square {
    my $x = shift;
    return $x * $x;
}
sub distance {
    my ($oldlat, $oldlon, $lat, $lon) = @_;
    my $distance = sqrt(square($lat-$oldlat) + square($lon-$oldlon));
    #print STDERR "Distance: $distance\n";
    return $distance;
}

sub generate_gpx {
    print <<EOF;
<?xml version="1.0" encoding="utf-8"?>
<gpx version="1.1" creator="Free Street Map GPS collector"
        xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
	xmlns="http://www.topografix.com/GPX/1/1"
        xsi:schemaLocation="http://www.topografix.com/GPX/1/1
        http://www.topografix.com/GPX/1/1/gpx.xsd">
  <metadata>
    <name>filename</name>
    <author>author</author>
    <copyright>copyright</copyright>
  </metadata>
  <trk>
    <trkseg>
EOF

    # This SQL code is PostgreSQL specific
    my $sql = "SELECT stamp, g, altitude, satellites, hdilution ".
      "FROM $dbtable ".
#      "WHERE satellites > 0 ". # Is this a sensible limit?
      "ORDER BY stamp, oid";
    my $sth = $dbh->prepare($sql) || die "Bad SQL '$sql'";
    $sth->execute();
    my $lastsecs = 0;
    my ($stamp, $g, $ele, $satellites, $hdop);
    my ($lat, $lon);
    my ($oldlat, $oldlon) = (0,0);
    while ( ($stamp, $g, $ele, $satellites, $hdop) = $sth->fetchrow_array() ) {
	$oldlat = $lat; 
        $oldlon = $lon;
	# Extract lat/lon from $g
	($lat, $lon) = $g =~ m/\((.+),(.+)\)/;

	unless (valid_point($lat, $lon)) {
	    print STDERR "Ignoring invalid point ($lat,$lon) hdop=$hdop sat=$satellites\n";
	    next;
	}

	my $date = &ParseDateString($stamp);
	my $secs = &UnixDate($date,"%s");
	# If the signal was lost or the jump was too long, make new track
	if (($lastsecs > 0 && $secs - $lastsecs > $timelimit)
	    || distance($oldlat, $oldlon, $lat, $lon) > $jumplimit) {
	    # qgis 0.5.0 was unable to handle several track segments
	    # in one track block.  Split each segments into separate
	    # blocks.
	    # Add waypint at the start of each track, to make it
	    # easier to find the tracks in qgis.
	    print <<EOF;
    </trkseg>
  </trk>

  <wpt lat="$lat" lon="$lon">
    <ele>$ele</ele>
  </wpt>

  <trk>
    <trkseg>
EOF
        }
	$lastsecs = $secs;

	# Convert $stamp to ISO date format.  Example: '2004-11-11T00:20:06Z'
	$date =~ s/^(....)(..)(..)(.+)/$1-$2-$3T$4Z/;

	print <<EOF;
      <trkpt lat="$lat" lon="$lon">
        <ele>$ele</ele>
        <time>$date</time>
EOF
        print <<EOF if (defined $hdop);
        <hdop>$hdop</hdop>
EOF
        print <<EOF if (defined $satellites);
        <sat>$satellites</sat>
EOF
	print <<EOF;
      </trkpt>
EOF
    }

    print <<EOF;
    </trkseg>

  </trk>
</gpx>
EOF
}

sub createTable {
    my $dbh = shift;
    # Postgresql table
    $dbh->do("CREATE TABLE $dbtable (".
	     "g          point, ".
	     "altitude   float4,".
	     "hdilution  float4, ".
	     "satellites integer, ".
	     "quality    integer".
	     "stamp      datetime default now()".
	     ");");
    #$dbh->do("CREATE INDEX temppoints_stamp_idx ON temppoints(stamp);");
}