5. Global map of water indices with Geemap

This notebook is optimized to run in Google Colab. An interactive demo of this example is deployed here.

In [1]:

#!pip install geemap spyndex fast-dash

#!pip install geemap spyndex fast-dash

In [1]:

import ee
import geemap.foliumap as geemap
import spyndex

import datetime
import pandas as pd
from fast_dash import FastDash, Fastify, html, dbc, dmc, dash

ee.Authenticate()
ee.Initialize()

import ee import geemap.foliumap as geemap import spyndex import datetime import pandas as pd from fast_dash import FastDash, Fastify, html, dbc, dmc, dash ee.Authenticate() ee.Initialize()

Read GAUL data values¶

In [3]:

data = pd.read_csv("assets/gaul_data.csv")

data = pd.read_csv("assets/gaul_data.csv")

In [4]:

data.head()

data.head()

Out[4]:

	ADM0_NAME	ADM1_NAME	ADM2_NAME
0	United Republic of Tanzania	Arusha	Karatu
1	United Republic of Tanzania	Iringa	Iringa Rural
2	United Republic of Tanzania	Iringa	Kilolo
3	United Republic of Tanzania	Manyara	Mbulu
4	United Republic of Tanzania	Pwani	Kisarawe

Get a list of water indices¶

In [5]:

columns = ["application_domain", "contributor", "date_of_addition", "long_name", "platforms", "reference", "short_name"]
track_index = []

for index in spyndex.indices.keys():
    i = spyndex.indices[index]
    track_index.append([i.application_domain,
                        i.contributor,
                        i.date_of_addition,
                        i.long_name,
                        i.platforms,
                        i.reference,
                        i.short_name])
                        
index_df = pd.DataFrame(track_index, columns=columns)

water_indices = index_df[(index_df.application_domain == "water") & 
                    (index_df.platforms.apply(lambda x: "Landsat-OLI" in x))]

# Add label column
water_indices["label"] = water_indices.apply(lambda row: f"{row.short_name} ({row.long_name})", axis=1)

columns = ["application_domain", "contributor", "date_of_addition", "long_name", "platforms", "reference", "short_name"] track_index = [] for index in spyndex.indices.keys(): i = spyndex.indices[index] track_index.append([i.application_domain, i.contributor, i.date_of_addition, i.long_name, i.platforms, i.reference, i.short_name]) index_df = pd.DataFrame(track_index, columns=columns) water_indices = index_df[(index_df.application_domain == "water") & (index_df.platforms.apply(lambda x: "Landsat-OLI" in x))] # Add label column water_indices["label"] = water_indices.apply(lambda row: f"{row.short_name} ({row.long_name})", axis=1)

/tmp/ipykernel_368/3195300281.py:20: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

In [6]:

water_indices.head()

water_indices.head()

Out[6]:

	application_domain	contributor	date_of_addition	long_name	platforms	reference	short_name	label
2	water	https://github.com/davemlz	2022-09-22	Augmented Normalized Difference Water Index	[Sentinel-2, Landsat-OLI, Landsat-TM, Landsat-...	https://doi.org/10.1016/j.envsoft.2021.105030	ANDWI	ANDWI (Augmented Normalized Difference Water I...
8	water	https://github.com/davemlz	2021-09-18	Automated Water Extraction Index	[Sentinel-2, Landsat-OLI, Landsat-TM, Landsat-...	https://doi.org/10.1016/j.rse.2013.08.029	AWEInsh	AWEInsh (Automated Water Extraction Index)
9	water	https://github.com/davemlz	2021-09-18	Automated Water Extraction Index with Shadows ...	[Sentinel-2, Landsat-OLI, Landsat-TM, Landsat-...	https://doi.org/10.1016/j.rse.2013.08.029	AWEIsh	AWEIsh (Automated Water Extraction Index with ...
69	water	https://github.com/davemlz	2022-04-20	Land Surface Water Index	[Sentinel-2, Landsat-OLI, Landsat-TM, Landsat-...	https://doi.org/10.1016/j.rse.2003.11.008	LSWI	LSWI (Land Surface Water Index)
71	water	https://github.com/davemlz	2022-01-17	Multi-Band Water Index	[Sentinel-2, Landsat-OLI, Landsat-TM, Landsat-...	https://doi.org/10.1016/j.jag.2018.01.018	MBWI	MBWI (Multi-Band Water Index)

