Categories
GIS

Bixi comes to Toronto!

BIXI Toronto is (nearly) here!

Here are the proposed station locations: Bixi_Toronto_shp.zip (Shapefile) or Bixi_Toronto-kml_csv.zip (KML and CSV).

Categories
GIS

More radio amateur grid squares

Toronto, as understood by the Maidenhead Locator system

After yesterday’s post, I went a bit nuts with working out the whole amateur radio grid locator thing (not that I’m currently likely to use it, though). I’d hoped to provide a shapefile of the entire world, but that would be too big for the format’s 2GB file size limit.

What I can give you, though, is:

  • A Perl program that will generate a shapefile of an entire Maidenhead grid field, down to the subsquare level: make_grid.pl. You’ll need Geo::Shapelib to make this work. 324 (= 182) of these files would cover the whole world, and at 8MB or so a pop, things get unwieldy quickly.
  • A Google Earth KML file covering the whole world in 20° by 10° grid fields: Maidenhead_Locator_World_Grid. (If you’re feeling nerdy, here it is in Shapefile format: Maidenhead_Locator_World_Grid-shp).

If anyone would like their grid square in Google Earth format, let me know, or read on …

Making KML Files

Several people have asked, so here’s how you convert to KML. You’ll need the OGR toolkit installed, which comes in several open-source geo software bundles: FWTools/osgeo4w/QGis. Let’s assume we want to make the grid square ‘EN’.

    1. Run make_grid.pl:
      make_grid.pl en
    2. Convert to KML using ogr2ogr:
      ogr2ogr -f KML EN-maidenhead_grid.kml EN-maidenhead_grid.shp
    3. Alternatively, if you just want to extract a square (say EN82), you can use ogr2ogr’s ‘where’ clause to select just the geometry you want:
      ogr2ogr -f KML -where "Square='82'" EN82-maidenhead_grid.kml EN-maidenhead_grid.shp

 

Categories
GIS

Maidenhead Grid Locator, in Perl

Amateur radio operators sometimes use the Maidenhead Grid Locator Square system to identify location. Here’s my attempt at a translator, written in slightly squirrelly Perl:

#!/usr/bin/perl -w
# maidenhead - 6 character locator
# created by scruss/VA3PID on 02011/04/01
# RCS/CVS: $Id: maidenhead,v 1.3 2011/04/03 11:04:38 scruss Exp $

use strict;
sub latlong2maidenhead;
sub maidenhead2latlong;

if ( $#ARGV >= 1 ) {    # two or more args
                        # convert lat/long to grid
  print latlong2maidenhead( $ARGV[0], $ARGV[1] ), "\n";
}
elsif ( $#ARGV == 0 ) {    # one arg
                           # convert grid to lat/long
  print join( ", ", maidenhead2latlong( $ARGV[0] ) ), "\n";
}
else {                     # no args

  # print usage
  print 'Usage: ', $0, ' dd.dddddd ddd.dddddd', "\n",
    ' or', "\n",
    $0, ' grid_location', "\n";
}
exit;

