Part 2. GLOBE Analysis with Environmental Data#
If you have not used Google Earth Engine before, follow this tutorial to set up your Earth Engine project. This will let you load environmental data for your GLOBE analysis!
Introduction to the Tools#
Python (a popular programming knowledge for analyzing data, automating tasks, creating software, training machine learning models, and more) to analyze GLOBE Mosquito Habitat Mapper and Land Cover data.
Google Colab (a free, cloud-based coding environment) to run Python on the browser. No setup or downloads necessary!
Google Earth Engine (a cloud-based platform that includes satellite imagery and other geospatial datasets) to load environmental data.
To run the code:
Each block of code is called a cell. To run a cell, hover over it and click the arrow in the top left of the cell, or click inside of the cell and press Shift + Enter.
Note: When you run a block of code for the first time, Google Colab will say Warning: This notebook was not authored by Google. Please click Run Anyway.
# Write your project ID for Google Earth Engine below
projectID = "Write Your Project ID"
# For example,
# projectID = "emerge-lessons-23148"
state_name = "Your State Name"
county_name = "Your County Name"
# For example,
# state_name = "Florida"
# county_name = "Broward County"
# Make sure the county ends in the word "County"
# Import libraries
import pandas as pd # For working with data
pd.set_option("display.max_columns", None) # Lets us see all columns of the data instead of just a preview
import geopandas as gpd # For working with spatial data
import numpy as np # For working with numbers
from datetime import date # For working with dates
import matplotlib.pyplot as plt # For making graphs
from matplotlib.colors import to_rgb # For getting colors
import branca.colormap as cm # For creating color scales
import seaborn as sns # For more options to load colors and create plots
import folium # For creating interactive maps
import ee # For loading environmental data with Google Earth Engine
import geemap # For creating interactive maps for Google Earth Engine data
# Initialize Earth Engine to load environmental data
ee.Authenticate()
ee.Initialize(project = projectID)
# Load data
mosquito = gpd.read_file('https://github.com/geo-di-lab/emerge-lessons/raw/refs/heads/main/docs/data/globe_mosquito.zip')
land_cover = gpd.read_file('https://github.com/geo-di-lab/emerge-lessons/raw/refs/heads/main/docs/data/globe_land_cover.zip')
counties = gpd.read_file('https://github.com/geo-di-lab/emerge-lessons/raw/refs/heads/main/docs/data/us_counties.zip').to_crs('EPSG:4326')
# Get boundary for chosen county
county = counties.loc[(counties['STATE_NAME'] == state_name) &
(counties['NAMELSAD'] == county_name)]
# Convert county boundaries to Earth Engine format
region = geemap.geopandas_to_ee(county)
# Create an empty map zoomed to the county
county_center = county.iloc[0].geometry.centroid
map = geemap.Map(center=[county_center.y, county_center.x],
zoom=10)
map.add_basemap('CartoDB.DarkMatter')
Now, we will map land surface temperature across the county. First, define a date range which would be used to find the average temperature across those dates.
start_date = "2025-06-01"
end_date = "2025-07-01"
We will visualize NASA remote sensing data from Google Earth Engine: MOD11A1.061 Terra Land Surface Temperature and Emissivity Daily Global 1km. The land surface temperature data are scaled, so to return to the original temperature in Kelvin, the value needs to be multipled by 0.02. Then, we can convert this to Fahrenheit.
def to_fahrenheit(lst):
celsius = lst * 0.02 - 273.15
fahrenheit = celsius * 1.8 + 32
return fahrenheit
def to_lst(fahrenheit):
celsius = (fahrenheit - 32) / 1.8
lst = (celsius + 273.15) / 0.02
return lst
# Load Google Earth Engine data
lst = (
ee.ImageCollection('MODIS/061/MOD11A1')
.filterDate(start_date, end_date)
.select('LST_Day_1km')
.median()
.clip(region)
)
# Define color range
colors = [
'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'
]
# Define minimum and maximum values and color scale
lst_vis = {
'min': to_lst(50),
'max': to_lst(100),
'palette': colors,
}
# Add data to map
map.addLayer(lst, lst_vis, "LST")
display(map)
Below is a colorbar that shows what the colors on the map represent: land surface temperature in degrees Fahrenheit.
plt.figure(figsize=(len(colors), 1))
plt.imshow([ [to_rgb(f'#{c}') for c in colors] ])
plt.text(-0.6, 0.1, '50 °F', va='center', ha='right', fontsize=24)
plt.text(len(colors) - 0.4, 0.1, '100 °F', va='center', ha='left', fontsize=24)
plt.axis('off')
plt.show()
Next, find the average temperature across the county.
mean_lst = lst.reduceRegion(reducer=ee.Reducer.mean(), geometry=region, scale=1000,
maxPixels=1e12).get('LST_Day_1km').getInfo()
print(f"The average temperature for {county_name} from {start_date} to {end_date} is {round(to_fahrenheit(mean_lst), 2)}°F")
Now, we’ll visualize the land use and land cover in the county: USFS Landscape Change Monitoring System v2024.10 (CONUS and OCONUS).
## Land Use
# Reset the map
map = geemap.Map(center=[county_center.y, county_center.x],
zoom=10)
map.add_basemap('CartoDB.DarkMatter')
# Define land cover color palette
palette = ['fbff97', 'e6558b', '004e2b', '9dbac5', 'a6976a', '1b1716']
visual = {'min': 1, 'max': 6, 'palette': palette}
# Add 1985 land use
landcover_1985 = (
ee.ImageCollection('USFS/GTAC/LCMS/v2024-10')
.filterDate('1985', '1986')
.filter('study_area == "CONUS"')
.first()
.clip(region)
)
map.addLayer(landcover_1985.select('Land_Use'), visual, '1985 Land Use')
# Add 2024 land use
landcover_2024 = (
ee.ImageCollection('USFS/GTAC/LCMS/v2024-10')
.filterDate('2024', '2025')
.filter('study_area == "CONUS"')
.first()
.clip(region)
)
map.addLayer(landcover_2024.select('Land_Use'), visual, '2024 Land Use')
display(map)
We will create a function that calculates the percent of each type of land use and land cover in the county.
def land_stats(image, name, labels):
"""
Calculates the percent of each type of land use or land cover in the county.
Inputs:
- image: Earth Engine image
- name: Land_Use or Land_Cover
- labels: dictionary of labels
"""
count = image.select(name).reduceRegion(ee.Reducer.frequencyHistogram(), geometry=region, scale=30, maxPixels=1e12).getInfo()
land_cover_counts = {}
for key, value in labels.items():
if key in count[name]:
land_cover_counts[value] = count[name][key]
total_pixels = sum([land_cover_counts[i] for i in land_cover_counts])
for label, num in land_cover_counts.items():
percent = round(100 * num / total_pixels, 1)
if percent > 0:
print(f"{label}: {percent}%")
Get the land use for 1985.
land_use_labels = {'1': 'Agriculture',
'2': 'Developed',
'3': 'Forest',
'4': 'Other',
'5': 'Rangeland or Pasture',
'6': 'Non-Processing Area Mask'}
land_stats(landcover_1985, 'Land_Use', land_use_labels)
Now, we repeat for 2024.
land_stats(landcover_2024, 'Land_Use', land_use_labels)
Visualize land cover on a map.
## Land Cover
# Reset the map
map = geemap.Map(center=[county_center.y, county_center.x],
zoom=10)
map.add_basemap('CartoDB.DarkMatter')
# Define color palette
palette = ['004e2b', '009344', '61bb46', 'acbb67', '8b8560', 'cafd4b', 'f89a1c', '8fa55f', 'bebb8e', 'e5e98a', 'ddb925', '893f54', 'e4f5fd', '00b6f0', '1b1716']
visual = {'min': 1, 'max': 15, 'palette': palette}
# Add land cover for 1985 and 2024
map.addLayer(landcover_1985.select('Land_Cover'), visual, '1985 Land Cover')
map.addLayer(landcover_2024.select('Land_Cover'), visual, '2024 Land Cover')
display(map)
Find the percent of each type of land cover
land_cover_labels = {
'1': 'Trees',
'2': 'Tall Shrubs & Trees Mix (AK Only)',
'3': 'Shrubs & Trees Mix',
'4': 'Grass/Forb/Herb & Trees Mix',
'5': 'Barren & Trees Mix',
'6': 'Tall Shrubs (AK Only)',
'7': 'Shrubs',
'8': 'Grass/Forb/Herb & Shrubs Mix',
'9': 'Barren & Shrubs Mix',
'10': 'Grass/Forb/Herb',
'11': 'Barren & Grass/Forb/Herb Mix',
'12': 'Barren or Impervious',
'13': 'Snow or Ice',
'14': 'Water',
'15': 'Non-Processing Area Mask'
}
land_stats(landcover_1985, 'Land_Cover', land_cover_labels)
land_stats(landcover_2024, 'Land_Cover', land_cover_labels)
Tying it all together: Getting Remote Sensing Data for a GLOBE point#
View the list of observations for your chosen county.
# Show mosquito observations for chosen county
data_county = mosquito.sjoin(county, how="inner", predicate="within")
data_county
Choose one of these obsarvations and note the index (the left-most number, like 6643 or 29024) of the observation you chose. Write it below.
index =
# For example,
# index = 29024
Show the data for the point at this index:
point = data_county[data_county.index == index]
point
Get environmental data (land surface temperature, precipitation, and vegetation) for the point you chose. Read more about this code in chapter 4 of the digital textbook.
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)
def get_data_for_point(point):
# Extract the point's date property
point_date = ee.Date(point.get('MeasuredDate'))
day_before = point_date.advance(-1, 'day')
month_before = point_date.advance(-30, 'day')
#### Land Surface Temperature from NASA's MODIS ####
lst = (
ee.ImageCollection('MODIS/061/MOD11A1')
.filterDate(month_before, point_date)
.select(['LST_Day_1km', 'QC_Day', 'LST_Night_1km', 'QC_Night', 'Clear_day_cov', 'Clear_night_cov'])
.median()
# Sample LST at the point’s location
.reduceRegion(
reducer=ee.Reducer.median(),
geometry=point.geometry(),
scale=1000
)
)
point = point.set(lst)
point = point.set('LST_Day_1km_F', to_fahrenheit(point.get('LST_Day_1km').getInfo()))
#### Sentinel-2, Median of Previous Month (used to calculate NDVI) ####
sentinel2 = (
ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterDate(month_before, point_date)
.map(mask_s2_clouds)
.median()
# Sample median band values at the point’s location
.reduceRegion(
reducer=ee.Reducer.median(),
geometry=point.geometry(),
scale=10
)
)
point = point.set(sentinel2)
#### Precipitation from NASA's Daymet V4 ####
rain = (
ee.ImageCollection('NASA/ORNL/DAYMET_V4')
.filterDate(day_before, point_date)
.select('prcp') # Daily total precipitation
.sum()
# Sample median band values at the point’s location
.reduceRegion(
reducer=ee.Reducer.median(),
geometry=point.geometry(),
scale=1000
)
)
point = point.set(rain)
# Add additional information
b8 = ee.Number(point.get('B8'))
b4 = ee.Number(point.get('B4'))
point = point.set('NDVI', b8.subtract(b4).divide(b8.add(b4)))
return point
Use this function to get precipitation, temperature, and vegetation data for your chosen point.
point_fc = geemap.geopandas_to_ee(point).first()
result = get_data_for_point(point_fc)
Print the data:
The Normalized Difference Vegetation Index (NDVI) is a value from -1 to 1 representing vegetation levels. High values close to 1 indicate areas of high vegetation like forests, whereas low values close to -1 indicate areas of low vegetation.
print(f"GLOBE Mosquito observation in {county_name}, {state_name}")
print(f"Date: {result.get('MeasuredDate').getInfo()}")
print(f"Water source: {result.get('WaterSource').getInfo()}")
print(f"Precipitation: {round(result.get('prcp').getInfo(), 2)} mm")
print(f"Land Surface Temperature: {round(result.get('LST_Day_1km_F').getInfo(), 2)} \N{degree sign}F")
print(f"Normalized Difference Vegetation Index (NDVI): {round(result.get('NDVI').getInfo(), 2)}")
Below are more resources for GLOBE and EMERGE: