First we will use cartopy's shapereader to download (and cache) states shapefile with 50 meters resolution from the NaturalEarth.

In [2]:

from cartopy.io import shapereader


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

states_shp = shapereader.natural_earth(**kw)
shp = shapereader.Reader(states_shp)

Now let's create a list of a few states that we want to highlight.

(The unicode_literals is needed to compare the Unicode names with the decoded names from latin-1. Argh!!! One day all of python will be Python3 and this ordeal will be over.)

In [3]:

from __future__ import unicode_literals

states = ('Minas Gerais', 'Mato Grosso', 'Goiás',
          'Bahia', 'Rio Grande do Sul', 'São Paulo')

Let's plot it all with cartopy.

In [4]:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt


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

fig, ax = plt.subplots(figsize=(7, 11),
                       subplot_kw=subplot_kw)
ax.set_extent([-82, -32, -45, 10])

ax.background_patch.set_visible(False)
ax.outline_patch.set_visible(False)

for record, state in zip(shp.records(), shp.geometries()):
    name = record.attributes['name'].decode('latin-1')
    if name in states:
        facecolor = 'DarkOrange'
    else:
        facecolor = 'LightGray'
    ax.add_geometries([state], ccrs.PlateCarree(),
                      facecolor=facecolor, edgecolor='black')

/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/matplotlib/artist.py:221: MatplotlibDeprecationWarning: This has been deprecated in mpl 1.5, please use the
axes property.  A removal date has not been set.
  warnings.warn(_get_axes_msg, mplDeprecation, stacklevel=1)

Cartopy's maps are great, but they are not interactive. We can fix that by plotting the same data over a folium Map instance.

The trick is to save the shapefile as a GeoJSON and plot it with folium's .geo_json.

In [5]:

import os
import shapefile
from json import dumps


def shape2json(fname, outfile="states.json", country='Brazil'):
    reader = shapefile.Reader(fname)
    fields = reader.fields[1:]
    field_names = [field[0] for field in fields]

    data = []
    for sr in reader.shapeRecords():
        atr = dict(zip(field_names, sr.record))
        geom = sr.shape.__geo_interface__
        if country in sr.record[field_names.index('admin')]:
            name = sr.record[field_names.index('name')].decode('latin-1')
            if name in states:
                data.append(dict(type="Feature", geometry=geom, properties=atr))
            
    keys = ['abbrev', 'name', 'name_alt']
    for b in data:
        for key in keys:
            b['properties'][key] = b['properties'][key].decode('latin-1')

    with open(outfile, "w") as geojson:
        geojson.write(dumps({"type": "FeatureCollection",
                             "features": data}, indent=2) + "\n")


shape = '50m_admin_1_states_provinces.shp'
cartopy_cache = '.local/share/cartopy/shapefiles/natural_earth/cultural/'
fname = os.path.join(os.path.expanduser('~'), cartopy_cache, shape)

shape2json(fname, outfile="states.json", country='Brazil')

In [6]:

import folium
import numpy as np
from IPython.display import IFrame


def inline_map(m, width=650, height=500):
    """Takes a folium instance and embed HTML."""
    m._build_map()
    srcdoc = m.HTML.replace('"', '&quot;')
    embed = HTML('<iframe srcdoc="{}" '
                 'style="width: {}px; height: {}px; '
                 'border: none"></iframe>'.format(srcdoc, width, height))
    return embed

bbox = [-82, -32, -45, 10]

lon_center, lat_center = np.array(bbox).reshape(2, 2).mean(axis=0)
mapa = folium.Map(width=650, height=500, zoom_start=4,
                  location=[-15, -50])

mapa.geo_json(geo_path='states.json', fill_color='none', line_color='Orange')
inline_map(mapa)

Out[6]:

That was easy! Check here to learn more about GeoJSONs plotting with folium.

In [7]:

HTML(html)

Out[7]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

Plotting shapefiles with cartopy and folium

Comments