[GIS] Shapefiles won’t display with cartopy

cartopypythonshapefile

I'm trying my first GIS project with cartopy and matplotlib and I'm having some trouble displaying some shapefiles.

The shapefiles were downloaded from here: http://data.vancouver.ca/datacatalogue/localAreaBoundary.htm

The code I'm trying to use has worked for other shapefiles but not for this set.

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt

fname = "local_area_boundary.shp"
adm1_shapes = list(shpreader.Reader(fname).geometries())

ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_geometries(adm1_shapes, ccrs.PlateCarree(),
                  edgecolor='black', facecolor='gray', alpha=0.5)
east = -123.0
west = -123.3
north = 49.1
south = 49.4
ax.set_extent([west, east, south, north], ccrs.PlateCarree())
plt.show()

Again, I'm 100% new to GIS and cartopy (but fairly experienced in Python) and I'm sure I'm missing something easy. Any help debugging this would be great.

Edit: I've modified the code so the shapefile is displayed, but it now has the wrong lat/lon coordinates

fname = 'local_area_boundary.shp'
adm1_shapes = list(shpreader.Reader(fname).geometries())

f, ax = plt.subplots()
ax = plt.axes(projection=ccrs.Mercator())

ax.add_geometries(adm1_shapes, ccrs.Mercator(),
                  edgecolor='black', facecolor='white', alpha=0.5)
a = 483644.82757398
b = 5449579.62327763
c = 498313.01571369
d = 5460349.22564139
ax.set_extent([a,c,b,d ], ccrs.Mercator())

Shapes are plotted but lat/lon coordinates are wrong
I would expect the lat/lon to be centred around 49.2N/-123.1E, but instead it's 44.1N/4.3E. Can you tell if this is just bad data, or if I'm still doing the transformation wrong?

Best Answer

That download page is lying when it says Coordinate system N/A, the data is in fact projected as it has a .prj file which contains:

PROJCS["NAD83 / UTM zone 10N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-123],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","26910"]]

So your data is in metres while your map is in degrees. So you will need to change the following lines:

ax = plt.axes(projection=ccrs.PlateCarree())
east = -123.0
west = -123.3
north = 49.1
south = 49.4
ax.set_extent([west, east, south, north], ccrs.PlateCarree())

I assume cartopy has some way to set an ESPG code, in which case you want EPSG:26910.

Related Question