Share @ LinkedIn Facebook  cartopy, maps
Cartopy - Basic Maps [Scatter Map, Bubble Map and Connection Map]

Cartopy - Basic Maps [Scatter Map, Bubble Map & Connection Map]

Table of Contents

Introduction

Many data analysis tasks nowadays involve working with geo spatial data and plot it on maps to get useful insights. For this purpose, we need a trustworthy and reliable plotting library that can be used for plotting maps from a different perspective. Python has a list of libraries (geopandas, bokeh, holoviews, plotly, bqplot, etc) which provides functionalities to plot different types of maps like scatter maps, choropleth map, bubble map, etc. We have already covered tutorials on creating maps using geopandas, holoviews and bqplot.

In this tutorial, we are going to introduce one more python map plotting library named cartopy. Canopy is a map plotting library in python which is based on matplotlib for plotting. It uses PROJ.4, numpy and shapely for handling data conversions between cartographic projection and handling shape files. Cartopy can be very useful to generate a high-quality static map chart that has high publication quality. If you want a high publication chart for your research paper and looking for plotting library with various projections available then cartopy is for you.

The main plus point of cartopy is its object-oriented projection system which lets us plot maps with points, lines polygons, etc using one projection and then can easily let us transfer it to another projection system.

As a part of this tutorial, we'll explore the API of cartopy in order to generate various maps. We'll be using various datasets available from kaggle for our map examples.

We'll start by loading necessary libraries for our plotting purpose.

In [1]:
import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
import cartopy

import sys
from itertools import combinations
from datetime import datetime
import random
import warnings

warnings.filterwarnings("ignore")

print("Python Version : ", sys.version)
print("Cartopy Version : ", cartopy.__version__)

%matplotlib inline
Python Version :  3.7.3 (default, Mar 27 2019, 22:11:17)
[GCC 7.3.0]
Cartopy Version :  0.17.0

1. First Simple Map

The generation of maps using cartopy can be done by following below mentioned common steps:

  • Create Matplotlib figure using plt.figure().
  • Add subplot to figure with projection attribute set as one of the projections available from cartopy.cars module.
  • Add features to map like Land, Ocean, Coastline, Borders, etc using cartopy.feature module.
  • Add your functionalities like points, lines, etc using common matplotlib API.
  • Show plot using plt.show().

Above mentioned steps are common steps when you want to generate world maps but it also provides other functionalities which can let you restrict map area to a particular region. We'll use the above-mentioned steps to plot a simple map.

Below we are plotting our first simple plot with projection PlateCarree. We have also added the COASTLINE feature so that we can see all continental borders.

In [ ]:
import cartopy.crs as crs
import cartopy.feature as cfeature

fig = plt.figure(figsize=(8,6))

ax = fig.add_subplot(1,1,1, projection=crs.PlateCarree())

ax.add_feature(cfeature.COASTLINE)

plt.show()

Cartopy

2. Trying Different Projections

We'll now try a few graph projections available in cartopy. The list of projections as available in crs module of the cartopy library.

Full Projections List: https://scitools.org.uk/cartopy/docs/latest/crs/projections.html#cartopy-projections

2.1 PlateCarree Projection

Below we are again trying the same PlateCarree projection which we had tried in the previous figure but we have also introduced natural earth look to the plot. We can add it by simply calling the stock_img() method on the axes object. We have also introduced a new way to add new coastline features to our plot.

In [ ]:
fig = plt.figure(figsize=(8,6))

ax = fig.add_subplot(1,1,1, projection=crs.PlateCarree())

ax.stock_img()
ax.coastlines()

plt.show()

Cartopy

2.2 AlbersEqualArea Projection

Below we have introduced another projection named AlbersEqualArea. We also have introduced grid lines on the plot as a light grey color. The graph generally keeps center at 0 longitudes and 0 latitudes which we can change by setting central_latitude and central_longitude parameter of AlbersEqualArea object.

In [ ]:
fig = plt.figure(figsize=(6,6))

ax = fig.add_subplot(1,1,1, projection=crs.AlbersEqualArea())

ax.stock_img()
ax.coastlines()
ax.gridlines()

plt.show()

Cartopy

2.3 Orthographic Projection

Below we have introduced another projection named Orthographic which has a shape of the globe. We also have passed its parameter central_latitude and central_longitude which centers glob location on that point which is the center of India.

In [ ]:
fig = plt.figure(figsize=(5,5))

ax = fig.add_subplot(1,1,1, projection=crs.Orthographic(central_latitude=20, central_longitude=80))

ax.stock_img()
ax.coastlines()

