Steps to create Landsat timeseries:
Submit
button to create timeseries of Landsat imagery and normalized difference indices.Web Apps: https://gishub.org/water-app, https://gishub.org/water-ngrok
Contact: Dr. Qiusheng Wu (Website, LinkedIn, Twitter, YouTube)
# Check geemap installation
import subprocess
try:
import geemap
except ImportError:
print("geemap package is not installed. Installing ...")
subprocess.check_call(["python", "-m", "pip", "install", "geemap"])
# Import libraries
import os
import ee
import geemap
import ipywidgets as widgets
from bqplot import pyplot as plt
from ipyleaflet import WidgetControl
# Create an interactive map
Map = geemap.Map(center=[37.71, 105.47], zoom=4, add_google_map=False)
Map.add_basemap("HYBRID")
Map.add_basemap("ROADMAP")
# Add Earth Engine data
fc = ee.FeatureCollection("users/giswqs/public/chn_admin_level2")
Map.addLayer(fc, {}, "二级行政区")
Map
# Designe interactive widgets
style = {"description_width": "initial"}
output_widget = widgets.Output(layout={"border": "1px solid black"})
output_control = WidgetControl(widget=output_widget, position="bottomright")
Map.add_control(output_control)
admin1_widget = widgets.Text(
description="一级行政区:", value="广东省", width=200, style=style
)
admin2_widget = widgets.Text(
description="二级行政区:", value="广州市", width=300, style=style
)
band_combo = widgets.Dropdown(
description="显示影像波段组合:",
options=[
"Red/Green/Blue",
"NIR/Red/Green",
"SWIR2/SWIR1/NIR",
"NIR/SWIR1/Red",
"SWIR2/NIR/Red",
"SWIR2/SWIR1/Red",
"SWIR1/NIR/Blue",
"NIR/SWIR1/Blue",
"SWIR2/NIR/Green",
"SWIR1/NIR/Red",
],
value="NIR/Red/Green",
style=style,
)
year_widget = widgets.IntSlider(
min=1984,
max=2020,
value=2010,
description="显示指定年份数据:",
width=400,
style=style,
)
fmask_widget = widgets.Checkbox(
value=True, description="应用fmask(去除云, 阴影, 雪)", style=style
)
# Normalized Satellite Indices: https://www.usna.edu/Users/oceano/pguth/md_help/html/norm_sat.htm
nd_options = [
"Vegetation Index (NDVI)",
"Water Index (NDWI)",
"Modified Water Index (MNDWI)",
"Snow Index (NDSI)",
"Soil Index (NDSI)",
"Burn Ratio (NBR)",
"Customized",
]
nd_indices = widgets.Dropdown(
options=nd_options,
value="Modified Water Index (MNDWI)",
description="归一化指数:",
style=style,
)
first_band = widgets.Dropdown(
description="波段1:",
options=["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"],
value="Green",
style=style,
)
second_band = widgets.Dropdown(
description="波段2:",
options=["Blue", "Green", "Red", "NIR", "SWIR1", "SWIR2"],
value="SWIR1",
style=style,
)
nd_threshold = widgets.FloatSlider(
value=0,
min=-1,
max=1,
step=0.01,
description="阈值:",
orientation="horizontal",
style=style,
)
nd_color = widgets.ColorPicker(
concise=False, description="颜色:", value="blue", style=style
)
def nd_index_change(change):
if nd_indices.value == "Vegetation Index (NDVI)":
first_band.value = "NIR"
second_band.value = "Red"
elif nd_indices.value == "Water Index (NDWI)":
first_band.value = "NIR"
second_band.value = "SWIR1"
elif nd_indices.value == "Modified Water Index (MNDWI)":
first_band.value = "Green"
second_band.value = "SWIR1"
elif nd_indices.value == "Snow Index (NDSI)":
first_band.value = "Green"
second_band.value = "SWIR1"
elif nd_indices.value == "Soil Index (NDSI)":
first_band.value = "SWIR1"
second_band.value = "NIR"
elif nd_indices.value == "Burn Ratio (NBR)":
first_band.value = "NIR"
second_band.value = "SWIR2"
elif nd_indices.value == "Customized":
first_band.value = None
second_band.value = None
nd_indices.observe(nd_index_change, names="value")
submit = widgets.Button(
description="Submit", button_style="primary", tooltip="Click me", style=style
)
full_widget = widgets.VBox(
[
widgets.HBox([admin1_widget, admin2_widget]),
widgets.HBox([band_combo, year_widget, fmask_widget]),
widgets.HBox([nd_indices, first_band, second_band, nd_threshold, nd_color]),
submit,
]
)
full_widget
# Capture user interaction with the map
def handle_interaction(**kwargs):
latlon = kwargs.get("coordinates")
if kwargs.get("type") == "click":
Map.default_style = {"cursor": "wait"}
xy = ee.Geometry.Point(latlon[::-1])
selected_fc = fc.filterBounds(xy)
with output_widget:
output_widget.clear_output()
try:
admin1_id = selected_fc.first().get("ADM1_ZH").getInfo()
admin2_id = selected_fc.first().get("ADM2_ZH").getInfo()
admin1_widget.value = admin1_id
admin2_widget.value = admin2_id
Map.layers = Map.layers[:4]
geom = selected_fc.geometry()
layer_name = admin1_id + "-" + admin2_id
Map.addLayer(
ee.Image().paint(geom, 0, 2), {"palette": "red"}, layer_name
)
print(layer_name)
except:
print("找不到相关行政区")
Map.layers = Map.layers[:4]
Map.default_style = {"cursor": "pointer"}
Map.on_interaction(handle_interaction)
# Click event handler
def submit_clicked(b):
with output_widget:
output_widget.clear_output()
print("Computing...")
Map.default_style = {"cursor": "wait"}
try:
admin1_id = admin1_widget.value
admin2_id = admin2_widget.value
band1 = first_band.value
band2 = second_band.value
selected_year = year_widget.value
threshold = nd_threshold.value
bands = band_combo.value.split("/")
apply_fmask = fmask_widget.value
palette = nd_color.value
roi = fc.filter(
ee.Filter.And(
ee.Filter.eq("ADM1_ZH", admin1_id),
ee.Filter.eq("ADM2_ZH", admin2_id),
)
)
Map.layers = Map.layers[:4]
geom = roi.geometry()
layer_name = admin1_id + "-" + admin2_id
Map.addLayer(ee.Image().paint(geom, 0, 2), {"palette": "red"}, layer_name)
images = geemap.landsat_timeseries(
roi=roi,
start_year=1984,
end_year=2020,
start_date="01-01",
end_date="12-31",
apply_fmask=apply_fmask,
)
nd_images = images.map(lambda img: img.normalizedDifference([band1, band2]))
result_images = nd_images.map(lambda img: img.gt(threshold))
selected_image = ee.Image(
images.toList(images.size()).get(selected_year - 1984)
)
selected_result_image = ee.Image(
result_images.toList(result_images.size()).get(selected_year - 1984)
).selfMask()
vis_params = {"bands": bands, "min": 0, "max": 3000}
Map.addLayer(selected_image, vis_params, "Landsat " + str(selected_year))
Map.addLayer(
selected_result_image,
{"palette": palette},
"Result " + str(selected_year),
)
def cal_area(img):
pixel_area = img.multiply(ee.Image.pixelArea()).divide(1e6)
img_area = pixel_area.reduceRegion(
**{
"geometry": roi.geometry(),
"reducer": ee.Reducer.sum(),
"scale": 1000,
"maxPixels": 1e12,
"bestEffort": True,
}
)
return img.set({"area": img_area})
areas = result_images.map(cal_area)
stats = areas.aggregate_array("area").getInfo()
x = list(range(1984, 2021))
y = [item.get("nd") for item in stats]
fig = plt.figure(1)
fig.layout.height = "270px"
plt.clear()
plt.plot(x, y)
plt.title("Temporal trend (1984-2020)")
plt.xlabel("Year")
plt.ylabel("Area (km2)")
output_widget.clear_output()
plt.show()
except Exception as e:
print(e)
print("An error occurred during computation.")
Map.default_style = {"cursor": "default"}
submit.on_click(submit_clicked)