Define utility functions¶

In [7]:

def get_index_function(index_name, image_collection):
    # Get the index function from Spyndex.
    return spyndex.computeIndex(index=index_name,
                     params={"B": image_collection.select("B2"),
                            "G": image_collection.select("B3"),
                            "R": image_collection.select("B4"),
                            "N": image_collection.select("B5"),
                            "S1": image_collection.select("B6"),
                            "S2": image_collection.select("B7"),
                            "T1": image_collection.select("B10"),
                            "T2": image_collection.select("B11"),
                            "gamma": 1,
                            "alpha": 1})

def calculate_mean_index(image, region):
    mean_value = image.reduceRegion(
        reducer=ee.Reducer.mean(),
        geometry=region,
        scale=30,  # Adjust the scale according to your dataset and accuracy requirements
        maxPixels=1e9
    )
    return mean_value.getInfo()

# Define a color palette manually.
palette = [
    '000033',  # Very dark blue
    '000066',  # Dark blue
    '000099',  # Medium dark blue
    '0000CC',  # Moderate blue
    '0000FF',  # Blue
    '3399FF',  # Lighter blue
    '66CCFF',  # Light blue
    '99CCFF',  # Very light blue
    'CCE6FF'   # Extremely light blue
]

def get_index_function(index_name, image_collection): # Get the index function from Spyndex. return spyndex.computeIndex(index=index_name, params={"B": image_collection.select("B2"), "G": image_collection.select("B3"), "R": image_collection.select("B4"), "N": image_collection.select("B5"), "S1": image_collection.select("B6"), "S2": image_collection.select("B7"), "T1": image_collection.select("B10"), "T2": image_collection.select("B11"), "gamma": 1, "alpha": 1}) def calculate_mean_index(image, region): mean_value = image.reduceRegion( reducer=ee.Reducer.mean(), geometry=region, scale=30, # Adjust the scale according to your dataset and accuracy requirements maxPixels=1e9 ) return mean_value.getInfo() # Define a color palette manually. palette = [ '000033', # Very dark blue '000066', # Dark blue '000099', # Medium dark blue '0000CC', # Moderate blue '0000FF', # Blue '3399FF', # Lighter blue '66CCFF', # Light blue '99CCFF', # Very light blue 'CCE6FF' # Extremely light blue ]

Define the main function and build Fast Dash app¶

In [8]:

country_component = dmc.Select(data=data.ADM0_NAME.unique().tolist(), 
                               value="United States of America",
                               label="Country of interest",
                               searchable=True)

state_component = dmc.Select(label="State, provience or an equivalent administrative unit",
                             value="New York",
                             searchable=True)

county_component = dmc.Select(label="County, district or an equivalent administrative unit",
                              value="Albany",
                             searchable=True)

index_component = dmc.Select(data=water_indices.label.unique().tolist(),
                             value="NDVIMNDWI (NDVI-MNDWI Model)",
                             searchable=True,
                             clearable=True)