plt.show()

Cartopy

3. Adding More Features to Maps

We'll now explore further adding more futures to our map plots. Till we only had added coastline features to our maps.

3.1 Adding Coastline, Land, Lakes, Borders, Ocean, Rivers and States Features to Map

Below we are first creating a map with projection as Mercator. We are then adding features by calling method add_feature() on-axis object passing it feature object available from the cartopy.feature module. The list of features available through cartopy is available from the cartopy.feature module.

We are adding below feature to our map:

  • COASTLINE feature - Adds coastline around all continents. We can change a coastline color, style, etc.
  • LAND feature - Adds land on top of world map. We can change it's color and opacity as well.
  • LAKES feature - Adds big lakes of the world. We can change its configuration as well.
  • BORDERS feature - Adds country borders for the whole world. We can change border style and color, etc.
  • OCEAN feature - Adds ocean with a color that we can change by setting our own configuration as well.
  • RIVERS feature - Adds big rivers of the world.
  • STATES feature - Adds state boundary for the US.

There are few more classes available from cartopy which lets us add our own designed features but we have not discussed them here as they are not part of the scope of this tutorial.

In [ ]:
fig = plt.figure(figsize=(8,8))

ax = fig.add_subplot(1,1,1, projection=crs.Mercator())

ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.LAND, color="lightgrey", alpha=0.5)
ax.add_feature(cfeature.LAKES, color="lime")
ax.add_feature(cfeature.BORDERS, linestyle="--")
ax.add_feature(cfeature.OCEAN, color="skyblue", alpha=0.4)
ax.add_feature(cfeature.RIVERS, edgecolor="red")
ax.add_feature(cfeature.STATES)

plt.show()

Cartopy

3.2 Adding NightShade Feature to Map

Below we are adding one more feature to cartopy map called nightshade feature. We need to create a datetime object with datetime data and need to pass it to Nightshade object so that it can create nightshade on a map at this time. It'll create a shade effect on a map covering countries where it'll be night during that time.

In [ ]:
from cartopy.feature.nightshade import Nightshade

fig = plt.figure(figsize=(8, 5))

ax = fig.add_subplot(1, 1, 1, projection=crs.PlateCarree())

ax.stock_img()

current_time = datetime.now()
ax.add_feature(Nightshade(current_time, alpha=0.3))

ax.set_title('Night time shading for {}'.format(current_time))

plt.show()

Cartopy

4. Using Latitude & Longitude Box to Show a Particular Area of Map

Below we are explaining how to show only a particular part of the map instead of showing the whole world map. It happens many times that we only need to show a particular country, continent, and region of interest than showing the whole world map. We can do this in two ways:

  • We can call method named set_extent() on axes object passing it latitude and longitudes of the area to be covered.
  • We can also call set_xlim() and set_ylim() methods passing them longitude and latitude limits respectively.

We have used the below set_extent() method passing it projection as PlateCarree so that it projects box from that projection. We have tried to cover Asia and Europe continents on a map.

In [ ]:
fig = plt.figure(figsize=(6,6))

ax = fig.add_subplot(1,1,1, projection=crs.Mercator())

ax.stock_img()
ax.coastlines()
ax.add_feature(cfeature.BORDERS)

lat1, lon1, lat2, lon2 = -10, 180, 0, 80
ax.set_extent([lat1, lon1, lat2, lon2], crs=crs.PlateCarree())
#ax.set_xlim(65, 95)
#ax.set_ylim(5, 40)

plt.show()

Cartopy

We have now explained how to create various maps using cartopy. We'll now explain how to add further more data on maps like plots, lines, etc below. For that, we'll be using the Starbucks store locations dataset and brazil flight dataset.

5. Load Starbucks Locations Dataset

The first dataset that we'll load is Starbucks store locations dataset available from kaggle. We suggest that you download a dataset from kaggle to follow along.

This dataset has information about Starbucks store locations worldwide. It has also geospatial information (longitude & latitude) about each store location. Apart from geospatial information, it has data about store address, name, number, ownership type, country, state, city, postcode, etc.

