Categories
GIS

three norths to choose from

North is a slippery thing, a trickster. There are at least three norths to choose from:

  • Magnetic north – where the compass needle points. Currently way in the northwest of Canada, it’s scooting rapidly towards Siberia. Come baaack, magnetic north! Was it something we said?
  • True north – straight up and down on the axis of this planet.
  • Grid north – the vertical lines on your map. Due to the earth’s curvature, there’s only one line on a grid that exactly coincides with true north. Away from that line, there’s a correction that you need to apply.

That grid to true correction is what I’m looking at today. In my work, I need to work out where shadows of wind turbines fall. I work in grid coordinates, so I need to correct for true north sun angles. Here’s a shell script that uses familiar programs invproj and geod to work out the local true north correction:

#!/bin/sh
# gridnorth.sh - show the grid shift angle at a certain point
# created by scruss on Sun 4 Apr 2010 18:19:29 EDT
# $Id: gridnorth.sh,v 1.2 2010/04/05 12:23:57 scruss Exp $

# you might want to change this
srid=2958

# two arguments mandatory, third optional
if
 [ $# -lt 2 ]
then
 echo Usage: $0 x y [srid]
fi

# truncate arguments to integers so the shell maths will work
x=`echo $1 | sed 's/\.[0-9]*//;'`
y=`echo $2 | sed 's/\.[0-9]*//;'`
# calculate points 1000 units (grid) south and north of the chosen point
lowy=$((y-1000))
highy=$((y+1000))

if
 [ $# -eq 3 ]
then
 srid=$3
fi

min=`echo $lowy  $x foo | invproj -r -f "%.6f" +init=EPSG:$srid`
max=`echo $highy $x foo | invproj -r -f "%.6f" +init=EPSG:$srid`
minlon=`echo $min | awk '{print $1;}'`
minlat=`echo $min | awk '{print $2;}'`
maxlon=`echo $max | awk '{print $1;}'`
maxlat=`echo $max | awk '{print $2;}'`
echo $minlat $minlon $maxlat $maxlon |\
 geod +ellps=WGS84 -f "%.3f" -p -I +units=m |\
 awk '{print $1;}'

The script is called with two arguments (the easting and the northing) and an optional third argument, the SRID for the current projection:

$ ./gridnorth.sh 639491 4843599 26717
1.198

The value returned is the angle (clockwise) that grid north differs from true north. Unless you’re using the same SRID that I’ve put in the script as default, you probably want to specify (or change) it.

So how do I know if this is even remotely correct? If you look at any CanMatrix scanned topographic map (either Print Ready or Georeferenced). there’s a declination indicator shown:

That’s the map eastern Toronto is in; its sheet reference is 030M11. We’ll ignore the magnetic declination for now, as 1984 is way out of date, and there is much to be said on geomagnetism that I don’t understand.  The map says its centre is 1° 12′ clockwise of true north. The centre of the map is approximately 641500E, 4831000N (UTM Zone 17N NAD27; EPSG 26717; or in real terms, in the lake, a few klicks SE of Ashbridges Bay) so:

$ ./gridnorth.sh 641500 4831000 26717
1.209

1° 12′ is 1.2, so I think I’m officially pretty durn close.

To show that grid skew changes across the map, the upper right hand corner of the map is 660000E, 4845000N, which is off by 1.375°.

Categories
GIS

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: