Why Use Interactive Maps

Interactive Maps are useful for earth data science because they:

• Clearly convey complex information.
• Are more engaging for viewers than static maps.
• Can be seamlessly integrated into Jupyter Notebooks.

There are two great Python packages for creating interactive maps: folium and mapboxgl. Both of these packages are build on top off the JavaScript library called leaflet.js.

This lesson will focus on folium, which has been around longer than mapboxgl and thus, is well-documented by the Python community.

One major difference between the two packages is that mapboxgl requires a MapBox Access Token. If you are interested in exploring mapboxgl on your own, note that the MapBox Access token is free to use, but does require making an account with MapBox.

You can find more information on the Github page for this package.

What is an API?

An API (or application programming interface) is an interface that opens a computer-based system to external requests and simplifies certain tasks, such as extracting subsets of data from a large repository or database.

For example, web-based APIs allow you to access data available using web-based interfaces that are separate from the API that you are accessing. Thus, web APIs are a way to avoid the extraneous visual interfaces that you do not need and get the desired data into the tool that you want to use.

Often, you access data from web-based APIs using a URL that contains sets of parameters that specifies the type and particular subset of data that you are interested in. You will learn more about using APIs later in this course.

For this lesson, you simply need to know that the basemaps that you will access to create your interactive maps come from APIs that are provided by various organizations such as OpenStreetMap, MapBox, Stamen, Google, etc.

# Import necessary packages
import os 
import folium
from folium import plugins
import rioxarray as rxr
import earthpy as et
import earthpy.spatial as es

# Import data from EarthPy
data = et.data.get_data('colorado-flood')

# Set working directory to earth-analytics
os.chdir(os.path.join(et.io.HOME, 'earth-analytics', 'data'))

Simple Basemap

You can make an interactive map with folium using just one line of code!

You can use the Map() function from folium and providing a latitude and longitude to center the map. The map is created using the default basemap from OpenStreetMap.

# Create a map using the Map() function and the coordinates for Boulder, CO
m = folium.Map(location=[40.0150, -105.2705])

# Display the map
m

Change Basemap

You can change the basemap for the map by providing a value for the tiles parameter of the Map() function.

There are many different options including Stamen Terrain, Stamen Toner and cartodbpositron. More details and basemaps names available on the Folium Documentation for the Map() function.

# Create a map using Stamen Terrain as the basemap
m = folium.Map(location=[40.0150, -105.2705],
               tiles = 'Stamen Terrain')

# Display map
m

Add Markers

You can also add markers to label specific points on top of a folium basemap, such as the coordinates that are being used to center the map. You can even add a pop-up label for the marker that is activated when you click on it.

# Create a map using Stamen Terrain, centered on Boulder, CO
m = folium.Map(location=[40.0150, -105.2705], 
               tiles = 'Stamen Terrain')

# Add marker for Boulder, CO
folium.Marker(
    location=[40.009515, -105.242714], # coordinates for the marker (Earth Lab at CU Boulder)
    popup='Earth Lab at CU Boulder', # pop-up label for the marker
    icon=folium.Icon()
).add_to(m)

# Display m
m

Raster Overlay on Interactive Map

You can also overlay rasters on folium basemaps.

The default coordinate system and projection for web-based basemaps is WGS84 Web Mercator. To overlay data on web-based basemaps, the overlay data needs to be in the WGS84 coordinate system (see this link for more information on this coordinate system).

Thus, to overlay a raster on a basemap, you first need to project the raster to WGS84 (EPSG 4326).

Project Raster

You can use the rioxarray package, which you imported as rxr, to project a raster. In this example, you will use a raster for a post-flood digital terrain model (DTM) in the Northwest Boulder area: post_DTM.tif.

Prepare Raster for Plotting

Folium is powerful, but you still need to prepare the raster before you can plot it.

The first thing you must do is replace all the nan values that rioxarray used to mask out areas with no data. These values don’t work with folium. In this case, you can replace those values with the minimum value in the image, but you can replace them with whatever values you see fit to best visualize your data.

The next thing we must do is scale the data to only contain values from 0 to 255. This makes the data of type uint8, which is necessary to plot properly in folium. Luckily, earthpy has a function to do just this! You can put the array in es.bytescale and it will take the values and scale them from 0 to 255 for you.

# Create a variable for destination coordinate system 
dst_crs = 'EPSG:4326' 

# Path to raster
in_path = os.path.join("colorado-flood", 
                       "spatial", 
                       "boulder-leehill-rd", 
                       "post-flood", 
                       "lidar", 
                       "post_DTM.tif")

# Open the raster in rioxarray
img = rxr.open_rasterio(in_path, masked=True)

# Reproject the raster to be the correct crs
img = img.rio.reproject(dst_crs)

# Replace all null values with the minimum value in the array
img_plot = img.where(~img.isnull(), img.min())

# Scale the array from 0 to 255
scaled_img = es.bytescale(img_plot.values[0])

Overlay Raster

Now that the raster is in the correct coordinate system (WGS84), you can overlay it on the basemp using the add_child() function and specifying the image (e.g. img) and setting an opacity and bounding box, if desired.

# Create a map using Stamen Terrain, centered on study area with set zoom level
m = folium.Map(location=[40.06, -105.30],
               tiles='Stamen Terrain', 
               zoom_start = 13)

map_bounds = [[40.05577828237005, -105.32837712340124], 
              [40.073923431943214, -105.28139535136515]]
# Overlay raster called img using add_child() function (opacity and bounding box set)
m.add_child(folium.raster_layers.ImageOverlay(scaled_img, 
                                              opacity=.9,
                                              bounds=map_bounds))

# Display map 
m