In [9]:
starbucks_locations = pd.read_csv("datasets/starbucks_store_locations.csv")
starbucks_locations.head()
Out[9]:
Brand Store Number Store Name Ownership Type Street Address City State/Province Country Postcode Phone Number Timezone Longitude Latitude
0 Starbucks 47370-257954 Meritxell, 96 Licensed Av. Meritxell, 96 Andorra la Vella 7 AD AD500 376818720 GMT+1:00 Europe/Andorra 1.53 42.51
1 Starbucks 22331-212325 Ajman Drive Thru Licensed 1 Street 69, Al Jarf Ajman AJ AE NaN NaN GMT+04:00 Asia/Dubai 55.47 25.42
2 Starbucks 47089-256771 Dana Mall Licensed Sheikh Khalifa Bin Zayed St. Ajman AJ AE NaN NaN GMT+04:00 Asia/Dubai 55.47 25.39
3 Starbucks 22126-218024 Twofour 54 Licensed Al Salam Street Abu Dhabi AZ AE NaN NaN GMT+04:00 Asia/Dubai 54.38 24.48
4 Starbucks 17127-178586 Al Ain Tower Licensed Khaldiya Area, Abu Dhabi Island Abu Dhabi AZ AE NaN NaN GMT+04:00 Asia/Dubai 54.54 24.51

We'll be filtering dataset removing entries from all other countries present other than the United States.

In [10]:
starbucks_us = starbucks_locations[starbucks_locations["Country"] == "US"]
starbucks_us.head()
Out[10]:
Brand Store Number Store Name Ownership Type Street Address City State/Province Country Postcode Phone Number Timezone Longitude Latitude
11964 Starbucks 3513-125945 Safeway-Anchorage #1809 Licensed 5600 Debarr Rd Ste 9 Anchorage AK US 995042300 907-339-0900 GMT-09:00 America/Anchorage -149.78 61.21
11965 Starbucks 74352-84449 Safeway-Anchorage #2628 Licensed 1725 Abbott Rd Anchorage AK US 995073444 907-339-2800 GMT-09:00 America/Anchorage -149.84 61.14
11966 Starbucks 12449-152385 Safeway - Anchorage #1813 Licensed 1501 Huffman Rd Anchorage AK US 995153596 907-339-1300 GMT-09:00 America/Anchorage -149.85 61.11
11967 Starbucks 24936-233524 100th & C St - Anchorage Company Owned 320 W. 100th Ave, 100, Southgate Shopping Ctr ... Anchorage AK US 99515 (907) 227-9631 GMT-09:00 America/Anchorage -149.89 61.13
11968 Starbucks 8973-85630 Old Seward & Diamond Company Owned 1005 E Dimond Blvd Anchorage AK US 995152050 907-344-4160 GMT-09:00 America/Anchorage -149.86 61.14

6. Scatter Map of Starbucks Locations on Earth [Robinson Projection]

Below we are plotting scatter map of showing Starbucks locations on the whole Earth.

We have first created a map using a new projection named Robinson. We have then added scatter points on a map using commonly using plt.scatter() the method passing its longitude and latitude of locations. We have also passed it PlateCarree projection as a transform parameter value. We need to pass it so that it plots points on that projection first and then projects it on a map with Robinson projection. This is the convention that we'll be following in future maps as well.

In [ ]:
fig = plt.figure(figsize=(10,8))

ax = fig.add_subplot(1,1,1, projection=crs.Robinson())

ax.set_global()

ax.add_feature(cfeature.COASTLINE, edgecolor="tomato")
ax.add_feature(cfeature.BORDERS, edgecolor="tomato")
ax.gridlines()

plt.scatter(x=starbucks_locations.Longitude, y=starbucks_locations.Latitude,
            color="dodgerblue",
            s=1,
            alpha=0.5,
            transform=crs.PlateCarree()) ## Important

plt.show()

Cartopy

We can notice a high number of Starbucks locations in the whole US and Europe. There is also a sizable number of stores in Chine and Canada as well. All other locations of the world have very less or no Starbucks presence.

7. Scatter Bubble Map depicting Starbucks Store Count Per US States [PlateCarree Projection]

For creating a bubble map of Starbucks store count per US states, we need to prepare data by following the below-mentioned steps.

  • We'll first group dataset by State and then take average longitude, latitude for each state
  • We'll then group dataset by State and count store counts per state.
  • We'll then merge the dataframe created in the previous two steps.

Our final dataset will have an entry for each state along with its latitude, longitude, and Starbucks store count which we'll use for plotting bubble map.

In [12]:
## Below step will get us geospatial data for each state.
us_statewise_lat_lon = starbucks_us.groupby(by="State/Province").mean()[["Longitude", "Latitude"]]

## Below step will get us store count per state.
statewise_store_count = starbucks_us.groupby(by="State/Province").count()[["Store Number"]].rename(columns={"Store Number":"Count"})

## Below step will merge store counts dataframe with geospatial dataframe.
statewise_store_count = us_statewise_lat_lon.join(statewise_store_count).reset_index()

