Toronto’s Angles

As many Canadians will tell you, most Torontonians are a bit skewed. While others might have different explanations, I put it down to our road grid. What we call north in the city is actually some 16.7° west of north.

While we do have a perpendicular grid, it’s twisted to the left (steady, Edmontonians). I attempted to work out what the average angle is.

I eyeballed several long straight streets, and picked out representative intersections. Finidng the coordinates of these intersections took a surprisingly long time, as each street has several road sections. To find an intersection, each section on one road has to be queried against each section of the other. It gets there, eventually.

To save the tedium of showing the queries, here are the results. The angles are anticlockwise from the X-axis.

N-S	 Intersection1 Intersection2 Angle
======== ============= ============= ======
Bathurst King	       Dupont	     16.184
Yonge	 King	       Rosehill	     16.783
Woodbine Queen	       Plains	     17.049

E-W	 Intersection1 Intersection2 Angle
======== ============= ============= ======
Danforth Greenwood     Main	     16.736
Eglinton Dufferin      Bayview	     16.159
Steeles	 Dufferin      Kennedy	     17.194

I was surprised that SpatiaLite didn’t have any angle measurement functions built in, so I had to feed the output through geod. So for instance, for Steeles Avenue, the intersection of Steeles West and Dufferin is at 43.787229°N, 79.470285°W, and the intersection of Steeles East and Kennedy is 43.823636°N, 79.307265°W. Feeding these numbers to geod:

echo '43.787229 -79.470285 43.823636 -79.307265'| geod +ellps=WGS84 -f "%.3f" -p -I +units=m

This spits out three numbers: 72.806    252.918    13727.395. The first is the bearing from the first point to the second; this is the number we’re interested in, and it’s (90-82.806) = 17.194° clockwise from E. The second is the bearing from the second point back to the first. The last is the distance in metres; Steeles is a long, straight street.

Averaging the numbers for the streets I chose gets 16.684°; just the thing for Torontohenge. For a visual axis, maybe using Steeles’s 17.194° could be better, as it’s a line above the city. You can set the map rotation angle in QGIS’s map composer, and get a better aligned city:

but where am i, really?

In my first post I asked where am i? The gps in my phone said I was standing at 43.73066°N, 79.26482°W. In real life, I was standing at the junction of Kenmark Blvd and Chevron Cres.

With the Open Toronto centreline data, I can check the location of road intersections:

select distinct (
astext (
transform (
intersection ( r1.geometry, r2.geometry ), 4326
 ) ) )
from centreline as r1, centreline as r2
where r1.lf_name = 'KENMARK BLVD' and r2.lf_name
 = 'CHEVRON CRES'

which returns:

NULL
POINT(-79.262423 43.730019)
POINT(-79.264815 43.730591)

Hmm; three answers. NULL I can’t answer; I’ll put it down as an exhortation to Be Here Now. The second is given a clue by one of the street names: Chevron Crescent – first thing I learned when I had my newspaper round is that a crescent’s going to end you up back on the same road you started from. The last one, though, agrees to 5(ish) decimal places – well within the accuracy of my simple GPS.

Update: This query gets really, really slow on long streets. Both Chevron and Kenmark only have two road segments. Kennedy and Steeles East each have over 90.

Finding the exact geographic centre of Toronto

Spacing has alluded to it. The Ontario Science Centre makes bold (and incorrect) claims about it. But here’s the real deal.

  • If you consider Toronto to be defined by its city wards, the centre of Toronto lies at 43.725518°N, 79.390531°W.
  • If you consider Toronto to be defined by its neighbourhoods, the centre of Toronto lies at 43.726495°N, 79.390641°W.

You can work this out in one line of SQL. By combining all the wards or neighbourhoods into one union shape (SpatiaLite uses the GUnion() function), and then calculating the centroid, that’s the centre of the city:

select astext(transform(centroid(gunion(geometry)),4326)) from wards

To get the results in a more human-friendly format, I transformed it to WGS84 (EPSG SRID 4326), and used astext() to get it in something other than binary.

Update, 18 March: great minds, etc … Torontoist just posted the similar In Search of Toronto’s Geographic Centre.

using CSV as a virtual data source

While we already know how to make trivial shapefiles with shapelib, sometimes that’s too tedious. Very frequently I get data in Comma Separated Value (CSV) format, and reliably importing/converting it can be a pain.

Here’s our sample CSV file, library_test.csv:

"Easting","Northing","Library"
625539.6,4837170.9,"Dufferin St. Clair"
625862.0,4838241.1,"Oakwood Village"
626251.0,4835287.2,"Bloor Gladstone"
626671.7,4836922.6,"Davenport"
627227.2,4840006.4,"Forest Hill"

ogr has a CSV driver. In its documentation the Virtual Format driver is touched upon. This allows you to set up a data definition file, especially useful if you read the same format frequently.

Here’s the VRT file for that CSV:

<OGRVRTDataSource>
    <!-- note that OGRVRTLayer name must be basename of source file -->
    <OGRVRTLayer name="library_test">
        <SrcDataSource>library_test.csv</SrcDataSource>
        <GeometryType>wkbPoint</GeometryType>
        <!-- your SRS goes here; I used EPSG SRID -->
        <LayerSRS>EPSG:2958</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="Easting" y="Northing"/>
   </OGRVRTLayer>
</OGRVRTDataSource>

