Introduction to Risk Mapping

Combining factors to create simple potential risk area maps.

import folium
import ee
import geemap
import geopandas as gpd
ee.Authenticate()

# Write your project ID here, in quotes
ee.Initialize(project = "emerge-lessons")

def add_ee_layer(self, ee_image_object, vis_params, name):
    """Adds a method for displaying Earth Engine image tiles to folium map."""
    map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)
    folium.raster_layers.TileLayer(
        tiles=map_id_dict['tile_fetcher'].url_format,
        attr='Map Data &copy; <a href="https://earthengine.google.com/">Google Earth Engine</a>',
        name=name,
        overlay=True,
        control=True
    ).add_to(self)

folium.Map.add_ee_layer = add_ee_layer

start_date = '2025-06-01'
end_date = '2025-06-30'

def mask_s2_clouds(image):
  """Masks clouds in a Sentinel-2 image using the QA band.
  Args:
      image (ee.Image): A Sentinel-2 image.
  Returns:
      ee.Image: A cloud-masked Sentinel-2 image.
  """
  qa = image.select('QA60')

  # Bits 10 and 11 are clouds and cirrus, respectively.
  cloud_bit_mask = 1 << 10
  cirrus_bit_mask = 1 << 11

  # Both flags should be set to zero, indicating clear conditions.
  mask = (
      qa.bitwiseAnd(cloud_bit_mask)
      .eq(0)
      .And(qa.bitwiseAnd(cirrus_bit_mask).eq(0))
  )

  return image.updateMask(mask).divide(10000)

Define the region of interest (over Central Florida).

point = ee.Geometry.Point(-81.660044, 28.473813)
region = point.buffer(distance=100000)

Load in data from Google Earth Engine, similar to lesson 1 from this chapter.

lst = (
    ee.ImageCollection('MODIS/061/MOD11A1')
      .filterBounds(region)
      .filterDate(start_date, end_date)
      .select('LST_Day_1km')
      .median()
      .clip(region)
)

rgb = (
    ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
      .filterBounds(region)
      .filterDate(start_date, end_date)
      .map(mask_s2_clouds)
      .median()
      .clip(region)
)
ndvi = rgb.normalizedDifference(['B8', 'B4']).rename('NDVI')

rain = (
    ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY')
      .filterBounds(region)
      .filterDate(start_date, end_date)
      .select('precipitation')
      .sum()
      .clip(region)
)

Visualize the data on a map (select the layer to show using the tool in the top right).

map = folium.Map(location=[28.263363, -83.497652], tiles="Cartodb dark_matter", zoom_start=7)

map.add_ee_layer(lst,
                {'min': 13000.0, 'max': 16500.0,
                'palette': ['040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
                '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
                '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
                'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
                'ff0000', 'de0101', 'c21301', 'a71001', '911003']},
                 "LST")
map.add_ee_layer(ndvi, {'min': -1, 'max': 1, 'palette': ['blue', 'white', 'green']}, 'NDVI')
map.add_ee_layer(rain, {'min': 0, 'max': 200, 'palette': ['white', 'blue']}, 'Precipitation')

folium.LayerControl(collapsed = False).add_to(map)
display(map)

Make this Notebook Trusted to load map: File -> Trust Notebook

Combine the layers together.

# Define a common projection (Sentinel-2)
s2_proj = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED') \
    .first().select('B1').projection()

# Resample LST and Precipitation to Sentinel-2 resolution
lst_resampled = lst.resample('bilinear').reproject(crs=s2_proj)
rain_resampled = rain.resample('bilinear').reproject(crs=s2_proj)

# Combine layers to one image
combined = ndvi.addBands(lst).addBands(rain)

# Calculate means
mean_lst = lst.reduceRegion(reducer=ee.Reducer.mean(), geometry=region, scale=1000,
                              maxPixels=1e12).get('LST_Day_1km').getInfo()
mean_rain = rain.reduceRegion(reducer=ee.Reducer.mean(), geometry=region, scale=5566,
                              maxPixels=1e12).get('precipitation').getInfo()

# Define risk classification function
def classify_risk(image):
    ndvi = image.select('NDVI')
    lst = image.select('LST_Day_1km')
    rain = image.select('precipitation')

    low_risk = ndvi.lt(0).And(lst.lt(mean_lst)).And(rain.lt(mean_rain))
    med_risk = ndvi.lte(0.3).Or(lst.eq(mean_lst)).Or(rain.eq(mean_rain))
    high_risk = ndvi.gt(0.3).And(lst.gt(mean_lst)).And(rain.gt(mean_rain))

    # Assign values: 1 = Low, 2 = Medium, 3 = High
    risk = low_risk.multiply(1).add(med_risk.multiply(2)).add(high_risk.multiply(3)).rename('Risk_Level')
    return image.addBands(risk)

# Apply classification
classified = classify_risk(combined)

# Refresh map
map = folium.Map(location=[28.263363, -83.497652], tiles="Cartodb dark_matter", zoom_start=7)

# Add risk layer
map.add_ee_layer(classified.select('Risk_Level'), {'min': 1, 'max': 3, 'palette': ['green', 'yellow', 'red']}, 'Risk Level')
display(map)

Make this Notebook Trusted to load map: File -> Trust Notebook

This is a very simple map showing low, medium, or high risk based on just a few variables (NDVI, LST, Precipitation). Previous research have found other variables that relate to mosquito risk, where it is important to consider local environmental conditions.