Categories
GIS

Not really getting the Azimuthal Equidistant projection right

Update: thanks to André, this works! See the updated ogr2ogr line.

All I really wanted to do is make a map like this:

This is an Azimuthal Equidistant projection, with me at the centre (of course) and the rest of the world spread out in a fan by distance and bearing. It’s somewhat surprising to find that South Africa is almost directly east of Toronto, and New Zealand to the southwest.

If I had a directional antenna and a rotor, this map would show me where I would have to point the antenna to contact that part of the world. I can’t rotate my dipole (unless I commit some unauthorized local plate tectonics) so I’m stuck with where my antenna transmits and receives best.

The above map was made with AZ_PROJ, a PostScript program of some complexity for plotting world maps for radio use. The instructions for installing and running AZ_PROJ are complex and slightly dated. I got the above output running it through Ghostscript like this:

gs -sDEVICE=pdfwrite -sOutputFile=va3pid_world.pdf az_ini.ps -- az_proj.ps world.wdb

The format of the az_ini.ps file is complex, and I’m glad I’m an old PS hacker to be able to make head or tail of it.

For all its user-hostility, AZ_PROJ is powerful. Here’s a version of the map I wanted all along:

This shows my furthest QSO in each of the 16 compass directions. (You might note that North is empty: my furthest contact in that direction is some 13km away, whether by lack of folks in that sector or dodginess of my antenna.) Contrast that with my Mercator QSO map, and you’ll see that Azimuthal Equidistant is a much better projection for this application.

To show how radically different the world looks to different people, here’s the world according to my mate Rob in Hamilton, NZ:

I’d been trying to use OGR to transform arbitrary shapefiles into this projection. For maps entirely contained within the same hemisphere (so having extent less than ±90° in any cardinal direction), this works:

ogr2ogr -t_srs "+proj=aeqd  +lat_0=43.7308 +lon_0=-79.2647" out.shp in.shp

The lat_0 and lon_0 parameters are just where you want the centre of the map to be. Things get a bit odd if you try to plot the whole world:

The antipodes get plotted underneath, and everything looks messed up. I may have to take my question to GIS – Stack Exchange to see if I can find an answer. Still, for all its wrongness, you can make something pretty, like my whole world Maidenhead locator grid projected this way turns into a rose:

Update: André advised on Stack Exchange that you need to specify the correct radius:

ogr2ogr -t_srs "+proj=aeqd  +lat_0=43.7308 +lon_0=-79.2647 +R=6371000" va3pid-world.shp TM_WORLD_BORDERS_SIMPL-0.3.shp

This works rather well:

aeqd-qgis

 

 

Categories
GIS

Toronto’s Human Centre, part 2: by neighbourhood

Beware throwaway comments; that way overanalysis lies. This was a challenge.

Taking the 2006 neighbourhoods population into account, the human centre of Toronto is at 43.717955°N, 79.389828°W …

… pretty close to the one I’ve already worked out by ward.

Scraping the neighbourhood populations was hard. For the 140 neighbourhoods, the data is stored in a pdf with the URL like http://www.toronto.ca/demographics/cns_profiles/2006/pdf1/cpa124.pdf (in this case, 124; Kennedy Park, represent!). The population number is stored in a table on page 2 of each pdf. I used pdf2xml to convert the files into something parseable.

Of course, the tables weren’t exactly in the same place in every file, so I took a sample of 10% of the files, and worked out the X & Y coordinates of the population box. pdf2xml spits out elements like

<TOKEN sid="p2_s427" id="p2_w417" font-name="arialmt" serif="yes" fixed-width="yes" bold="no" italic="no" font-size="7.47183" font-color="#000000" rotation="0" angle="0" x="299.739" y="122.117" base="129.634" width="22.7692" height="9.94501">17,050</TOKEN>

Yes, I should have used an XML parser, but a Small Matter of Perl got me 126 out of the 140 matching. The rest I keyed in by hand …

Table after the jump.