sub maidenhead2latlong {

  # convert a Maidenhead Grid location (eg FN03ir)
  #  to decimal degrees
  # this code could be cleaner/shorter/clearer
  my @locator =
    split( //, uc(shift) );    # convert arg to upper case array
  my $lat      = 0;
  my $long     = 0;
  my $latdiv   = 0;
  my $longdiv  = 0;
  my @divisors = ( 72000, 36000, 7200, 3600, 300, 150 )
    ;                          # long,lat field size in seconds
  my $max = ( $#locator > $#divisors ) ? $#divisors : $#locator;

  for ( my $i = 0 ; $i <= $max ; $i++ ) {
    if ( int( $i / 2 ) % 2 ) {    # numeric
      if ( $i % 2 ) {             # lat
        $latdiv = $divisors[$i];    # save for later
        $lat += $locator[$i] * $latdiv;
      }
      else {                        # long
        $longdiv = $divisors[$i];
        $long += $locator[$i] * $longdiv;
      }
    }
    else {                          # alpha
      my $val = ord( $locator[$i] ) - ord('A');
      if ( $i % 2 ) {               # lat
        $latdiv = $divisors[$i];    # save for later
        $lat += $val * $latdiv;
      }
      else {                        # long
        $longdiv = $divisors[$i];
        $long += $val * $longdiv;
      }
    }
  }
  $lat  += ( $latdiv / 2 );         # location of centre of square
  $long += ( $longdiv / 2 );
  return ( ( $lat / 3600 ) - 90, ( $long / 3600 ) - 180 );
}

sub latlong2maidenhead {

  # convert a WGS84 coordinate in decimal degrees
  #  to a Maidenhead grid location
  my ( $lat, $long ) = @_;
  my @divisors =
    ( 72000, 36000, 7200, 3600, 300, 150 );    # field size in seconds
  my @locator = ();

  # add false easting and northing, convert to seconds
  $lat  = ( $lat + 90 ) * 3600;
  $long = ( $long + 180 ) * 3600;
  for ( my $i = 0 ; $i < 3 ; $i++ ) {
    foreach ( $long, $lat ) {
      my $div  = shift(@divisors);
      my $part = int( $_ / $div );
      if ( $i == 1 ) {    # do the numeric thing for 2nd pair
        push @locator, $part;
      }
      else {              # character thing for 1st and 3rd pair
        push @locator,
          chr( ( ( $i < 1 ) ? ord('A') : ord('a') ) + $part );
      }
      $_ -= ( $part * $div );    # leaves remainder in $long or $lat
    }
  }
  return join( '', @locator );
}

Grid square to lat/long and grid square shapefiles to follow.

Categories
GIS

image georeferencing with QGIS

Quantum GIS (QGIS) has a very powerful image georeferencing module. What that allows you to do is convert screen pixels to a map of an area. The pixels could be from a scanned map or from a screen image. Scanned maps can sometimes be a little distorted, but QGIS’s georeferencer can handle that, within limits.

For an example, say I live in Redickville, ON (I don’t, but some folks do). I’ve heard from North Dufferin Agricultural & Community Taskforce (NDACT) that a large quarry is planned in my area. How close am I going to be from it?

NDACT has a helpful map (2.4 MB PNG image) had a helpful map (which I’ve kept here so you can work through this: MelancthonAerial), but it has no scale. It’s clearly derived from something like Google Maps, so it’s not exactly a free image I can throw about. One thing about georeferencing is that both the source image and the map from which you take you control points affect the licensing of the final map. You’re ending up with a derived work.

There are lots of sources of coordinates for your control points. You could always use Google Maps, but then you’re well down the derived work sinkhole. GeoGratis has a bunch of good data sources, and they are free to use. I’m going to use Toporama‘s digital image maps, as they’re clear and fairly accurate.

Redickville is on Toporama sheet 041A01. I’ve downloaded it as UTM, as it’s easier to measure distances that way. You want to set up your project so it uses the coordinate reference system (CRS) that the control point map uses. Toporama uses EPSG 26917 in the area (easily checked with gdalinfo; it’ll come up with something like AUTHORITY[“EPSG”,”26917″]]), so you should set the project to use that CRS:

And here’s the raster map loaded into QGIS;

Opening up the georeferencer plugin (which is now in the Raster menu, as of QGIS 1.8) gives you a whole lot of blank:

If we open the raster, you can zoom into the area into which you want to put control points. You want to have the highest zoom that the map’s still clear, as the accuracy of your final map depends on how well you placed control points. I’ve chosen road intersections as my control points. Here are the ones I’m going to use:

Point   E-W Road        N-S Road        Note
1       County Rd 21    2nd Line W      Honeywood
2       County Rd 21    4th Line       
3       20th Side Rd    5th Line       
4       15th Side Rd    County Rd 124   just S of where 124 straightens
5       15th Side Rd    Mulmur Townline
6       4th Line        5th Line        4th Line really runs SE-NW

These intersections are fairly well spaced apart, and are clear on both maps. So I choose a pixel on the map in the georeferencer, and a dialogue comes up:

(If you’ve downloaded the MelancthonAerial archive from here, the georeferencing points should load automatically from the “Melancthon Aerial July 22 09.png.points” file. If you want to go through the exercise of manually adding reference points, delete the points file.)

Here I’ve already clicked on the corresponding point on the Toporama map, and the coordinates have been filled in. If you know the coordinates, you can enter them in the boxes. Once you click OK, you have your first control point:

We can add the rest one by one. QGIS’s “Zoom Previous” is really useful for flipping between scales on both the plugin window and the main map. It’s probably a good idea to “Save GCP Points As …” every now and again so you don’t lose your work. Here are all of the points in the georeferencer plugin:

Now you want to modify the Transformation Settings; it’s the little wrench/spanner in the toolbar:

There are lots of options here:

  • Transformation Type: Linear is useful if you’re just adding georeference data to an already computer generated map. Helmert is a simple shear/scale/rotate transformation. The various Polynomial types will correct more gross distortions, but need many control points and can be processor intensive. Thin Plate Spline will distort your map locally to match GCPs; this can work well if your map’s a bit “vernacular”, but if your control points are wrong, your map will end up hilariously melty.
  • Resampling Method: This controls how the output pixels are calculated. Nearest Neighbour is quick and blocky, but useful if you’re mapping an image that has sharp transitions. Cubic maintains more detail, while applying some smoothing at the cost of some detail loss and a fair bit of processing power. The other options can look nice, but eat CPU. This is worth experimenting with many options, as there’s no one solution for all maps.
  • Compression: This controls the file size of the output GeoTiff file. JPEG and Deflate can result in small files, but there’s a chance that other GIS systems can’t read the data. Note that JPEG is lossy, and will lose some detail.

Don’t forget to set the output raster file, and make sure that the target SRS matches the CRS (EPSG 26917) you chose earlier.

Helmert is a useful transformation type to check your work. The plugin plots the residual difference between the two sets of map control points as red lines. Here, I’ve clearly made an error in point ID 2:

If you unselect the bad point, the plugin quickly calculates the residuals again. When this one bad point is removed, the residuals drop down to less than 1.0 for each point. I’ve got enough points for a Helmert transformation with 5 GCPs, but I’d probably want some more points for more complex transformations.

Once you hit “Start Transformation”, QGIS will create your referenced image. If you’ve chosen a simple linear transformation, it won’t create a GeoTiff file, but just a world file for your image.

So here’s the georeferenced image overlaid on the Toporama map, with a bit of transparency. It’s quite a good match:

(the above’s an image saved straight from QGIS. It helpfully creates a world file too, so here it is: redickville.jpgw)

So to answer the hypothetical question, Redickville’s pretty well surrounded by this proposed quarry.

Categories
GIS

Garmin vs Arduino GPS dance-off

So I built a GPS logger. It’s an Arduino Uno with the SparkFun microSD and GPS shields. No display, massively unwieldy, but it does spit out a track point every second and logs it fairly reliably to microSD. It works as well as my expensive old Garmin GPSMap 60Csx, as a walk around the block shows:

Excuse the noise at the top of the trace; both units were brought inside. While they keep a tenuous lock, they certainly don’t give much accuracy. I think the Arduino did rather well; certainly better than my BlackBerry

Now all I need is a case and a more useful power supply. I was disappointed that a USB MintyBoost charger didn’t seem to work for long, causing the GPS to lock up.

Categories
GIS

fixing garmin file dates

While the Garmin GPSMap 60csx is a lovely unit, it saves its tracks on the card with a date just slightly younger than I am. The following Unix one-liner will correct the file dates to the actual dates the data were collected:

for f in 20??????.gpx; do touch -t ${f%.gpx}2359.59 $f; done

I remember having a really awesome reason for making the time for each file 23:59:59, but I’ve completely forgotten what it was. Since all I remember was the awesomeness, I see no reason to change it …

Update 2013-09-25: changed the call to basename with a Bash internal string function.

Categories
GIS

The SpatiaLite Cookbook

Alessandro‘s new SpatiaLite Cookbook is rather good.

Categories
GIS

Interfacing the Pharos GPS-500 to Mac OS X

Update 2015: This is very old. I don’t know the status of a decent OS X 10.10 driver for the Prolific PL-2303 in this GPS.

I’m stopping short of calling this “Using the Pharos GPS-500 to Mac OS X” because all I’ve been able to do is read raw NMEA sentences from the device. But that might be of use to you.

The Source is clearing out copies of Microsoft Street & Trips with GPS 2008 for $20. The GPS is a very simple USB Pharos GPS-500, which uses the SiRF III chipset. Between the USB cable and the GPS is a small black box which looks suspiciously like a serial to USB converter. I have no use for the software, but the GPS is a bargain, considering a similar bare unit costs $60.

Plugging it into a Mac does nothing beyond being recognized as a “USB-Serial Controller D” from Prolific Technology, Inc. The ancient driver on Pharos’ website identifies the serial chipset as the Prolific PL2303. The only driver I could get to work with Snow Leopard was the failberg/osx-pl2303, a fork of an earlier project from Sourceforge. You’ll know if it’s working when you get a device called /dev/tty.PL2303-something appear.

Reading the data’s pretty simple if you have GNU Screen installed. I entered the following command:

screen /dev/tty.PL2303-12345678^XX^D?^XX 4800

and very quickly started to get NMEA data scrolling in the terminal:

$GPRMC,003322.000,A,4343.8349,N,07915.8845,W,0.38,112.13,211210,,,A*7B
$GPGGA,003323.000,4343.8351,N,07915.8838,W,1,05,2.8,174.3,M,-35.1,M,,0000*65
$GPGSA,A,3,09,18,14,21,22,,,,,,,,6.9,2.8,6.3*34
$GPRMC,003323.000,A,4343.8351,N,07915.8838,W,0.89,122.04,211210,,,A*76
$GPGGA,003324.000,4343.8351,N,07915.8843,W,1,06,1.5,176.0,M,-35.1,M,,0000*62
$GPGSA,A,3,27,09,18,14,21,22,,,,,,,3.4,1.5,3.1*30
$GPGSV,3,1,12,18,78,082,28,22,61,307,27,09,56,076,30,27,40,053,19*7F
$GPGSV,3,2,12,14,38,254,18,21,32,180,21,15,27,062,16,26,12,059,*7E
$GPGSV,3,3,12,19,10,302,,12,07,119,,06,04,261,,03,02,272,*7E

To stop screen, type Control-A K. Do not just unplug the GPS, as you risk your machine crashing.

These NMEA sentences can easily be converted to GPX.

Categories
GIS

Magellan Triton 400 micro-review

My office’s Garmin was stuck in a branch office last week, and we needed a GPS for the next morning, so we got a cheapo Magellan Triton 400 for $90 at Future Shop. I think that’s a clearance price, and none of the big-box retailers still carry it.

The Triton 400’s a Windows CE unit with a surprisingly good display for the price. I only had a little time to set it up and test it briefly outside the office, so all I can do is give you is first impressions.

Pros:

  • Cheap!
  • SIRFstarIII chipset for reasonably fast/accurate acquisition
  • SDHC card (worked with my 4GB card)
  • Bright display
  • Works with GPSBabel (as Magellan’s VantagePoint obviously installs a copy)
  • Doesn’t route

Cons:

  • Wouldn’t acquire any position until I updated the firmware (at which time I discovered that it basically updates a full Windows CE image from an archive)
  • Weird proprietary USB cable
  • Tiny buttons that aren’t very positive
  • Eats batteries
  • Overly simplistic menu structure makes it hard to set up
  • VantagePoint is buggy, and will repeatably crash under certainly (admittedly rare) menu items
  • Only works under Windows; the USB protocol is proprietary
  • Hard limit of 5000 points per track, and track logging stops when this limit is reached
  • SD card is only usable for maps and geotagged photos, not track storage.

I should be able to play with it in more detail in the new year. It was cheap, but not all cheap things are good.

More on the Magellan – OpenStreetMap Wiki page.

Categories
GIS

Got a receiver inside my head – writing Spectrum Direct data as CSV

Industry Canada publishes the locations of all licensed radio spectrum users on Spectrum Direct. You can find all the transmitters/receivers near you by using its Geographical Area Search. And there are a lot near me:

While Spectrum Direct’s a great service, it has three major usability strikes against it:

  1. You can’t search by address or postal code; you need to know your latitude and longitude. Not just that, it expects your coordinates as a integer of the format DDMMSS.
  2. It’s very easy to overwhelm the system. Where I live, I can pretty much search for only 5km around me before the system times out.
  3. The output formats aren’t very useful. You can either get massively verbose XML, or very long line undelimited text, and neither of these are very easy to work with.

Never fear, Perl is here! I wrote a tiny script that glues together Dave O’Neill‘s Parse::SpectrumDirect::RadioFrequency module (which I wonder if you can guess what it does?) to Robbie Bow‘s  Text::CSV::Slurp module. The latter is used to blort out the former’s results to a CSV file that you can load into any GIS/mapping system.

Here’s the code:

#!/usr/bin/perl -w
# spectest.pl - generate CSV from Industry Canada Spectrum Direct data
# created by scruss on 02010/10/29 - for https://glaikit.org/

# usage: spectest.pl geographical_area.txt &gt; outfile.csv

use strict;
use Parse::SpectrumDirect::RadioFrequency;
use Text::CSV::Slurp;
use constant MINLAT =&gt; 40.0;    # all of Canada is &gt;40 deg N, for checking

my $prefetched_output = '';

# get the whole file as a string
while (&lt;&gt;) {
 $prefetched_output .= $_;
}

my $parser = Parse::SpectrumDirect::RadioFrequency-&gt;new();

# magically parse Spectrum Direct file
$parser-&gt;parse($prefetched_output) or die &quot;$!\n&quot;;
my $legend_hash = $parser-&gt;get_legend();    # get column descriptions
my @keys        = ();
foreach (@$legend_hash) {

 # retrieve column keys in order so the output will resemble input
 push @keys, $_-&gt;{key};
}

# get the data in a ref to an array of hashes
my $stations = $parser-&gt;get_stations();

my @good_stations = ();

# clean out bad values
foreach (@$stations) {
 next if ( $_-&gt;{Latitude} &lt; MINLAT );
 push @good_stations, $_;
}

# create csv file in memory then print it
my $csv = Text::CSV::Slurp-&gt;create(
 input       =&gt; \@good_stations,
 field_order =&gt; \@keys
);
print $csv;
exit;

The results aren’t perfect; QGis boaked on a file it made where one of the records appeared to have line breaks in it. It could filter out multiple pieces of equipment at the same call sign location. But it works, mostly, which is good enough for me.