Your CSV file will now behave like a shapefile, or indeed any other geo-format that OGR understands. QGIS is a bit picky – it doesn’t seem to always work out the path of the source file.

To prove these are real coordinates, here’s what I did to make a Google Earth KML file:

ogr2ogr -f KML -t_srs EPSG:4326 library_test.kml library_test.vrt -dsco NameField=Library

Technically, you don’t need to specify the SRS for KML output as it only supports EPSG:4326, but I found you got trivially different results if it was omitted.

Try this in Google Earth: library_test.kml

ward maps: kinda working, sorta

Now I’ve sorted out formatting the labels and scraping the data, I should be almost ready to produce a pretty map.

Well, almost. The DBF component of a shapefile seems somewhat resistant to adding a column, and SQLite doesn’t seem very happy with its ALTER TABLE ADD COLUMN ... syntax.

As usual, I needed to create the database table from the shapefile. I’m not bothered about CRS, so I used -1.


.read init_spatialite-2.3.sql ASCII

.loadshp TCL3_ICITW Wards CP1252 -1

alter table wards add column candidates integer

I had mixed success getting data to load into this new column. So I improvised.

!!! WARNING: EGREGIOUS MISUSE OF DATA FOLLOWS !!!

(Sensitive readers are advised to look away)

There’s a seeming unused numeric column SHAPE_LEN in the table. As my new candidates column was coming up with occasional nulls, I cheated:


UPDATE Wards set shape_len=3 where scode_name="1"

UPDATE Wards set shape_len=1 where scode_name="2"

UPDATE Wards set shape_len=0 where scode_name="3"

...

UPDATE Wards set shape_len=3 where scode_name="44";

I then added SHAPE_LEN as the label, and defined a range based colour gradient for the wards in QGIS’s layer properties:

And this is how it looks:

Another partial success, as Professor Piehead would say.

closer to ward maps: scraping the data

Toronto publishes its candidates here  http://app.toronto.ca/vote2010/findByOffice.do?officeType=2&officeName=Councillor in a kind of tabular format. All I want to do is count the number of candidates per ward, remembering that some wards have no candidates yet.

Being lazy, I’d far rather have another program parse the HTML, so I work from the formatted output of W3M. It’s relatively easy to munge the output using Perl. From there, I hope to stick the additional data either into a new column in the shapefile, or use SpatiaLite. I’m undecided.

My dubious Perl script:


#!/usr/bin/perl -w
# ward_candidates - mimic mez ward map
# created by scruss on 02010/03/01
# RCS/CVS: $Id$

use strict;
my $URL =
'http://app.toronto.ca/vote2010/findByOffice.do?officeType=2&officeName=Councillor';
my $stop = 1;

my %wards;
for ( 1 .. 44 ) {
 $wards{$_} = 0;    # initialise count to zero for each ward
}

open( IN, "w3m -dump \"$URL\" |" );
while (<IN>) {
 chomp;
 s/^\s+//;
 next if (/^$/);
 $stop = 1 if (/^Withdrawn Candidate/);
 unless ( 1 == $stop ) {
 my ($ward) = /(\d+)$/;
 $wards{$ward}++;    # increment candidate for this ward
 }
 $stop = 0 if (/^City Councillor/);
}
close(IN);

foreach ( sort { $a <=> $b } ( keys(%wards) ) ) {
 printf( "%2d\t%2d\n", $_, $wards{$_} );
}

exit;

which outputs the following (header added for clarity):

Ward Candidates
==== ==========
 1     3
 2     1
 3     0
 4     0
 5     1
 6     1
 7     7
 8     3
 9     2
10     3
11     2
12     3
13     1
14     4
15     3
16     1
17     2
18     4
19     6
20     2
21     1
22     1
23     1
24     0
25     2
26     3
27    12
28     3
29     6
30     3
31     3
32     2
33     1
34     0
35     5
36     2
37     2
38     2
39     1
40     2
41     1
42     5
43     3
44     3

Labelling: harder than it looks

I’m rather taken with Mez’s rather neat Toronto ward candidate maps. I wonder if I could reproduce them (semi-)automatically?

As a start, here’s the Toronto Wards layer, rendered in QGIS with the ward number as a label:

You’ll notice that something is quite off. It looks like QGIS uses the centre of the minimum bounding rectangle of a polygon as the label point. While this is okay for nice regular shapes, weird glaikit shapes end up with the label outside the boundary. Not good.

I was about to give up on this completely, when I saw QGIS’s “Labeling” [sic] plugin. What it does is work out a variety of better visual positions for your labels. Here’s the setting I chose:

The result is much more pleasing:

Much better.

making trivial shapefiles with shapelib

While most GIS applications work with delimited text inputs, sometimes you just have to have a shapefile. Amongst many other things, Frank Warmerdam wrote the Shapefile C Library, which comes with a few simple tools. I suspect Frank meant the little utilities to be code samples that wouldn’t see much use, but they do the job.

Let’s take the coordinates 43.73066°N, 79.26482°W from my first entry. I will make a single point shapefile with this coordinate.

First you have to make the SHP file and the DBF database:

dbfcreate junction -s Name 16
shpcreate junction point

This makes an empty shapefile for storing points, with one string field ‘Name’ of width 16 characters.

Now you have to add your point – this takes two stages, adding the database row, and then adding the geometry:

dbfadd junction.dbf 'Chevron/Kenmark'
shpadd junction.shp -79.26482 43.73066

And that’s it – you’ve made a trivial shapefile.