Skip navigation
NASA Logo, National Aeronautics and Space Administration
Currently Being Moderated

Reading Shapefile Dataset with Cartopy

VERSION 2  Click to view document history
Created on: Oct 14, 2015 1:18 PM by Jules Kouatchou - Last Modified:  Oct 14, 2015 3:11 PM by Jules Kouatchou

What is Cartopy?

  1. Package for drawing maps for for data analysis and visualization
  2. Relies on PROJ.4, numpy and shapely libraries
  3. Has a simple and intuitive drawing interface to Matplotlib

 

What Does Cartopy Provide?

  1. Facilities to transform coordinates to different map projections
  2. Matplotlib is used to plot contours, images, vectors, lines or points in the transformed coordinates.
  3. Shorelines, river and political boundary datasets.
  4. Facilities for reading shapefiles.

 

What is Shapefile?

The shapefile format:

  1. Is a digital vector storage format for storing geometric location and associated attribute information.
  2. Geographic features in a shapefile can be represented by points, lines, or polygons (areas)
  3. Is non-topological. It does not maintain spatial relationship information such as connectivity, adjacency, and area definition.
  4. Was introduced with ArcView GIS version 2 in the early 1990s.

 

Content of Shapefile Files

Every shapefile data set includes at least three files:

  • .shp: The main file that contains the primary geographic reference data and records of various shape types included, such as points, polygons, or multipatches.
  • .dbf: The dBase file that stores attributes for each shape. It alows quicker access to the spatial features of the data.
  • .shx: Organize the records of a shapefile for reference.

 

 

Cartopy and the Natural Earth Dataset

Cartopy provides an interface for access to frequently used data such as the GSHHS dataset and from the NaturalEarthData website. These interfaces allow the user to define the data programmatically, and if the data does not exist on disk, it will be retrieved from the appropriate source (normally by downloading the data from the internet).

 

To acquire the countries dataset from Natural Earth, we may use:

import cartopy.io.shapereader as shpreader

 

shpfilename = shpreader.natural_earth(resolution='110m', category='cultural', name='admin_0_countries')

 

We can use the Reader function to get the first country in the shapefile:

 

reader = shpreader.Reader(shpfilename)

countries = reader.records()

country = next(countries)

 

Accessing Cartopy

To have access to cartopy on discover, you need to load the modules:

 

other/comp/gcc-4.6.3-sp3

lib/mkl-15.0.2.164

other/SIVO-PyD/spd_1.23.0_gcc-4.6.3-sp3_mkl-15.0.2.164

 

Example 1: Map the Globe and Color the United States

 

#!/usr/bin/env python

 

import matplotlib.pyplot as plt

import cartopy

import cartopy.io.shapereader as shpreader

import cartopy.crs as ccrs

 

# Select the map projection

#----------------------

ax = plt.axes(projection=ccrs.PlateCarree())

ax.add_feature(cartopy.feature.OCEAN)

 

# Select the area of interest

#-----------------------

ax.set_extent([-150, 60, -25, 60])

 

# Read the Natural Earth shapefile dataset

#----------------------------------

shpfilename = shpreader.natural_earth(resolution='110m', category='cultural', name='admin_0_countries')

reader = shpreader.Reader(shpfilename)

countries = reader.records()

 

for country in countries:

    if country.attributes['adm0_a3'] == 'USA':

        ax.add_geometries(country.geometry, ccrs.PlateCarree(), facecolor=(0, 0, 1),

                          label=country.attributes['adm0_a3'])

    else:

        ax.add_geometries(country.geometry, ccrs.PlateCarree(),

                          facecolor=(0, 1, 0), label=country.attributes['adm0_a3'])

 

plt.show()

 

Example 2: Select the country Cameroon and color each of its administrative region with a different color

#!/usr/bin/env python

 

from matplotlib.colors import cnames

import cartopy.crs as ccrs

import cartopy.feature as cfeature

from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

from cartopy.io import shapereader

 

# Read the Natural Earth shapefile dataset

#----------------------------------

kw = dict(resolution='10m', category='cultural', name='admin_1_states_provinces')

states_shp = shapereader.natural_earth(**kw)

shp = shapereader.Reader(states_shp)

 

# Select the map projection

#----------------------

subplot_kw = dict(projection=ccrs.PlateCarree())

 

fig, ax = plt.subplots(figsize=(7, 11), subplot_kw=subplot_kw)

 

# Select the area that includes Cameroon

#----------------------------------

ax.set_extent([7.5, 17.5, 1.25, 15])

 

ax.background_patch.set_visible(False)

ax.outline_patch.set_visible(False)

 

 

# Get from Matplotlib a list of colors

#------------------------------

colors = list(cnames.keys())

len_colors = len(colors)

 

k = 0

for record, state in zip(shp.records(), shp.geometries()):

    if record.attributes['admin'] == 'Cameroon':

        if k+1 == len_colors:

            k = 0

        else:

            k += 1

        facecolor = colors[k]

    else:

        facecolor = 'LightGray'

    ax.add_geometries([state], ccrs.PlateCarree(), facecolor=facecolor, edgecolor='black')

 

plt.show()

 

Comments (0)
USAGov logo NASA Logo - nasa.gov