statewise_store_count.head()
Out[12]:
State/Province Longitude Latitude Count
0 AK -148.591224 61.376939 49
1 AL -86.855529 32.995294 85
2 AR -92.829091 35.240364 55
3 AZ -111.879877 33.379385 488
4 CA -119.505207 35.510578 2821

Once we are done with creating a dataset, we'll plot the bubble map exactly the same way as the previous map. We'll call plt.scatter() the method passing it latitude and longitude of each state of the US. We also have passed the size of the marker to be used for each step which is count of stores available in that state. This time we are also restricting graph by setting a limit on it using set_extent() so that it only displays US map. We also have added for states of the US.

In [ ]:
fig = plt.figure(figsize=(12,10))

ax = fig.add_subplot(1,1,1, projection=crs.PlateCarree())

ax.stock_img()
ax.coastlines()
ax.add_feature(cfeature.STATES)

ax.set_extent([-135, -66.5, 20, 55],
              crs=crs.PlateCarree()) ## Important


plt.scatter(x=statewise_store_count.Longitude, y=statewise_store_count.Latitude,
            color="orangered",
            s=statewise_store_count.Count,
            alpha=0.8,
            transform=crs.PlateCarree()) ## Important

plt.show()

Cartopy

We can notice from the above chart that there is a very high presence of Starbucks in California, Washington, Texas, Florida, Illinois, and New York.

8. Loading Brazil Flights Dataset

The second dataset that we'll be using is flight data for a number of flights from and to Brazil. The dataset is available at kaggle.

The dataset has information about flight time, type, airline name, source city, destination city, source country, destination country, source state, destination state, source latitude & longitude, and destination latitude & longitude.

We suggest that you download the dataset to follow along with us. Please make a note that the dataset is quite big and has nearly 2.5 Mn rows. If your computer has less memory and wants to follow along with tutorial then we suggest that you load the first few rows only using nrows parameter of read_csv().

In [14]:
brazil_flights = pd.read_csv("brazil_flights_data.csv", nrows=100000, encoding="latin")
brazil_flights.head()
Out[14]:
Voos Companhia.Aerea Codigo.Tipo.Linha Partida.Prevista Partida.Real Chegada.Prevista Chegada.Real Situacao.Voo Codigo.Justificativa Aeroporto.Origem ... Estado.Origem Pais.Origem Aeroporto.Destino Cidade.Destino Estado.Destino Pais.Destino LongDest LatDest LongOrig LatOrig
0 AAL - 203 AMERICAN AIRLINES INC Internacional 2016-01-30T08:58:00Z 2016-01-30T08:58:00Z 2016-01-30T10:35:00Z 2016-01-30T10:35:00Z Realizado NaN Afonso Pena ... PR Brasil Salgado Filho Porto Alegre RS Brasil -51.175381 -29.993473 -49.172481 -25.532713
1 AAL - 203 AMERICAN AIRLINES INC Internacional 2016-01-13T12:13:00Z 2016-01-13T12:13:00Z 2016-01-13T21:30:00Z 2016-01-13T21:30:00Z Realizado NaN Salgado Filho ... RS Brasil Miami Miami N/I Estados Unidos -80.287046 25.795865 -51.175381 -29.993473
2 AAL - 203 AMERICAN AIRLINES INC Internacional 2016-01-29T12:13:00Z 2016-01-29T12:13:00Z 2016-01-29T21:30:00Z 2016-01-29T21:30:00Z Realizado NaN Salgado Filho ... RS Brasil Miami Miami N/I Estados Unidos -80.287046 25.795865 -51.175381 -29.993473
3 AAL - 203 AMERICAN AIRLINES INC Internacional 2016-01-19T12:13:00Z 2016-01-18T12:03:00Z 2016-01-19T21:30:00Z 2016-01-18T20:41:00Z Realizado LIBERACAO SERV. TRAFEGO AEREO/ANTECIPACAO Salgado Filho ... RS Brasil Miami Miami N/I Estados Unidos -80.287046 25.795865 -51.175381 -29.993473
4 AAL - 203 AMERICAN AIRLINES INC Internacional 2016-01-30T12:13:00Z 2016-01-30T12:13:00Z 2016-01-30T21:30:00Z 2016-01-30T21:30:00Z Realizado NaN Salgado Filho ... RS Brasil Miami Miami N/I Estados Unidos -80.287046 25.795865 -51.175381 -29.993473

5 rows × 21 columns

9. Connection Map depicting Flights from Brazil to Other Countries [Mollweide Projection]

