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)