How do I match the projection of my cartopy map with that of a shapefile?

Question

I am trying to synthesise the projections of a coastlines() map with that of a shapefile, whose .prj file says:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",
SPHEROID["WGS_1984",6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

My attempt is:

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

# set up a map with coastlines around Auckland:
plt.figure(figsize=(10, 10))
platecarree = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84'))

ax = plt.axes(projection=platecarree)
extent = [174.25, 175.25, -37.5, -36.5]
ax.set_extent(extent)
ax.coastlines('10m',color='red')

# read in shapefile and plot the polygons:
shp2 = shapereader.Reader('auckland_geology_wgs84gcs.shp')
formations = shp2.records()

for formation in formations:
     # plot water blue, and all other rocks yellow
     if formation.attributes['MAIN_ROCK'] == b'                                ':
         ax.add_geometries(formation.geometry, ccrs.PlateCarree(),facecolor='blue',alpha=.1)
     else:
         ax.add_geometries(formation.geometry, ccrs.PlateCarree(), facecolor='yellow',alpha=.1)
plt.show()

I tried giving the globe parameter in my platecarree definition the radius and inverse flattening from the prj file, but I didn't see any change to the output if I set or even varied those numbers.

In addition, with the defined "platecarree" projection (with the call to the globe with WGS84) as the crs in the add_geometries calls, my output is blank.

As is, the result looks to me like a projection mismatch

Thomas Lecocq Thomas Lecocq · Accepted Answer · 2016-07-04T10:01:45

I've tried to reproduce your problem using QGIS and data downloaded from Natural Earth (10m coastlines) and from GADM (NZ adm0 level). It looks like the NE10m coastlines are the culprit ! The GADM aligns perfectly with your geology layer, while the NE10m is off (and deformed). screenshot of QGIS with Geological map & coastlines

How do I match the projection of my cartopy map with that of a shapefile?

1 Answers