11. Earth Engine Applications#
11.1. Introduction#
11.2. Technical requirements#
conda install -n base mamba -c conda-forge
mamba create -n gee -c conda-forge geemap pygis
conda activate gee
jupyter lab
# %pip install pygis
import ee
import geemap
geemap.ee_initialize()
11.3. Analyzing surface water dynamics#
11.3.1. Surface water occurrence#
dataset = ee.Image('JRC/GSW1_4/GlobalSurfaceWater')
dataset.bandNames()
Map = geemap.Map()
Map.add_basemap('HYBRID')
Map
image = dataset.select(['occurrence'])
region = Map.user_roi # Draw a polygon on the map
if region is None:
region = ee.Geometry.BBox(-99.957, 46.8947, -99.278, 47.1531)
vis_params = {'min': 0.0, 'max': 100.0, 'palette': ['ffffff', 'ffbbbb', '0000ff']}
Map.addLayer(image, vis_params, 'Occurrence')
Map.addLayer(region, {}, 'ROI', True, 0.5)
Map.centerObject(region)
Map.add_colorbar(vis_params, label='Water occurrence (%)', layer_name='Occurrence')
df = geemap.image_histogram(
image,
region,
scale=30,
return_df=True,
)
df
hist = geemap.image_histogram(
image,
region,
scale=30,
x_label='Water Occurrence (%)',
y_label='Pixel Count',
title='Surface Water Occurrence',
layout_args={'title': dict(x=0.5)},
return_df=False,
)
hist
# !pip install -U kaleido
hist.update_layout(
autosize=False,
width=800,
height=400,
margin=dict(l=30, r=20, b=10, t=50, pad=4)
)
hist.write_image('water_occurrence.jpg', scale=2)
11.3.2. Surface water monthly history#
dataset = ee.ImageCollection('JRC/GSW1_4/MonthlyHistory')
dataset.size()
dataset.aggregate_array("system:index")
Map = geemap.Map()
Map
image = dataset.filterDate('2020-08-01', '2020-09-01').first()
region = Map.user_roi # Draw a polygon on the map
if region is None:
region = ee.Geometry.BBox(-99.957, 46.8947, -99.278, 47.1531)
vis_params = {'min': 0.0, 'max': 2.0, 'palette': ['ffffff', 'fffcb8', '0905ff']}
Map.addLayer(image, vis_params, 'Water')
Map.addLayer(region, {}, 'ROI', True, 0.5)
Map.centerObject(region)
geemap.jrc_hist_monthly_history(
region=region, scale=30, frequency='month', denominator=1e4, y_label='Area (ha)'
)
geemap.jrc_hist_monthly_history(
region=region,
start_month=6,
end_month=9,
scale=30,
frequency='month',
denominator=1e4,
y_label='Area (ha)',
)
geemap.jrc_hist_monthly_history(
region=region,
start_month=6,
end_month=9,
scale=30,
frequency='year',
reducer='mean',
denominator=1e4,
y_label='Area (ha)',
)
geemap.jrc_hist_monthly_history(
region=region,
start_month=6,
end_month=9,
scale=30,
frequency='year',
reducer='max',
denominator=1e4,
y_label='Area (ha)',
)
11.4. Mapping flood extent#
11.4.1. Create an interactive map#
Map = geemap.Map(center=[29.3055, 68.9062], zoom=6)
Map
11.4.2. Search datasets#
country_name = 'Pakistan'
pre_flood_start_date = '2021-08-01'
pre_flood_end_date = '2021-09-30'
flood_start_date = '2022-08-01'
flood_end_date = '2022-09-30'
11.4.3. Visualize datasets#
country = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017').filter(
ee.Filter.eq('country_na', country_name)
)
style = {'color': 'black', 'fillColor': '00000000'}
Map.addLayer(country.style(**style), {}, country_name)
Map.centerObject(country)
Map
11.4.4. Create Landsat composites#
landsat_col_2021 = (
ee.ImageCollection('LANDSAT/LC08/C02/T1')
.filterDate(pre_flood_start_date, pre_flood_end_date)
.filterBounds(country)
)
landsat_2021 = ee.Algorithms.Landsat.simpleComposite(landsat_col_2021).clipToCollection(
country
)
vis_params = {'bands': ['B6', 'B5', 'B4'], 'max': 128}
Map.addLayer(landsat_2021, vis_params, 'Landsat 2021')
landsat_col_2022 = (
ee.ImageCollection('LANDSAT/LC08/C02/T1')
.filterDate(flood_start_date, flood_end_date)
.filterBounds(country)
)
landsat_2022 = ee.Algorithms.Landsat.simpleComposite(landsat_col_2022).clipToCollection(
country
)
Map.addLayer(landsat_2022, vis_params, 'Landsat 2022')
Map.centerObject(country)
Map
11.4.5. Compare Landsat composites side by side#
Map = geemap.Map()
Map.setCenter(68.4338, 26.4213, 7)
left_layer = geemap.ee_tile_layer(
landsat_2021, vis_params, 'Landsat 2021'
)
right_layer = geemap.ee_tile_layer(
landsat_2022, vis_params, 'Landsat 2022'
)
Map.split_map(
left_layer, right_layer, left_label='Landsat 2021', right_label='Landsat 2022'
)
Map.addLayer(country.style(**style), {}, country_name)
Map
11.4.6. Compute Normalized Difference Water Index (NDWI)#
ndwi_2021 = landsat_2021.normalizedDifference(['B3', 'B5']).rename('NDWI')
ndwi_2022 = landsat_2022.normalizedDifference(['B3', 'B5']).rename('NDWI')
Map = geemap.Map()
Map.setCenter(68.4338, 26.4213, 7)
ndwi_vis = {'min': -1, 'max': 1, 'palette': 'ndwi'}
left_layer = geemap.ee_tile_layer(ndwi_2021, ndwi_vis, 'NDWI 2021')
right_layer = geemap.ee_tile_layer(ndwi_2022, ndwi_vis, 'NDWI 2022')
Map.split_map(left_layer, right_layer, left_label='NDWI 2021', right_label='NDWI 2022')
Map.addLayer(country.style(**style), {}, country_name)
Map
11.4.7. Extract Landsat water extent#
threshold = 0.1
water_2021 = ndwi_2021.gt(threshold).selfMask()
water_2022 = ndwi_2022.gt(threshold).selfMask()
Map = geemap.Map()
Map.setCenter(68.4338, 26.4213, 7)
Map.addLayer(landsat_2021, vis_params, 'Landsat 2021', False)
Map.addLayer(landsat_2022, vis_params, 'Landsat 2022', False)
left_layer = geemap.ee_tile_layer(
water_2021, {'palette': 'blue'}, 'Water 2021'
)
right_layer = geemap.ee_tile_layer(
water_2022, {'palette': 'red'}, 'Water 2022'
)
Map.split_map(
left_layer, right_layer, left_label='Water 2021', right_label='Water 2022'
)
Map.addLayer(country.style(**style), {}, country_name)
Map
11.4.8. Extract Landsat flood extent#
flood_extent = water_2022.unmask().subtract(water_2021.unmask()).gt(0).selfMask()
Map = geemap.Map()
Map.setCenter(68.4338, 26.4213, 7)
Map.addLayer(landsat_2021, vis_params, 'Landsat 2021', False)
Map.addLayer(landsat_2022, vis_params, 'Landsat 2022', False)
left_layer = geemap.ee_tile_layer(
water_2021, {'palette': 'blue'}, 'Water 2021'
)
right_layer = geemap.ee_tile_layer(
water_2022, {'palette': 'red'}, 'Water 2022'
)
Map.split_map(
left_layer, right_layer, left_label='Water 2021', right_label='Water 2022'
)
Map.addLayer(flood_extent, {'palette': 'cyan'}, 'Flood Extent')
Map.addLayer(country.style(**style), {}, country_name)
Map
11.4.9. Calculate Landsat flood area#
area_2021 = geemap.zonal_stats(
water_2021, country, scale=1000, statistics_type='SUM', return_fc=True
)
geemap.ee_to_df(area_2021)
area_2022 = geemap.zonal_stats(
water_2022, country, scale=1000, statistics_type='SUM', return_fc=True
)
geemap.ee_to_df(area_2022)
flood_area = geemap.zonal_stats(
flood_extent, country, scale=1000, statistics_type='SUM', return_fc=True
)
geemap.ee_to_df(flood_area)
11.4.10. Create Sentinel-1 SAR composites#
s1_col_2021 = (
ee.ImageCollection('COPERNICUS/S1_GRD')
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
.filterDate(pre_flood_start_date, pre_flood_end_date)
.filterBounds(country)
.select('VV')
)
s1_col_2021
s1_col_2022 = (
ee.ImageCollection('COPERNICUS/S1_GRD')
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
.filterDate(flood_start_date, flood_end_date)
.filterBounds(country)
.select('VV')
)
s1_col_2022
Map = geemap.Map()
Map.add_basemap('HYBRID')
sar_2021 = s1_col_2021.reduce(ee.Reducer.percentile([20])).clipToCollection(country)
sar_2022 = s1_col_2022.reduce(ee.Reducer.percentile([20])).clipToCollection(country)
Map.addLayer(sar_2021, {'min': -25, 'max': -5}, 'SAR 2021')
Map.addLayer(sar_2022, {'min': -25, 'max': -5}, 'SAR 2022')
Map.centerObject(country)
Map
11.4.11. Apply speckle filtering#
col_2021 = s1_col_2021.map(lambda img: img.focal_median(100, 'circle', 'meters'))
col_2022 = s1_col_2022.map(lambda img: img.focal_median(100, 'circle', 'meters'))
Map = geemap.Map()
Map.add_basemap('HYBRID')
sar_2021 = col_2021.reduce(ee.Reducer.percentile([20])).clipToCollection(country)
sar_2022 = col_2022.reduce(ee.Reducer.percentile([20])).clipToCollection(country)
Map.addLayer(sar_2021, {'min': -25, 'max': -5}, 'SAR 2021')
Map.addLayer(sar_2022, {'min': -25, 'max': -5}, 'SAR 2022')
Map.centerObject(country)
Map
11.4.12. Compare Sentinel-1 SAR composites side by side#
Map = geemap.Map()
Map.setCenter(68.4338, 26.4213, 7)
left_layer = geemap.ee_tile_layer(sar_2021, {'min': -25, 'max': -5}, 'SAR 2021')
right_layer = geemap.ee_tile_layer(sar_2022, {'min': -25, 'max': -5}, 'SAR 2022')
Map.split_map(
left_layer, right_layer, left_label='Sentinel-1 2021', right_label='Sentinel-1 2022'
)
Map.addLayer(country.style(**style), {}, country_name)
Map
11.4.13. Extract SAR water extent#
threshold = -18
water_2021 = sar_2021.lt(threshold)
water_2022 = sar_2022.lt(threshold)
Map = geemap.Map()
Map.setCenter(68.4338, 26.4213, 7)
Map.addLayer(sar_2021, {'min': -25, 'max': -5}, 'SAR 2021')
Map.addLayer(sar_2022, {'min': -25, 'max': -5}, 'SAR 2022')
left_layer = geemap.ee_tile_layer(
water_2021.selfMask(), {'palette': 'blue'}, 'Water 2021'
)
right_layer = geemap.ee_tile_layer(
water_2022.selfMask(), {'palette': 'red'}, 'Water 2022'
)
Map.split_map(
left_layer, right_layer, left_label='Water 2021', right_label='Water 2022'
)
Map.addLayer(country.style(**style), {}, country_name)
Map
11.4.14. Extract SAR flood extent#
flood_extent = water_2022.unmask().subtract(water_2021.unmask()).gt(0).selfMask()
Map = geemap.Map()
Map.setCenter(68.4338, 26.4213, 7)
Map.addLayer(sar_2021, {'min': -25, 'max': -5}, 'SAR 2021')
Map.addLayer(sar_2022, {'min': -25, 'max': -5}, 'SAR 2022')
left_layer = geemap.ee_tile_layer(
water_2021, {'palette': 'blue'}, 'Water 2021'
)
right_layer = geemap.ee_tile_layer(
water_2022, {'palette': 'red'}, 'Water 2022'
)
Map.split_map(
left_layer, right_layer, left_label='Water 2021', right_label='Water 2022'
)
Map.addLayer(flood_extent, {'palette': 'cyan'}, 'Flood Extent')
Map.addLayer(country.style(**style), {}, country_name)
Map
11.4.15. Calculate SAR flood area#
area_2021 = geemap.zonal_stats(
water_2021, country, scale=1000, statistics_type='SUM', return_fc=True
)
geemap.ee_to_df(area_2021)
area_2022 = geemap.zonal_stats(
water_2022, country, scale=1000, statistics_type='SUM', return_fc=True
)
geemap.ee_to_df(area_2022)
flood_area = geemap.zonal_stats(
flood_extent, country, scale=1000, statistics_type='SUM', return_fc=True
)
geemap.ee_to_df(flood_area)
11.5. Forest cover change analysis#
11.5.1. Forest cover mapping#
dataset = ee.Image('UMD/hansen/global_forest_change_2021_v1_9')
dataset.bandNames()
Map = geemap.Map()
first_bands = ['first_b50', 'first_b40', 'first_b30']
first_image = dataset.select(first_bands)
Map.addLayer(first_image, {'bands': first_bands, 'gamma': 1.5}, 'Year 2000 Bands 5/4/3')
last_bands = ['last_b50', 'last_b40', 'last_b30']
last_image = dataset.select(last_bands)
Map.addLayer(last_image, {'bands': last_bands, 'gamma': 1.5}, 'Year 2021 Bands 5/4/3')
treecover = dataset.select(['treecover2000'])
treeCoverVisParam = {'min': 0, 'max': 100, 'palette': ['black', 'green']}
name = 'Tree cover (%)'
Map.addLayer(treecover, treeCoverVisParam, name)
Map.add_colorbar(treeCoverVisParam, label=name, layer_name=name)
Map
threshold = 10
treecover_bin = treecover.gte(threshold).selfMask()
treeVisParam = {'palette': ['green']}
Map.addLayer(treecover_bin, treeVisParam, 'Tree cover bin')
11.5.2. Forest loss and gain mapping#
Map = geemap.Map()
treeloss_year = dataset.select(['lossyear'])
treeLossVisParam = {'min': 0, 'max': 21, 'palette': ['yellow', 'red']}
layer_name = 'Tree loss year'
Map.addLayer(treeloss_year, treeLossVisParam, layer_name)
Map.add_colorbar(treeLossVisParam, label=layer_name, layer_name=layer_name)
Map
treeloss = dataset.select(['loss']).selfMask()
treegain = dataset.select(['gain']).selfMask()
Map.addLayer(treeloss, {'palette': 'red'}, 'Tree loss')
Map.addLayer(treegain, {'palette': 'yellow'}, 'Tree gain')
Map
11.5.3. Zonal statistics by country#
Map = geemap.Map()
countries = ee.FeatureCollection(geemap.examples.get_ee_path('countries'))
style = {'color': '#000000ff', 'fillColor': '#00000000'}
Map.addLayer(countries.style(**style), {}, 'Countries')
Map
geemap.zonal_stats_by_group(
treecover_bin,
countries,
'forest_cover.csv',
statistics_type='SUM',
denominator=1e6,
scale=1000,
)
geemap.pie_chart(
'forest_cover.csv', names='NAME', values='Class_sum', max_rows=20, height=400
)
geemap.bar_chart(
'forest_cover.csv',
x='NAME',
y='Class_sum',
max_rows=20,
x_label='Country',
y_label='Forest area (km2)',
)
geemap.zonal_stats_by_group(
treeloss,
countries,
'treeloss.csv',
statistics_type='SUM',
denominator=1e6,
scale=1000,
)
geemap.pie_chart(
'treeloss.csv', names='NAME', values='Class_sum', max_rows=20, height=600
)
geemap.bar_chart(
'treeloss.csv',
x='NAME',
y='Class_sum',
max_rows=20,
x_label='Country',
y_label='Forest loss area (km2)',
)
11.6. Global land cover mapping#
11.6.1. Dynamic World#
Map = geemap.Map()
region = ee.Geometry.BBox(-89.7088, 42.9006, -89.0647, 43.2167)
start_date = '2021-01-01'
end_date = '2022-01-01'
image = geemap.dynamic_world_s2(region, start_date, end_date, clip=True)
vis_params = {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 3000}
Map.addLayer(image, vis_params, 'Sentinel-2 image')
landcover = geemap.dynamic_world(
region, start_date, end_date, return_type='hillshade', clip=True
)
Map.addLayer(landcover, {}, 'Land Cover')
Map.add_legend(
title="Dynamic World Land Cover",
builtin_legend='Dynamic_World',
layer_name='Land Cover',
)
Map.centerObject(region, 13)
Map
classes = geemap.dynamic_world(
region, start_date, end_date, return_type='class', clip=True
)
vis_params = {
"min": 0,
"max": 8,
"palette": [
"#419BDF",
"#397D49",
"#88B053",
"#7A87C6",
"#E49635",
"#DFC35A",
"#C4281B",
"#A59B8F",
"#B39FE1",
],
}
Map.addLayer(classes, vis_params, 'Class')
probability = geemap.dynamic_world(
region, start_date, end_date, return_type='visualize', clip=True
)
Map.addLayer(probability, {}, 'Visualize')
probability = geemap.dynamic_world(
region, start_date, end_date, return_type='probability', clip=True
)
Map.addLayer(probability, {}, 'Probability')
df = geemap.image_area_by_group(classes, region=region, scale=10, denominator=1e6)
df
Map = geemap.Map()
region = ee.Geometry.BBox(-89.7088, 42.9006, -89.0647, 43.2167)
start_date = '2017-01-01'
end_date = '2022-12-31'
images = geemap.dynamic_world_timeseries(
region, start_date, end_date, frequency='year', return_type="hillshade"
)
Map.ts_inspector(images, date_format='YYYY')
Map.add_legend(title="Dynamic World Land Cover", builtin_legend='Dynamic_World')
Map.centerObject(region)
Map
11.6.2. ESA WorldCover#
Map = geemap.Map()
esa_2020 = ee.ImageCollection("ESA/WorldCover/v100").first()
esa_2021 = ee.ImageCollection("ESA/WorldCover/v200").first()
Map.addLayer(esa_2020, {'bands': ['Map']}, 'ESA 2020')
Map.addLayer(esa_2021, {'bands': ['Map']}, 'ESA 2021')
Map.add_legend(builtin_legend='ESA_WorldCover', title='ESA Land Cover')
Map
df = geemap.image_area_by_group(
esa_2021, scale=1000, denominator=1e6, decimal_places=4, verbose=True
)
df
df.to_csv('esa_2021.csv')
Map = geemap.Map(center=[-88.9088, 43.0006], zoom=12)
# Create Dynamic World Land Cover 2021
region = ee.Geometry.BBox(-179, -89, 179, 89)
start_date = '2021-01-01'
end_date = '2022-01-01'
dw_2021 = geemap.dynamic_world(region, start_date, end_date, return_type='hillshade')
# Create ESA WorldCover 2021
esa_2021 = ee.ImageCollection("ESA/WorldCover/v200").first()
# Create a split map
left_layer = geemap.ee_tile_layer(esa_2021, {'bands': ['Map']}, "ESA 2021")
right_layer = geemap.ee_tile_layer(dw_2021, {}, "Dynamic World 2021")
Map.split_map(left_layer, right_layer)
# Add legends
Map.add_legend(
title="ESA WorldCover", builtin_legend='ESA_WorldCover', position='bottomleft'
)
Map.add_legend(
title="Dynamic World Land Cover",
builtin_legend='Dynamic_World',
position='bottomright',
)
Map
11.6.3. Esri global land cover#
def esri_annual_land_cover(year):
collection = ee.ImageCollection('projects/sat-io/open-datasets/landcover/ESRI_Global-LULC_10m_TS')
start_date = ee.Date.fromYMD(year, 1, 1)
end_date = start_date.advance(1, 'year')
image = collection.filterDate(start_date, end_date).mosaic()
return image.set('system:time_start', start_date.millis())
start_year = 2017
end_year = 2021
years = ee.List.sequence(start_year, end_year)
images = ee.ImageCollection(years.map(esri_annual_land_cover))
images
palette = [
"#1A5BAB",
"#358221",
"#000000",
"#87D19E",
"#FFDB5C",
"#000000",
"#ED022A",
"#EDE9E4",
"#F2FAFF",
"#C8C8C8",
"#C6AD8D",
]
vis_params = {"min": 1, "max": 11, "palette": palette}
Map = geemap.Map()
Map.ts_inspector(images, left_vis=vis_params, date_format='YYYY')
Map.add_legend(title="Esri Land Cover", builtin_legend='ESRI_LandCover')
Map
Map = geemap.Map(center=[43.0006, -88.9088], zoom=12)
# Create Dynamic World Land Cover 2021
dw_2021 = geemap.dynamic_world(region, start_date, end_date, return_type='hillshade')
# Create Esri Global Land Cover 2021
esri_2021 = esri_annual_land_cover(2021)
# Create a split map
left_layer = geemap.ee_tile_layer(esri_2021, vis_params, "Esri 2021")
right_layer = geemap.ee_tile_layer(dw_2021, {}, "Dynamic World 2021")
Map.split_map(left_layer, right_layer)
# Add legends
Map.add_legend(
title="Esri Land Cover", builtin_legend='ESRI_LandCover', position='bottomleft'
)
Map.add_legend(
title="Dynamic World Land Cover",
builtin_legend='Dynamic_World',
position='bottomright',
)
Map