EarthCube 2023

Interactive Geospatial Analysis and Data Visualization with Leafmap

This notebook provides an introduction to interactive geospatial analysis and data visualization with leafmap. It is designed for the notebook demo at the EarthCube 2023.

Installation¶

Uncomment and run the following cell to install necessary packages for this notebook. Restart the kernel/runtime after installation.

In [ ]:

# %pip install leafmap[raster] geopandas rasterstats

# %pip install leafmap[raster] geopandas rasterstats

Import libraries¶

In [ ]:

import os
import leafmap
import geopandas as gpd

import os import leafmap import geopandas as gpd

Create interactive maps¶

Specify the map center [lat, lon], zoom level, and height.

In [ ]:

m = leafmap.Map(center=[40, -100], zoom=4, height="600px")
m

m = leafmap.Map(center=[40, -100], zoom=4, height="600px") m

Add basemaps¶

Add OpenTopoMap, USGS 3DEP Elevation, and USGS Hydrography basemaps.

In [ ]:

m = leafmap.Map()
m.add_basemap("OpenTopoMap")
m.add_basemap("USGS 3DEP Elevation")
m.add_basemap("USGS Hydrography")
m.add_layer_manager()
m

m = leafmap.Map() m.add_basemap("OpenTopoMap") m.add_basemap("USGS 3DEP Elevation") m.add_basemap("USGS Hydrography") m.add_layer_manager() m

Add NLCD land cover map and legend.

In [ ]:

m = leafmap.Map(center=[40, -100], zoom=4)
m.add_basemap("HYBRID")
m.add_basemap("NLCD 2019 CONUS Land Cover")
m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type")
m

m = leafmap.Map(center=[40, -100], zoom=4) m.add_basemap("HYBRID") m.add_basemap("NLCD 2019 CONUS Land Cover") m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type") m

Add WMS layers.

In [ ]:

m = leafmap.Map(center=[40, -100], zoom=4)
m.add_basemap("Esri.WorldImagery")
url = "https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2019_Land_Cover_L48/wms?"
m.add_wms_layer(
    url,
    layers="NLCD_2019_Land_Cover_L48",
    name="NLCD 2019 CONUS Land Cover",
    format="image/png",
    transparent=True,
)
m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type")
m

m = leafmap.Map(center=[40, -100], zoom=4) m.add_basemap("Esri.WorldImagery") url = "https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2019_Land_Cover_L48/wms?" m.add_wms_layer( url, layers="NLCD_2019_Land_Cover_L48", name="NLCD 2019 CONUS Land Cover", format="image/png", transparent=True, ) m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type") m

Change basemaps interactively.

In [ ]:

m = leafmap.Map()
m.add_basemap_gui()
m

m = leafmap.Map() m.add_basemap_gui() m

Visualize raster datasets¶

Cloud Optimized GeoTIFF¶

A Cloud Optimized GeoTIFF (COG) is a regular GeoTIFF file that is optimized for serving on an HTTP file server, with an internal organization that enables more efficient workflows on a cloud environment. This is achieved by allowing clients to issue HTTP GET requests to ask for only the parts of a file that they need. For more information about COG, please visit https://www.cogeo.org.

In [ ]:

url = "https://open.gishub.org/data/raster/srtm90.tif"

url = "https://open.gishub.org/data/raster/srtm90.tif"

In [ ]:

leafmap.cog_stats(url)

leafmap.cog_stats(url)

In [ ]:

m = leafmap.Map()
m.add_cog_layer(url, name="Remote COG", colormap_name="terrain")
m.add_colormap(vmin=0, vmax=4000, cmap="terrain", label="Elevation (m)")
m.add_inspector_gui()
m

m = leafmap.Map() m.add_cog_layer(url, name="Remote COG", colormap_name="terrain") m.add_colormap(vmin=0, vmax=4000, cmap="terrain", label="Elevation (m)") m.add_inspector_gui() m

Download the COG file.

In [ ]:

leafmap.download_file(url, "srtm90.tif")

leafmap.download_file(url, "srtm90.tif")

In [ ]:

m = leafmap.Map()
m.add_raster("srtm90.tif", colormap="terrain", layer_name="Local COG")
m.add_colormap(vmin=0, vmax=4000, cmap="terrain", label="Elevation (m)")
m

m = leafmap.Map() m.add_raster("srtm90.tif", colormap="terrain", layer_name="Local COG") m.add_colormap(vmin=0, vmax=4000, cmap="terrain", label="Elevation (m)") m

SpatioTemporal Asset Catalog¶

The SpatioTemporal Asset Catalog (STAC) specification provides a common language to describe a range of geospatial information so that it can more easily be indexed and discovered. A SpatioTemporal Asset is any file that represents information about the earth captured in a certain space and time. STAC aims to enable that next generation of geospatial search engines, while also supporting web best practices so geospatial information is more easily surfaced in traditional search engines. More information about STAC can be found at the STAC website. The STAC Index website is a one-stop-shop for discovering STAC catalogs, collections, APIs, software and tools. In this example, we will use the STAC SPOT Orthoimages of Canada.

In [ ]:

url = "https://canada-spot-ortho.s3.amazonaws.com/canada_spot_orthoimages/canada_spot5_orthoimages/S5_2007/S5_11055_6057_20070622/S5_11055_6057_20070622.json"

url = "https://canada-spot-ortho.s3.amazonaws.com/canada_spot_orthoimages/canada_spot5_orthoimages/S5_2007/S5_11055_6057_20070622/S5_11055_6057_20070622.json"

In [ ]:

leafmap.stac_bands(url)

leafmap.stac_bands(url)

In [ ]:

m = leafmap.Map()
m.add_stac_layer(url, bands=["pan"], name="Panchromatic")
m.add_stac_layer(url, bands=["B3", "B2", "B1"], name="False color")
m.add_layer_manager()
m

m = leafmap.Map() m.add_stac_layer(url, bands=["pan"], name="Panchromatic") m.add_stac_layer(url, bands=["B3", "B2", "B1"], name="False color") m.add_layer_manager() m

Visualize vector datasets¶

GeoJSON¶

In [ ]:

m = leafmap.Map()
url = "https://open.gishub.org/data/world/continents.geojson"
m.add_geojson(url, layer_name="Continents")
m

m = leafmap.Map() url = "https://open.gishub.org/data/world/continents.geojson" m.add_geojson(url, layer_name="Continents") m

Shapefile¶

In [ ]:

m = leafmap.Map()
url = "https://open.gishub.org/data/world/countries.zip"
m.add_shp(url, layer_name="Countries")
m

m = leafmap.Map() url = "https://open.gishub.org/data/world/countries.zip" m.add_shp(url, layer_name="Countries") m

Other vector formats¶

In [ ]:

m = leafmap.Map()
url = "https://github.com/opengeos/leafmap/raw/master/examples/data/cable_geo.geojson"
m.add_vector(url, layer_name="Cable lines")
m

m = leafmap.Map() url = "https://github.com/opengeos/leafmap/raw/master/examples/data/cable_geo.geojson" m.add_vector(url, layer_name="Cable lines") m

Vector styling¶

In [ ]:

import leafmap

m = leafmap.Map()
m.add_basemap("CartoDB.DarkMatter")
url = "https://github.com/opengeos/leafmap/raw/master/examples/data/cable_geo.geojson"
callback = lambda feat: {"color": feat["properties"]["color"]}
m.add_vector(url, layer_name="Cable lines", style_callback=callback)
m

import leafmap m = leafmap.Map() m.add_basemap("CartoDB.DarkMatter") url = "https://github.com/opengeos/leafmap/raw/master/examples/data/cable_geo.geojson" callback = lambda feat: {"color": feat["properties"]["color"]} m.add_vector(url, layer_name="Cable lines", style_callback=callback) m

XY coordinates¶

In [ ]:

url = "https://github.com/opengeos/leafmap/raw/master/examples/data/us_cities.csv"
leafmap.csv_to_df(url).head()

url = "https://github.com/opengeos/leafmap/raw/master/examples/data/us_cities.csv" leafmap.csv_to_df(url).head()

In [ ]:

m = leafmap.Map(center=[40, -100], zoom=4)
m.add_points_from_xy(url, x="longitude", y="latitude")
m

m = leafmap.Map(center=[40, -100], zoom=4) m.add_points_from_xy(url, x="longitude", y="latitude") m

Split Map¶

In [ ]:

m = leafmap.Map(center=[36.1, -114.9], zoom=10)
m.add_basemap("HYBRID")
m.split_map(
    left_layer="NLCD 2001 CONUS Land Cover",
    right_layer="NLCD 2019 CONUS Land Cover",
    left_label="2001",
    right_label="2019",
)
m.add_layer_control()
m

m = leafmap.Map(center=[36.1, -114.9], zoom=10) m.add_basemap("HYBRID") m.split_map( left_layer="NLCD 2001 CONUS Land Cover", right_layer="NLCD 2019 CONUS Land Cover", left_label="2001", right_label="2019", ) m.add_layer_control() m

In [ ]:

m = leafmap.Map()
before = (
    "https://github.com/opengeos/data/releases/download/raster/Libya-2023-07-01.tif"
)
after = "https://github.com/opengeos/data/releases/download/raster/Libya-2023-09-13.tif"
m.split_map(before, after, left_label="Before", right_label="After")
m

m = leafmap.Map() before = ( "https://github.com/opengeos/data/releases/download/raster/Libya-2023-07-01.tif" ) after = "https://github.com/opengeos/data/releases/download/raster/Libya-2023-09-13.tif" m.split_map(before, after, left_label="Before", right_label="After") m

OpenAerialMap¶

OpenAerialMap (OAM) provides openly licensed satellite and unmanned aerial vehicle (UAV) imagery. Currently, it has over 12,400+ images around the globe. You can search and visualize OAM imagery interactively. You can also download images automatically with one line of code.

In [ ]:

bbox = [-79.6344, -0.9063, -77.3383, 1.0436]
start_date = "2016-04-01"
end_date = "2016-04-30"
gdf = leafmap.oam_search(
    bbox=bbox, start_date=start_date, end_date=end_date, return_gdf=True
)
print(f"Found {len(gdf)} images")

bbox = [-79.6344, -0.9063, -77.3383, 1.0436] start_date = "2016-04-01" end_date = "2016-04-30" gdf = leafmap.oam_search( bbox=bbox, start_date=start_date, end_date=end_date, return_gdf=True ) print(f"Found {len(gdf)} images")

In [ ]:

gdf.head()

gdf.head()

The tile URLs are stored in the tms column of the GeoDataFrame.

In [ ]:

tiles = gdf["tms"].tolist()
tiles[:5]

tiles = gdf["tms"].tolist() tiles[:5]

The image sources (downloadable URLs) are stored in the uuid column of the GeoDataFrame.

In [ ]:

images = gdf["uuid"].tolist()
images[:5]

images = gdf["uuid"].tolist() images[:5]

Download all images using the download_files() function.

In [ ]:

leafmap.download_files(images[:2])

leafmap.download_files(images[:2])

Add the image footprints to the map.

In [ ]:

import leafmap

m = leafmap.Map()
m.add_gdf(gdf, layer_name="Footprints")
m.add_cog_layer(images[0], nodata=0, name="OpenAerialMap")
# m.add_tile_layer(tiles[0], attribution='OpenAerialMap', name='OpenAerialMap')
m

import leafmap m = leafmap.Map() m.add_gdf(gdf, layer_name="Footprints") m.add_cog_layer(images[0], nodata=0, name="OpenAerialMap") # m.add_tile_layer(tiles[0], attribution='OpenAerialMap', name='OpenAerialMap') m

Search OAM imagery interactively using the interactive GUI.

In [ ]:

m = leafmap.Map(center=[4.7955, -75.6899], zoom=15)
m.add_oam_gui()
m

m = leafmap.Map(center=[4.7955, -75.6899], zoom=15) m.add_oam_gui() m

AWS Open Data¶

The AWS Open Data Program hosts a lots of open and public datasets on AWS, including Landsat 8, Sentinel-2, NAIP, and many more. Check out https://github.com/opengeos#data-catalogs for a list of open and public datasets on AWS.

In [ ]:

url = "https://earth-search.aws.element84.com/v1/"
collection = "sentinel-2-l2a"
time_range = "2021-12-01/2021-12-31"
bbox = [-122.2751, 47.5469, -121.9613, 47.7458]

url = "https://earth-search.aws.element84.com/v1/" collection = "sentinel-2-l2a" time_range = "2021-12-01/2021-12-31" bbox = [-122.2751, 47.5469, -121.9613, 47.7458]

In [ ]:

search = leafmap.stac_search(
    url=url,
    max_items=10,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    query={"eo:cloud_cover": {"lt": 10}},
    sortby=[{"field": "properties.eo:cloud_cover", "direction": "asc"}],
)
search

search = leafmap.stac_search( url=url, max_items=10, collections=[collection], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 10}}, sortby=[{"field": "properties.eo:cloud_cover", "direction": "asc"}], ) search

In [ ]:

search = leafmap.stac_search(
    url=url,
    max_items=10,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    get_collection=True,
)
search

search = leafmap.stac_search( url=url, max_items=10, collections=[collection], bbox=bbox, datetime=time_range, get_collection=True, ) search

In [ ]:

search = leafmap.stac_search(
    url=url,
    max_items=10,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    get_gdf=True,
)
search.head()

search = leafmap.stac_search( url=url, max_items=10, collections=[collection], bbox=bbox, datetime=time_range, get_gdf=True, ) search.head()

