Tag: crs

  • tale of two cities: coordinate reference systems, and what on earth is the maywood tot lot?

    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.

    Once we do that, things align as they should. Here’s my neighbourhood, with its parks

    But things are still off if you’re querying the SQL directly:

    
    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:

    src/2019/address-points
    src/2019/business-improvement-areas
    src/2019/city-wards
    src/2019/parks
    src/2019/solid-waste-management-districts
    src/2019/toronto-centreline
    src/26717/food-banks
    src/26717/neighbourhoods
    src/26717/places-of-worship
    src/26717/priority-invest-neighbourhoods
    src/26717/rent-bank-zones
    src/26717/rent-banks
    src/26717/transit-city
    dest/2958/address-points
    dest/2958/business-improvement-areas
    dest/2958/city-wards
    dest/2958/food-banks
    dest/2958/neighbourhoods
    dest/2958/parks
    dest/2958/places-of-worship
    dest/2958/priority-invest-neighbourhoods
    dest/2958/rent-bank-zones
    dest/2958/rent-banks
    dest/2958/solid-waste-management-districts
    dest/2958/toronto-centreline
    dest/2958/transit-city

    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
    
    NAME Area_ha
    MAYWOOD TOT LOT 0.7
    KITCHENER PARK 0.8
    GLEN SHEPPARD PARK 1.0
    GREYSTONE PARK 1.1
    CORVETTE PARK 2.5
    MID-SCARBOROUGH C.C. & ARENA 2.9
    FARLINGER RAVINE 2.9

    (note how I sneakily used the round() function to avoid too many decimal places?)