None

Note

This tutorial was generated from an IPython notebook that can be downloaded here.

Working with geopandas (shapefiles)

regionmask includes support for regions defined as geopandas GeoDataFrame. These are often shapefiles, which can be opened in the formats .zip or .shp with geopandas.read_file(url_or_path).

There are two possibilities:

  1. Directly create a mask from a geopandas GeoDataFrame or GeoSeries using mask_geopandas or mask_3D_geopandas.

  2. Convert a GeoDataFrame to a Regions object (regionmask’s internal data container) using from_geopandas.

As always, start with the imports:

import cartopy.crs as ccrs
import geopandas as gp
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pooch

import regionmask

regionmask.__version__

'0.6.2'

Opening an example shapefile

The U.S. Geological Survey (USGS) offers a shapefile containing the outlines of continens [1]. We use the library pooch to locally cache the file:

file = pooch.retrieve(
    "https://pubs.usgs.gov/of/2006/1187/basemaps/continents/continents.zip", None
)

continents = gp.read_file("zip://" + file)

display(continents)

Downloading data from 'https://pubs.usgs.gov/of/2006/1187/basemaps/continents/continents.zip' to file '/home/docs/.cache/pooch/7dd514faeaa71efe73294dece9245e99-continents.zip'.
SHA256 hash of downloaded file: af0ba524a62ad31deee92a9700fc572088c2b93a39ba66f320677dd8dacaaaaf
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.

       CONTINENT                                           geometry
0           Asia  MULTIPOLYGON (((93.27554 80.26361, 93.31304 80...
1  North America  MULTIPOLYGON (((-25.28167 71.39166, -25.32889 ...
2         Europe  MULTIPOLYGON (((58.06138 81.68776, 57.98055 81...
3         Africa  MULTIPOLYGON (((0.69465 5.77337, 0.66667 5.803...
4  South America  MULTIPOLYGON (((-81.71306 12.49028, -81.72014 ...
5        Oceania  MULTIPOLYGON (((-177.39334 28.18416, -177.3958...
6      Australia  MULTIPOLYGON (((142.27997 -10.26556, 142.21053...
7     Antarctica  MULTIPOLYGON (((51.80305 -46.45667, 51.72139 -...

Create a mask from a GeoDataFrame

mask_geopandas and mask_3D_geopandas allow to directly create a mask from a GeoDataFrame or GeoSeries:

lon = np.arange(-180, 180)
lat = np.arange(-90, 90)

mask = regionmask.mask_geopandas(continents, lon, lat)

Let’s plot the new mask:

f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))
mask.plot(
    ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False,
)

ax.coastlines(color="0.1");

<cartopy.mpl.feature_artist.FeatureArtist at 0x7f57f5405750>
../_images/geopandas_9_1.png

Similarly a 3D boolean mask can be created from a GeoDataFrame:

mask_3D = regionmask.mask_3D_geopandas(continents, lon, lat)

and plotted:

from matplotlib import colors as mplc

cmap1 = mplc.ListedColormap(["none", "#9ecae1"])

f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))

mask_3D.sel(region=0).plot(
    ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap=cmap1,
)

ax.coastlines(color="0.1");

<cartopy.mpl.feature_artist.FeatureArtist at 0x7f57fd8fe590>
../_images/geopandas_13_1.png

2. Convert GeoDataFrame to a Regions object

Creating a Regions object with regionmask.from_geopandas requires a GeoDataFrame:

continents_regions = regionmask.from_geopandas(continents)
continents_regions

<regionmask.Regions>
Name:     unnamed

Regions:
0 r0 Region0
1 r1 Region1
2 r2 Region2
3 r3 Region3
4 r4 Region4
5 r5 Region5
6 r6 Region6
7 r7 Region7

[8 regions]

This creates default names ("Region0", …, "RegionN") and abbreviations ("r0", …, "rN").

However, it is often advantageous to use columns of the GeoDataFrame as names and abbrevs. If no column with abbreviations is available, you can use abbrevs='_from_name', which creates unique abbreviations using the names column.

continents_regions = regionmask.from_geopandas(
    continents, names="CONTINENT", abbrevs="_from_name", name="continent"
)
continents_regions

<regionmask.Regions>
Name:     continent

Regions:
0    Asi          Asia
1 NorAme North America
2    Eur        Europe
3    Afr        Africa
4 SouAme South America
5    Oce       Oceania
6    Aus     Australia
7    Ant    Antarctica

[8 regions]

As usual the newly created Regions object can be plotted on a world map:

continents_regions.plot(label="name", coastlines=False);

<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x7f57ea1e4890>
../_images/geopandas_19_1.png

And to create mask a mask for arbitrary latitude/ longitude grids:

lon = np.arange(0, 360)
lat = np.arange(-90, 90)

mask = continents_regions.mask(lon, lat)

which can then be plotted

f, ax = plt.subplots(subplot_kw=dict(projection=ccrs.PlateCarree()))
h = mask.plot(
    ax=ax,
    transform=ccrs.PlateCarree(),
    cmap="Reds",
    add_colorbar=False,
    levels=np.arange(-0.5, 8),
)

cbar = plt.colorbar(h, shrink=0.625, pad=0.025, aspect=12)
cbar.set_ticks(np.arange(8))
cbar.set_ticklabels(continents_regions.names)

ax.coastlines(color="0.2")

continents_regions.plot_regions(add_label=False);

<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x7f57ea3105d0>
../_images/geopandas_23_1.png

References

[1] Environmental Systems Research , Inc. (ESRI), 20020401, World Continents: ESRI Data & Maps 2002, Environmental Systems Research Institute, Inc. (ESRI), Redlands, California, USA.