In [ ]:

search = leafmap.stac_search(
    url=url,
    max_items=10,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    get_assets=True,
)
search

search = leafmap.stac_search( url=url, max_items=10, collections=[collection], bbox=bbox, datetime=time_range, get_assets=True, ) search

In [ ]:

search = leafmap.stac_search(
    url=url,
    max_items=10,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    get_info=True,
)
search

search = leafmap.stac_search( url=url, max_items=10, collections=[collection], bbox=bbox, datetime=time_range, get_info=True, ) search

In [ ]:

search = leafmap.stac_search(
    url=url,
    max_items=10,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    get_links=True,
)
search

search = leafmap.stac_search( url=url, max_items=10, collections=[collection], bbox=bbox, datetime=time_range, get_links=True, ) search

In [ ]:

m = leafmap.Map(center=[37.7517, -122.4433], zoom=8)
m.add_stac_gui()
m

m = leafmap.Map(center=[37.7517, -122.4433], zoom=8) m.add_stac_gui() m

In [ ]:

# m.stac_gdf

# m.stac_gdf

In [ ]:

# m.stac_dict

# m.stac_dict

In [ ]:

# m.stac_item

# m.stac_item

Split raster into tiles¶

In [ ]:

m = leafmap.Map()
url = "https://open.gishub.org/data/raster/srtm90.tif"
m.add_cog_layer(url, name="Remote COG")
m

m = leafmap.Map() url = "https://open.gishub.org/data/raster/srtm90.tif" m.add_cog_layer(url, name="Remote COG") m

In [ ]:

leafmap.split_raster(url, out_dir="tiles", tile_size=(1000, 1000), overlap=0)

leafmap.split_raster(url, out_dir="tiles", tile_size=(1000, 1000), overlap=0)

In [ ]:

tiles = leafmap.find_files("tiles", ext=".tif")
tiles

tiles = leafmap.find_files("tiles", ext=".tif") tiles

In [ ]:

for tile in tiles[:6]:
    name = os.path.basename(tile).replace(".tif", "")
    m.add_raster(
        tile, cmap="terrain", vmin=0, vmax=4000, layer_name=name, zoom_to_layer=False
    )

for tile in tiles[:6]: name = os.path.basename(tile).replace(".tif", "") m.add_raster( tile, cmap="terrain", vmin=0, vmax=4000, layer_name=name, zoom_to_layer=False )

In [ ]:

m.add_layer_manager()
m

m.add_layer_manager() m

Merge tiles into a single raster¶

In [ ]:

leafmap.merge_rasters("tiles", output="merged.tif", input_pattern="*.tif")

leafmap.merge_rasters("tiles", output="merged.tif", input_pattern="*.tif")

In [ ]:

m = leafmap.Map()
m.add_raster("merged.tif", colormap="terrain", layer_name="Merged raster")
m

m = leafmap.Map() m.add_raster("merged.tif", colormap="terrain", layer_name="Merged raster") m

Zonal statistics¶

In [ ]:

dsm = "https://open.gishub.org/data/elevation/dsm.tif"
hag = "https://open.gishub.org/data/elevation/hag.tif"
buildings = "https://open.gishub.org/data/elevation/buildings.geojson"

dsm = "https://open.gishub.org/data/elevation/dsm.tif" hag = "https://open.gishub.org/data/elevation/hag.tif" buildings = "https://open.gishub.org/data/elevation/buildings.geojson"

In [ ]:

m = leafmap.Map()
m.add_cog_layer(dsm, name="DSM", palette="terrain")
m.add_cog_layer(hag, name="Height Above Ground", palette="magma")
m.add_geojson(buildings, layer_name="Buildings")
m

m = leafmap.Map() m.add_cog_layer(dsm, name="DSM", palette="terrain") m.add_cog_layer(hag, name="Height Above Ground", palette="magma") m.add_geojson(buildings, layer_name="Buildings") m

In [ ]:

gdf = gpd.read_file(buildings)
len(gdf)

gdf = gpd.read_file(buildings) len(gdf)

In [ ]:

gdf.head()

gdf.head()

In [ ]:

types = ["min", "max", "mean", "std", "count"]
gdf = leafmap.zonal_stats(gdf, hag, stats=types, gdf_out=True)
gdf

types = ["min", "max", "mean", "std", "count"] gdf = leafmap.zonal_stats(gdf, hag, stats=types, gdf_out=True) gdf

In [ ]:

m = leafmap.Map()
m.add_gdf(gdf, layer_name="Zonal Stats")
m

m = leafmap.Map() m.add_gdf(gdf, layer_name="Zonal Stats") m