Category: GIS

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°.