For reasons that are not particularly clear, the Toronto.ca|Open data is in two different coordinate reference systems (CRS), MTM 3 Degree Zone 10, NAD27 (EPSG 2019) and UTM 6 Degree Zone 17N NAD27 (EPSG 26717). This confuses QGIS even if you’ve input the proper SRIDs into SpatiaLite. The image above shows two apparent Torontos, one in each of the CRSs.
What you have to do is go to to the Project Properties, select the Coordinate Reference System (CRS) tab, and “Enable ‘on the fly’ CRS transformation”. This will line those city layers right back up.
select Distance(Parks.geometry, Neighbourhoods.geometry)/1000 as Distance_km from Parks, Neighbourhoods where Parks.name='CORVETTE PARK' and Neighbourhoods.hood='Kennedy Park'
which returns a distance of over 314 km. That’s not right.
So we need to transform the geometries to the same CRS.
!!! NB: I might be doing the next bit wrong. CRS transformation is subtle. I’m not, particularly.
The OGR Simple Feature Library is your friend. It can convert pretty much any geo format to another, and can transform coordinates between systems. In exchange for this power, it wants your soul is rather complex.
I’ve chosen to use NAD83(CSRS) / UTM zone 17N (EPSG 2958) for my Toronto maps. It’s fairly accurate over the whole city area. To convert the Parks and Neighbourhoods shape files:
ogr2ogr -s_srs EPSG:2019 -t_srs EPSG:2958 dest/2958/parks/TCL3_UPARK.shp src/2019/parks/TCL3_UPARK.shp ogr2ogr -s_srs EPSG:26717 -t_srs EPSG:2958 dest/2958/neighbourhoods/Neighbourhoods.shp src/26717/neighbourhoods/Neighbourhoods.shp
Note that it wants the destination file first, then the source. Haven’t seen that order since PIP under CPM/2.2. I was also a bit nerdly, and arranged the files in directories by SRID:
If we load the transformed shapefiles into Spatialite, and run that query again, it comes out with the correct distance: 0.0 km, as Corvette Park is in the Kennedy Park Neighbourhood.
Now we can run a proper query: what parks are in Kennedy Park, and what are their areas?
select tp.name, round(Area(tp.geometry)/10000,1) as Area_ha from Parks as tp, Neighbourhoods as tn where tn.hood='Kennedy Park' and within(tp.geometry, tn.geometry) order by Area_ha
|MAYWOOD TOT LOT||0.7|
|GLEN SHEPPARD PARK||1.0|
|MID-SCARBOROUGH C.C. & ARENA||2.9|
(note how I sneakily used the round() function to avoid too many decimal places?)