2

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?

1

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.

  • This is very helpful, I didn't realize the data in the shapefile is in meters not degrees. – Mike J Oct 28 '17 at 16:43
  • I've added an edit with code that fixes my original problem but the coordinates are still not correct. – Mike J Oct 28 '17 at 17:21
  • Mercator is often epsg:3857 which is not what you are looking for. – Ian Turton Oct 28 '17 at 18:34
  • Ah ok thanks got it! For anyone who finds this in the future, the cartopy call is cartopy.crs.epsg(26910) to get the projections. – Mike J Oct 29 '17 at 0:16
0

I use QGIS Explorer to open the file and read the metadata of the shapefile and found this information:

Extents In layer spatial reference system units xMin,yMin 483644.83,5449579.62 : xMax,yMax 498313.02,5460349.23

Layer Spatial Reference System +proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

It is clear that the projection is UTM zone 10. And the extents in UTM grid units are as mentioned above.

With this information, the code is written as:

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

fname = r"Path_to\local_area_boundary.shp"
adm1_shapes = list(shpreader.Reader(fname).geometries())
f, ax = plt.subplots()
#f.set_size_inches([10, 6])

# projection from the metadata of shapefile
# UTM(zone=10)
useproj = ccrs.UTM(10)

ax = plt.axes(projection=useproj)
ax.add_geometries(adm1_shapes, useproj,
              edgecolor='black', facecolor='gray', alpha=0.5)

choi = 1  # use either 1 or 2
if choi==1:
    # extents from the metadata of the shapefile
    # use UTM(10) grid coordinates (meters)
    west,south = 483644.83, 5449579.62
    east,north = 498313.02, 5460349.23
    ax.set_extent([west, east, south, north], useproj)
else:
    # or extents with lat long in degrees
    east = -123.0
    west = -123.3
    north = 49.1
    south = 49.4
    ax.set_extent([west, east, south, north], ccrs.Geodetic())

plt.show()

And the resulting plot is:

output plot

Edit1 Note: the coordinate display under the map is for location on the map center.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.