Contextily 1: Intro

Overview:¶

Read in the NYS Mesonet sites table with GeoPandas
Plot the NYS Mesonet sites on a map
Use the contextily package to overlay an image as a basemap, served via a remote tile server

Prerequisites¶

Concepts	Importance	Notes
(Geo)Pandas	Necessary
Cartopy	Necessary

Time to learn: 10 minutes

Imports¶

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import geopandas as gp
import pandas as pd
from metpy.plots import USCOUNTIES
import contextily as ctx
import matplotlib.pyplot as plt

Read in the NYS Mesonet sites table with Geopandas¶

First, open the NYSM sites table using Pandas, as we’ve done before.

nysm_sites = pd.read_csv('/spare11/atm533/data/nysm_sites.csv')

nysm_sites.head()

In the Shapefiles Intro notebook, we used Geopandas to read in a shapefile as a Dataframe. Geopandas is not limited to shapefiles; it can also georeference existing Pandas dataframes that contain relevant columns or rows. The NYSM dataframe has latitude and longitude columns; we will pass these in to the GeoDataFrame method which creates the geo-referenced dataframe.

Note: A GeoDataFrame requires a column named geometry. Besides polygons and multipolygons, another valid geometry is Points. We thus create a Geometry series using Geopandas' points_from_xy method. This method accepts two series... in this case, the longitude and latitude columns from the NYSM sites dataframe.

gdf = gp.GeoDataFrame(nysm_sites,geometry=gp.points_from_xy(nysm_sites.lon,nysm_sites.lat))

Examine the `GeoDataFrame` and its geometry column¶

gdf

As we did in the first Shapefiles notebook, let’s take a look at the geometry attribute of one of the sites.

gdf.loc[gdf.stid == 'VOOR'].geometry

112    POINT (-73.97562 42.65242)
Name: geometry, dtype: geometry

It’s a POINT ... whose x- and y- coordinates are the corresponding longitude and latitude in degrees.

Also as in the previous notebook, we can use Geopandas `plot` method to get a quick visualization.¶

gdf.plot(figsize=(12,9));

Make a plot using Matplotlib and Cartopy.¶

To plot the points corresponding to NYS Mesonet sites, we use the same Geopandas plot method, but need to specify the axes. In this case, we first define a set of GeoAxes, as we always do when using Cartopy, and then specify the name of that GeoAxes object as the value for the ax argument.

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

res='10m'
plt.title('NYS Mesonet Sites')
gdf.plot(ax=ax,color='tab:red')

ax.coastlines(resolution=res)
ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)
ax.add_feature (cfeature.STATES.with_scale(res))
ax.add_feature(USCOUNTIES,edgecolor='grey');
ax.set_extent([-80, -72,40.5,45.2], ccrs.PlateCarree())

Add a background GeoTIFF raster image.¶

We’ll use the contextily package’s add_basemap method for that. It will request GeoTIFFs from a remote server that not only provides a source of background maps, but also makes them available at various zoom levels. Depending on how we define the extent of our map, one or more tiles will be downloaded from the remote image server. The default basemap comes from Openstreetmap.

Note: The default coordinate reference system (crs) in Contextily is WebMercator (EPSG 3857). We can either transform the remote images to our desired projection, or set up our map using WebMercator.

We specify a zoom level here, although if you do not specify one, Contextily will determine an appropriate default value.

Tip: The larger the zoom value, the greater detail in the background image (but larger values will take longer to download and render!).

fig = plt.figure(figsize=(15,10))
ax = plt.axes(projection=ccrs.epsg('3857'))

res='10m'
plt.title('NYS Mesonet Sites')
gdf.plot(ax=ax,color='tab:red',transform=ccrs.PlateCarree())
ax.coastlines(resolution=res)
ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)
ax.add_feature (cfeature.STATES.with_scale(res))
ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );
ax.set_extent([-80, -71.4,40.5,45.2], ccrs.PlateCarree())
ctx.add_basemap(ax,zoom=9,transform=ccrs.epsg('3857'))

Try a lower zoom value.

fig = plt.figure(figsize=(15,10))
ax = plt.axes(projection=ccrs.epsg('3857'))

res='10m'
plt.title('NYS Mesonet Sites')
gdf.plot(ax=ax,color='tab:red',transform=ccrs.PlateCarree())
ax.coastlines(resolution=res)
ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)
ax.add_feature (cfeature.STATES.with_scale(res))
ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );

ax.set_extent([-80, -71.4,40.5,45.2], ccrs.PlateCarree())
ctx.add_basemap(ax,zoom=7,transform=ccrs.epsg('3857'))

Plot a map whose extent covers just the Greater Capital District.¶

Do not specify a zoom level. Let contextily’s add_ basemap figure determine the optimum one.

fig = plt.figure(figsize=(15,10))
ax = plt.axes(projection=ccrs.epsg('3857'))

res='10m'
plt.title('NYS Mesonet Sites')
gdf.plot(ax=ax,color='tab:red',transform=ccrs.PlateCarree())
ax.coastlines(resolution=res)

ax.add_feature (cfeature.LAKES.with_scale(res), alpha = 0.5)
ax.add_feature (cfeature.STATES.with_scale(res))
ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1 );

ax.set_extent([-74.4, -73.4,42.3,43.1], ccrs.PlateCarree())
ctx.add_basemap(ax)

Question: How might you add the NYSM site ID's to the plot?

Summary¶

Besides shapefiles, Geopandas can geo-reference Pandas dataframes.
The contextily package allows for the plotting of remotely-served basemaps.
These basemap servers support tiling; thus, zooming in or out will reveal greater or less features.