Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

County boarders in Cartopy

Tags:

python

cartopy

How do you plot US county borders in Cartopy?

It's very straight forward to plot state and country boundaries

ax.add_feature(cfeature.BORDERS.with_scale('50m'))
ax.add_feature(cfeature.STATES.with_scale('50m'))

But I can't seem to find a similar method to add county boundaries. This was one of the nice things about Basemap.

like image 231
blaylockbk Avatar asked Jun 29 '18 17:06

blaylockbk


People also ask

What is Cartopy in Python?

Cartopy is a Python package designed for geospatial data processing in order to produce maps and other geospatial data analyses. Cartopy makes use of the powerful PROJ, NumPy and Shapely libraries and includes a programmatic interface built on top of Matplotlib for the creation of publication quality maps.


1 Answers

Given cartopy's ability to draw shapefiles, this question essentially boils down to "where can I find US county outlines?".

A similar question was asked on the Natural Earth forum at http://www.naturalearthdata.com/forums/topic/u-s-county-shape-file/. It pointed to a location at http://nationalatlas.gov/mld/countyp.html, which has unfortunately had a dose of bit-rot. A quick Google suggests this can now be found at:

https://nationalmap.gov/small_scale/atlasftp.html?openChapters=chpbound#chpbound

I decided to download one of those county shapefiles:

https://prd-tnm.s3.amazonaws.com/StagedProducts/Small-scale/data/Boundaries/countyl010g_shp_nt00964.tar.gz

With that in place, I used cartopy's shapereader to get the geometries out, and created a custom feature that can then be added to the axes:

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


reader = shpreader.Reader('countyl010g.shp')

counties = list(reader.geometries())

COUNTIES = cfeature.ShapelyFeature(counties, ccrs.PlateCarree())

plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())

ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
ax.add_feature(COUNTIES, facecolor='none', edgecolor='gray')

ax.coastlines('50m')

ax.set_extent([-83, -65, 33, 44])
plt.show()

US counties

This is derived from the example at https://scitools.org.uk/cartopy/docs/v0.14/examples/feature_creation.html which constructs a NaturalEarthFeature, rather than a ShapelyFeature, but otherwise, the principle is pretty much the same.

Hope that is useful to you.

like image 63
pelson Avatar answered Oct 07 '22 09:10

pelson