练习1 基于GEE的哨兵1号数据获取及水体提取取

1. 数据获取

本案例需要科学上网，成功安装geemap后方可使用，slope_correction函数封装在utils中，有需要的话请联系作者yujunchuan@mail.cgs.gov.cn

import geemap #需安装geemap
import ee
# ee.Authenticate()
geemap.set_proxy(port=33210) #根据自己电脑配置修改
ee.Initialize()
import math
from utils import *

定义必要的函数

def clip_img(img):
    return img.clip(region)

def filtter_speckles(img):
    vv=img.select('VV')
    vh= img.select('VH')
    vv_smoothed = vv.focal_mean(50,'circle','meters').rename('VV_Filtered')
    vh_smoothed = vh.focal_mean(50,'circle','meters').rename('VH_Filtered')
    return img.addBands(vv_smoothed).addBands(vh_smoothed)

def SDWI(img):
    VH = img.select("VH_mean")
    VV = img.select("VV_mean")
    x = VV.multiply(VH).multiply(ee.Number(10))
    y=(x).log().rename('sdwi')
    # sdwi=y.subtract(8)
    file=img.addBands(y)
    return file.select("sdwi")
def addDate(image):
    date = ee.Date(image.get('system:time_start')).format('yyyyMMdd')
    return image.set("image_date",date)

设定鄱阳湖区域范围

region=ee.Geometry.Polygon([[[115.795898, 28.623104],[116.779175, 28.623104],[116.779175, 29.441989],[115.795898, 29.441989],[115.795898, 28.623104]]])
Map=geemap.Map()
Map.center_object(region)
Map.addLayer(region,{}, 'region')
Map

2. SDWI指数计算

编写sdwi计算函数，如非山区可以不进行地形校正

def sdwi_cal(data):
    corrected = slope_correction(data).main().map(clip_img)
    VV = corrected.mean().select('VV')
    VH = corrected.mean().select('VH')
    # Map.addLayer(VV, {'bands':['VV'], 'min':[-21], 'max':[0.5], 'gamma': 0.65}, 'SAR_RGB_VV')
    # Map.addLayer(VH, {'bands':['VH'], 'min':[-28], 'max':[4], 'gamma': 0.65}, 'SAR_RGB_VH')
    S1 = corrected.map(filtter_speckles)
    name=S1.select("VV_Filtered").first().get("image_date")
    SVV = S1.select("VV_Filtered").reduce(ee.Reducer.mean()).rename("VV_mean").set('image_date',name)
    SVH = S1.select("VH_Filtered").reduce(ee.Reducer.mean()).rename("VH_mean").set('image_date',name)
    Sen1=SVV.addBands(SVH)
    sdwi = SDWI(Sen1).select('sdwi')
    # Map.addLayer(sdwi,{'bands':['sdwi'], 'min':[8], 'max':[10], 'gamma': 0.65}, 'sdwi')
    return sdwi

根据阅读计算sdwi

def monthly_sdwi(mon):
    start=ee.Date.fromYMD(year=2022,month=mon,day=1)
    end=start.advance(1.0,'month')
    data = (ee.ImageCollection('COPERNICUS/S1_GRD')
            .filterBounds(region)
            .filterMetadata('transmitterReceiverPolarisation','equals',["VV", "VH"])
            .filterMetadata('instrumentMode','equals','IW')
            .filterDate(start, end)
            .map(addDate))
    sdwi=sdwi_cal(data)
    name=sdwi.get("image_date")
    return sdwi

monthly=ee.ImageCollection(ee.List.sequence(6,9).map(monthly_sdwi))
bandname=monthly.aggregate_array('image_date')
singlesdwi=monthly.toBands().rename(bandname) 
for i in range(4):
    Map.addLayer(singlesdwi.select(i).updateMask(singlesdwi.select(i).gte(8.5)),  {"palette":['6a5acd']}, "swdi"+str(i))
Map

3. 水面变化分析

通过多视角查看前后水体变化

Map=geemap.Map()
Map.center_object(region)
Map.split_map(left_layer='SATELLITE', right_layer='SATELLITE')
vis= {"palette":['6a5acd']}
left_layer=geemap.ee_tile_layer(singlesdwi.select(0).updateMask(singlesdwi.select(0).gte(8.5)), vis,name='2022-06')
right_layer=geemap.ee_tile_layer(singlesdwi.select(3).updateMask(singlesdwi.select(3).gte(8.5)), vis,name='2022-09')
Map.split_map(left_layer, right_layer)
Map

年度12个月变化监测，提取永久水体

monthly=ee.ImageCollection(ee.List.sequence(1,12).map(monthly_sdwi))
bandname=monthly.aggregate_array('image_date')
singlesdwi=monthly.toBands().rename(bandname) 
for i in range(0,12):
    if i ==0:
        initial=singlesdwi.select(0).where(singlesdwi.select(i).lte(8.5),0).where(singlesdwi.select(i).gt(8.5),1)
    else:
        final=singlesdwi.select(i).where(singlesdwi.select(i).lte(8.5),0).where(singlesdwi.select(i).gt(8.5),1)
        initial=initial.add(final)

Map.addLayer(initial, {"min":[1],"max":[12], "palette":['ffffff', '99d9ea', '0000ff']}, "water")
Map

4. 下载提取结果

for i in range(8):
    img=singlesdwi.select(i)
    img=img.setDefaultProjection('epsg:4326',None,10)
    name=singlesdwi.bandNames().getInfo()
    outname = './sentinel_'+name[i]+str(i)+".tif"
    geemap.download_ee_image(img, filename=outname,scale=10,region=region,crs='EPSG:4326')

想了解更多请关注[45度科研人]公众号，欢迎给我留言！