To plot a connection map depicting flights from Brazil to other countries, we need to first prepare the dataset by following a list of the below steps.

  • We'll filter and remove rows where the destination country is Brazil, as well as source and destination country, is Brazil. This way our dataset will have only entries from Brazil to other countries.
  • We'll then group dataset by origin & destination country and calculate a count of flights from Brazil to that country.
  • We'll then group dataset by origin & destination country, take mean of all columns and then filter the dataset to have only source & destination geospatial data.
  • We'll then merge dataset created in the previous two steps to create a final dataset.

Our final dataset will have information about flight count from source to destination countries as well as geospatial data for each country.

In [15]:
## Removing entry where source and destination are not same and destination country is not Brazil. 
## This will remove domestic entries as well as flights from other countries to brazil.
## It'll keep entries of flights from brazil to other countries.
brazil_flights = brazil_flights[(brazil_flights["Pais.Origem"]!=brazil_flights["Pais.Destino"]) & (brazil_flights["Pais.Destino"]!="Brasil")]

## This will get us count of flight from Brazil to each country.
brazil_to_other_country_count = brazil_flights.groupby(by=["Pais.Origem", "Pais.Destino"]).count()[["Voos"]].rename(columns={"Voos":"Count"})

## This will get us mean of all columns grouped according source and destination countries.
brazil_to_other_country_lat_lon = brazil_flights.groupby(by=["Pais.Origem", "Pais.Destino"]).mean()

## This will filter and remove columns other that geospatial data columns
brazil_to_other_country_lat_lon = brazil_to_other_country_lat_lon[["LongDest","LatDest","LongOrig","LatOrig"]]

## This will merge flight count dataframe with geospatial data dataframe.
brazil_to_other_country_flights = brazil_to_other_country_count.join(brazil_to_other_country_lat_lon).reset_index()

brazil_to_other_country_flights.head()
Out[15]:
Pais.Origem Pais.Destino Count LongDest LatDest LongOrig LatOrig
0 Brasil Africa Do Sul 45 28.241146 -26.136673 -46.478126 -23.434553
1 Brasil Alemanha 147 9.239697 49.682746 -44.737826 -21.822414
2 Brasil Angola 31 13.234862 -8.847951 -45.124154 -23.174074
3 Brasil Antilhas Holandesas 31 -69.571231 12.370281 -56.107170 -8.958052
4 Brasil Argentina 1796 -58.827371 -34.644030 -45.827534 -23.188851

Our last map is a connection map where we'll be displaying flight movement from Brazil to various countries of the earth using a line. We also show a thick line where there are more flights going on and a thin line where fewer flights are going.

We have used Mollweide projection this time for the creation of a world map. We are then looping through each flight entry taking the source and destination latitude, longitude to draw line between them using conventional plt.plot() method. We have used Geodetic projection for plotting lines so that it shows a curvy line instead of a straight line from source to destination. We are also plotting names of destination countries where flights are landing from Brazil to highlight them.

In [ ]:
fig = plt.figure(figsize=(15,8))

ax = fig.add_subplot(1,1,1, projection=crs.Mollweide())

ax.stock_img()
ax.set_global()
ax.coastlines()
ax.add_feature(cfeature.STATES)

## Looping thorugh each entry of flights
for idx in brazil_to_other_country_flights.index:

    ## Getting source and destination latitude as well as longitude
    lat1, lon1 = brazil_to_other_country_flights.loc[idx]["LatOrig"], brazil_to_other_country_flights.loc[idx]["LongOrig"]
    lat2, lon2 = brazil_to_other_country_flights.loc[idx]["LatDest"], brazil_to_other_country_flights.loc[idx]["LongDest"]

    ## Getting flights count
    num_flights = brazil_to_other_country_flights.loc[idx]["Count"]

    ## Plotting line from source country to destination country.
    plt.plot([lon1, lon2], [lat1, lat2],
             linewidth=2 if num_flights > 50 else 0.5,
             color="lime",
             transform=crs.Geodetic()) ## Important

    ## Plotting destination country name on map.
    plt.text(lon2, lat2,
             brazil_to_other_country_flights.loc[idx]["Pais.Destino"],
             color="red", fontsize=10,
             horizontalalignment="center", verticalalignment="center",
             transform=crs.Geodetic()) ## Important

plt.show()

Cartopy

We can sizable number of flight traveling from Brazil to Southern American countries, US, Canada, and European countries. There are some travels to African countries as well but none traveling to Asia.

This ends our small tutorial on covering the basics of cartopy. Please feel free to let us know your views in the comments section.

References



Sunny Solanki  Sunny Solanki