# Define a function that takes start_date, end_date, country, city, and index as arguments.
def water_spectral_indices(country: country_component,
                        state_or_province: state_component = None,
                        county_or_district: county_component = None,
                        water_index: index_component = None,
                        minimum_index_value: int = -1,
                        maximum_index_value: int = 1) -> html.Iframe(height="100%"):
    
    """
    Visualize the selected water spectral index. Select your geography from the inputs to get started. Map on the left displays the median index value \
    for the years 2013 - 2015 and the map on the right displays median values for the years 2021-2023. Compare the two maps to understand how these values \
    changed over the years.
    
    Learn more about the Spyndex Python library that enables this dynamic index computation [here](https://github.com/awesome-spectral-indices/spyndex).
    
    :param country: The country of interest.
    :type country: str
    
    :param state_or_province: State, provience or an equivalent administrative unit.
    :type state_or_province: str
    
    :param county_or_district: County, district or an equivalent administrative unit.
    :type county_or_district: str
    
    :param water_index: Name of the water index.
    :type water_index: str
    
    :return: HTML of the leafmap object.
    :rtype: str
    """
    
    if not state_or_province:
        raise Exception("Please select a state")
        
    if not county_or_district:
        raise Exception("Please select a county")
    
    # Load the Global Administrative Unit Layers (GAUL) dataset, which includes administrative boundaries.
    admin_boundaries = ee.FeatureCollection('FAO/GAUL/2015/level2')

    # Filter the GAUL dataset for the specified city and country.
    geometry = admin_boundaries \
        .filter(ee.Filter.eq('ADM0_NAME', country)) \
        .filter(ee.Filter.eq('ADM1_NAME', state_or_province)) \
        .filter(ee.Filter.eq('ADM2_NAME', county_or_district)) \
        .geometry()
        
    # Load a satellite image collection (e.g., Landsat 8 Surface Reflectance).
    left_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \
        .filterBounds(geometry) \
        .filterDate("2013-01-01", "2015-12-31") \
        .median() \
        .clip(geometry)
    
    right_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \
        .filterBounds(geometry) \
        .filterDate("2021-01-01", "2023-12-31") \
        .median() \
        .clip(geometry)
    
    # Get index name
    index_name = water_indices[water_indices.label == water_index].short_name.iloc[0]
    
    index_function_left = get_index_function(index_name, left_image)
    index_function_right = get_index_function(index_name, right_image)
    
    # Generate legend and visualization parameters
    # Define visualization parameters using the custom palette.
    vis_params = {
        'min': minimum_index_value,
        'max': maximum_index_value,
        'palette': palette
    }

    # Calculate the range interval for each color.
    num_colors = len(palette)
    interval = 2.0 / num_colors

    # Create the legend dictionary with range values as keys.
    legend_dict = {}
    for i, color in enumerate(palette):
        # Define the range for each color.
        range_min = round(-1 + i * interval, 2)
        range_max = round(-1 + (i + 1) * interval, 2)
        # Use the range as the key and the color as the value.
        legend_dict[f'{range_min}-{range_max}'] = color
    
    # Visualize the spectral index using Geemap.
    Map = geemap.Map(basemap="CartoDB.Positron")
    Map.centerObject(geometry, zoom=11)
    
    left_layer = geemap.ee_tile_layer(index_function_left, vis_params, f"{index_name} 2013-2015")
    right_layer = geemap.ee_tile_layer(index_function_right, vis_params, f"{index_name} 2021-2023")
    Map.split_map(left_layer, right_layer)

    # Add a legend to the map.
    Map.add_legend(title="Legend", legend_dict=legend_dict, position='bottomright')
    
    comparison_map = Map.to_html()

    return comparison_map

