Use ccrs.epsg() to plot zipcode perimeter shapefile with EPSG 4326 coordinate system

Question

I've obtained a shapefile of zipcode perimeters from here and would like to plot them on top of a Cartopy map, as I did in this example.

According to the source, this data is in EPSG 4326 coordinate system. When I attempt to plot the data

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt 
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

# Create a Stamen terrain background instance
stamen_terrain = cimgt.Stamen('terrain-background')
fig = plt.figure(figsize = (mapsize,mapsize))
ax = fig.add_subplot(1, 1, 1, projection=stamen_terrain.crs)

# Set range of map, stipulate zoom level
ax.set_extent([-122.7, -121.5, 37.15, 38.15], crs=ccrs.Geodetic())
ax.add_image(stamen_terrain, 8, zorder = 0)

# Add city borders
shape_feature = ShapelyFeature(Reader(shapefile).geometries(), ccrs.epsg(4326), 
                linewidth = 2, facecolor = (1, 1, 1, 0), 
                edgecolor = (0.3, 0.3, 0.3, 1))
ax.add_feature(shape_feature, zorder = 1)

I see the following error:

ValueError: EPSG code does not define a projection

Probably related -- when I look at the ccrs.epsg() function, it says that this EPSG code is not supported

help(ccrs.epsg)
Help on function epsg in module cartopy.crs:

epsg(code)
    Return the projection which corresponds to the given EPSG code.

    The EPSG code must correspond to a "projected coordinate system",
    so EPSG codes such as 4326 (WGS-84) which define a "geodetic coordinate
    system" will not work.

    Note
    ----
        The conversion is performed by querying https://epsg.io/ so a
        live internet connection is required.

Given this result, I also tried using ccrs.Geodetic():

# Add city borders
shape_feature = ShapelyFeature(Reader(shapefile).geometries(), ccrs.Geodetic(), 
                linewidth = 2, facecolor = (1, 1, 1, 0), 
                edgecolor = (0.3, 0.3, 0.3, 1))
ax.add_feature(shape_feature, zorder = 1)

But this also fails to show the zipcode perimeters, and shows this warning message:

UserWarning: Approximating coordinate system <cartopy._crs.Geodetic object at 0x1a2d2375c8> with the PlateCarree projection.
  'PlateCarree projection.'.format(crs))

I've tried ccrs.PlateCarree() too, but no luck. Please help!

From the content of `ZIPCODE.prj`, the projection is `NAD_1983_UTM_Zone_10N`. So it is not `EPSG 4326`. For detail, see https://epsg.io/26910 . — swatchai, Jul 15 '19 at 17:01
Thanks. I've also tried `ccrs.epsg(26910)` -- it does not return an error, but also does not plot the boundaries -- any thoughts on what to try next? — Michael Boles, Jul 15 '19 at 17:56

score 2 · Accepted Answer · answered Jul 18 '19 at 02:15

To plot data from different sources together, one must declare correct coordinate reference system for each data set. In the case of shapefile, you can find it in its accompanying xxx.prj file.

Here is the working code:

import cartopy.io.img_tiles as cimgt 
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

shapefile_name= "./data/ZIPCODE.shp"
mapwidth, mapheight = 8, 8
pad = 0.25

stamen_terrain = cimgt.Stamen('terrain-background')
stm_crs = stamen_terrain.crs

fig = plt.figure(figsize = (mapwidth, mapheight))
ax = fig.add_subplot(1, 1, 1, projection=stm_crs)  #world mercator

# Set extent of map
ax.set_extent([-123.3-pad, -121.5+pad, 37.05-pad, 38.75+pad], crs=ccrs.Geodetic())
# Plot base map
ax.add_image(stamen_terrain, 8, zorder=0)

# Add polygons from shapefile
# Note: the use of `ccrs.epsg(26910)`
shape_feature = ShapelyFeature(Reader(shapefile_name).geometries(), ccrs.epsg(26910))

# You can choose one of the 2 possible methods to plot
# ... the geometries from shapefile
# Styling is done here.
method = 1
if method==1:
    # iteration is hidden
    ax.add_feature(shape_feature, facecolor='b', edgecolor='red', alpha=0.4, zorder = 15)
if method==2:
    # iterate and use `.add_geometries()`
    # more flexible to manipulate particular items
    for geom in shape_feature.geometries():
        ax.add_geometries([geom], crs=shape_feature.crs, facecolor='b', edgecolor='red', alpha=0.4)

plt.show()

The output plot:

Your `method==2` example is quite helpful. With this, I was able to color polygons by value given in a separate dataframe: `for counter, geom in enumerate(shape_feature.geometries()):` `if data['Area'][counter] < 50:` `if data['Population'][counter] > 500:` `ax.add_geometries([geom], crs=shape_feature.crs,` `facecolor=color[counter], edgecolor='k', alpha=0.8)` `else:` `continue` ` — Michael Boles, Jul 27 '19 at 20:15

Use ccrs.epsg() to plot zipcode perimeter shapefile with EPSG 4326 coordinate system

1 Answers1