country_component = dmc.Select(data=data.ADM0_NAME.unique().tolist(), value="United States of America", label="Country of interest", searchable=True) state_component = dmc.Select(label="State, provience or an equivalent administrative unit", value="New York", searchable=True) county_component = dmc.Select(label="County, district or an equivalent administrative unit", value="Albany", searchable=True) index_component = dmc.Select(data=water_indices.label.unique().tolist(), value="NDVIMNDWI (NDVI-MNDWI Model)", searchable=True, clearable=True) # Define a function that takes start_date, end_date, country, city, and index as arguments. def water_spectral_indices(country: country_component, state_or_province: state_component = None, county_or_district: county_component = None, water_index: index_component = None, minimum_index_value: int = -1, maximum_index_value: int = 1) -> html.Iframe(height="100%"): """ Visualize the selected water spectral index. Select your geography from the inputs to get started. Map on the left displays the median index value \ for the years 2013 - 2015 and the map on the right displays median values for the years 2021-2023. Compare the two maps to understand how these values \ changed over the years. Learn more about the Spyndex Python library that enables this dynamic index computation [here](https://github.com/awesome-spectral-indices/spyndex). :param country: The country of interest. :type country: str :param state_or_province: State, provience or an equivalent administrative unit. :type state_or_province: str :param county_or_district: County, district or an equivalent administrative unit. :type county_or_district: str :param water_index: Name of the water index. :type water_index: str :return: HTML of the leafmap object. :rtype: str """ if not state_or_province: raise Exception("Please select a state") if not county_or_district: raise Exception("Please select a county") # Load the Global Administrative Unit Layers (GAUL) dataset, which includes administrative boundaries. admin_boundaries = ee.FeatureCollection('FAO/GAUL/2015/level2') # Filter the GAUL dataset for the specified city and country. geometry = admin_boundaries \ .filter(ee.Filter.eq('ADM0_NAME', country)) \ .filter(ee.Filter.eq('ADM1_NAME', state_or_province)) \ .filter(ee.Filter.eq('ADM2_NAME', county_or_district)) \ .geometry() # Load a satellite image collection (e.g., Landsat 8 Surface Reflectance). left_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \ .filterBounds(geometry) \ .filterDate("2013-01-01", "2015-12-31") \ .median() \ .clip(geometry) right_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \ .filterBounds(geometry) \ .filterDate("2021-01-01", "2023-12-31") \ .median() \ .clip(geometry) # Get index name index_name = water_indices[water_indices.label == water_index].short_name.iloc[0] index_function_left = get_index_function(index_name, left_image) index_function_right = get_index_function(index_name, right_image) # Generate legend and visualization parameters # Define visualization parameters using the custom palette. vis_params = { 'min': minimum_index_value, 'max': maximum_index_value, 'palette': palette } # Calculate the range interval for each color. num_colors = len(palette) interval = 2.0 / num_colors # Create the legend dictionary with range values as keys. legend_dict = {} for i, color in enumerate(palette): # Define the range for each color. range_min = round(-1 + i * interval, 2) range_max = round(-1 + (i + 1) * interval, 2) # Use the range as the key and the color as the value. legend_dict[f'{range_min}-{range_max}'] = color # Visualize the spectral index using Geemap. Map = geemap.Map(basemap="CartoDB.Positron") Map.centerObject(geometry, zoom=11) left_layer = geemap.ee_tile_layer(index_function_left, vis_params, f"{index_name} 2013-2015") right_layer = geemap.ee_tile_layer(index_function_right, vis_params, f"{index_name} 2021-2023") Map.split_map(left_layer, right_layer) # Add a legend to the map. Map.add_legend(title="Legend", legend_dict=legend_dict, position='bottomright') comparison_map = Map.to_html() return comparison_map

In [9]:

app = FastDash(water_spectral_indices, port=8001, theme="cosmo", mode="external")

app = FastDash(water_spectral_indices, port=8001, theme="cosmo", mode="external")

WARNING:root:Parsing function docstring is still an experimental feature. To reduce uncertainty, consider setting `about` to `False`.

Additional app functionality¶

In [10]:

# Keep only the states present in the selected country
@app.app.callback(dash.Output("state_or_province", "data"),
                 dash.Input("country", "value"))
def filter_states(country):
    return sorted(data[data.ADM0_NAME == country].ADM1_NAME.unique().tolist())

# Keep only the counties present in the selected country and state
@app.app.callback(dash.Output("county_or_district", "data"),
                 dash.Input("country", "value"),
                 dash.Input("state_or_province", "value"))
def filter_counties(country, state):
    filtered_data = data.copy()
    if country:
        filtered_data = filtered_data[filtered_data.ADM0_NAME == country]
        
    if state:
        filtered_data = filtered_data[filtered_data.ADM1_NAME == state]
        
    return sorted(filtered_data.ADM2_NAME.unique().tolist())

# Keep only the states present in the selected country @app.app.callback(dash.Output("state_or_province", "data"), dash.Input("country", "value")) def filter_states(country): return sorted(data[data.ADM0_NAME == country].ADM1_NAME.unique().tolist()) # Keep only the counties present in the selected country and state @app.app.callback(dash.Output("county_or_district", "data"), dash.Input("country", "value"), dash.Input("state_or_province", "value")) def filter_counties(country, state): filtered_data = data.copy() if country: filtered_data = filtered_data[filtered_data.ADM0_NAME == country] if state: filtered_data = filtered_data[filtered_data.ADM1_NAME == state] return sorted(filtered_data.ADM2_NAME.unique().tolist())

In [11]:

app.run()

app.run()

Dash app running on http://127.0.0.1:8001/