沐风-GEE系列专栏

涵盖生态环境，卫星遥感，气候变化，云计算。

• Google Earth Engine 系列学习笔记(一)
• Google Earth Engine 系列学习笔记(二)
• GEE-Scholars 可视化时间序列农业景观
• GEE-Scholars 机器学习之卫星图像监督分类

---

本章目录

• 沐风-GEE系列专栏
• 前言
• 一、Sentinel-2 光学数据动态监测水体变化
• • 1、定义监测区域
  • 2、定义监测时间
  • 3、哨兵数据去云
  • 4、NDWI指数提取水体
• 二、Sentinel-1 雷达数据动态监测统计水体变化
• • 1、定义监测区域
  • 2、定义数据监测时间
  • 3、计算水体的算法
  • 4、水体面积以及图表制作
• 总结

---

前言

今年7月份河南省出现历史罕见的持续性强降雨天气，对当地造成了重大的影响。据统计，农作物受灾面积75千公顷，成灾面积25.2千公顷，绝收面积4.7千公顷；直接经济损失54228.72万元。

面对自然灾害，核心城市地区更易受到媒体关注并报道，市政条件和治理条件等多种原因，可以相对更为有效的控制灾情以及开展灾后治理工作；而对于郊区乃至更大范围的农村村镇区域，地广人稀，地势相对复杂。人们利用自媒体发布的视频和照片，对于受灾位置等微观信息可以有直接的作用，但对于受灾的宏观范围无法做到有效反映。

利用无人机航拍以及卫星影像来做受灾范围的评估，可以为抢险救灾提供有效的决策依据，尤其是对更大范围的山区，可以做到全区域覆盖监测。然而传统手段利用客户端ENVI、Arcgis等软件进行图像解译过程较为复杂，尤其是数据下载和预处理几乎可以占据70%以上的处理时间，处理效率亟待提升。所以，一种整合数据和云计算的一体化的平台，尤为重要。

本篇文章基于GEE云平台的光学数据以及雷达数据快速获取河南特大水灾的水淹区域，为抢险救灾提供一手客观依据。

一、Sentinel-2 光学数据动态监测水体变化

1、定义监测区域

// var roi 
Map.centerObject(roi, 7);

2、定义监测时间

var start = ee.Date('2021-06-06')
var med = ee.Date('2021-07-20')
var end = ee.Date('2021-07-30')

3、哨兵数据去云

function mask_cloud_s2(start,end){
  var s2c = ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY');
  var s2Sr = ee.ImageCollection("COPERNICUS/S2")
  var s2c = s2c.filterDate(start, end).filterBounds(roi)
  var s2Sr = s2Sr.filterDate(start, end).filterBounds(roi).filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 40));
  var viz = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};

  function indexJoin(collectionA, collectionB, propertyName) {
    var joined = ee.ImageCollection(ee.Join.saveFirst(propertyName).apply({
      primary: collectionA,
      secondary: collectionB,
      condition: ee.Filter.equals({
        leftField: 'system:index',
        rightField: 'system:index'})
    }));
    return joined.map(function(image) {
      return image.addBands(ee.Image(image.get(propertyName)));
    });
  }
  
  function maskS2clouds01(image) {
    var s2c = image.select('probability');
    var isCloud = s2c.gte(50)
    return image.updateMask(isCloud.not())
      .select("B.*")
      .copyProperties(image, ["system:time_start"]);
  }
  
  var withCloudProbability = indexJoin(s2Sr, s2c, 'cloud_probability');
  var masked_01 = ee.ImageCollection(withCloudProbability.map(maskS2clouds01));
  var viz = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};
  Map.addLayer(withCloudProbability, viz, 'RGB_raw',false);
  var median01 = masked_01.median()
   return median01
}
var median01_before=mask_cloud_s2(start,med)
var viz = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};
Map.addLayer(median01_before, viz, 'RGB_before',false);

var median01_after=mask_cloud_s2(med,end)
var viz = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};
Map.addLayer(median01_after, viz, 'RGB_after',false);

4、NDWI指数提取水体

function NDWI(image){
  var ndwi=image.normalizedDifference(["B8","B3"])//.clipToCollection(vector)
  return ndwi
}
var ndwi_bf=median01_before.normalizedDifference(["B8","B3"])
Map.addLayer(ndwi_bf, {min:0,max:1,palette:['3005da','ffec02']}, 'waterIndex_before',false);
var ndwi_af_raw=median01_after.normalizedDifference(["B8","B3"])
Map.addLayer(ndwi_af_raw, {min:0,max:1,palette:['3005da','ffec02']}, 'waterIndex_after_raw',false);
var ndwi_after=ee.ImageCollection([ndwi_bf,ndwi_af_raw]).mosaic()
Map.addLayer(ndwi_after, {min:0,max:1,palette:['3005da','ffec02']}, 'waterIndex_after',false);
var water_af=ndwi_after.lte(0.15)
var water_bf=ndwi_bf.lte(0.1)
Map.addLayer(water_bf, {min:0,max:1,palette:['fdffed','0647ff']}, 'water_bf',false);
Map.addLayer(water_af, {min:0,max:1,palette:['fdffed','0647ff']}, 'water_af',true);
var image=water_bf.addBands(water_bf).addBands(water_af).rename(['b1','b2','b3'])

var viz = {bands: ['b1','b2','b3'], min: 0, max:1};
Map.addLayer(image, viz, 'RGB_water');

二、Sentinel-1 雷达数据动态监测统计水体变化

1、定义监测区域

// var roi Polygon
Map.centerObject(roi, 7);

2、定义数据监测时间

var s1Col = s1.filterBounds(roi)
              .filterDate("2021-7-01", "2021-8-10")
              .filter(ee.Filter.listContains("transmitterReceiverPolarisation", "VV"))
              .select("VV")
              .map(function(image) {
                var maskedImage = image.mask().and(image.gte(-30));
                return image.updateMask(maskedImage);
              })
              .map(function(image) {
                var time_start = image.get("system:time_start");
                return image.set("date", ee.Date(time_start).format("yyyyMMdd"));
              })
              .sort("system:time_start");

数据概况：处理后的影像

3、计算水体的算法

function OTSU(histogram) {
  var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
  var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
  var size = means.length().get([0]);
  var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
  var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);
  var mean = sum.divide(total);
  var indices = ee.List.sequence(1, size);
  var bss = indices.map(function(i) {
    var aCounts = counts.slice(0, 0, i);
    var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
    var aMeans = means.slice(0, 0, i);
    var aMean = aMeans.multiply(aCounts)
                      .reduce(ee.Reducer.sum(), [0]).get([0])
                      .divide(aCount);
    var bCount = total.subtract(aCount);
    var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
    return aCount.multiply(aMean.subtract(mean).pow(2)).add(
           bCount.multiply(bMean.subtract(mean).pow(2)));
  });
  return means.sort(bss).get([-1]);
}

4、水体面积以及图表制作

var waterImgCol=joinCol.map(function(image) {
  var his = image.reduceRegion({
                    reducer: ee.Reducer.histogram(), 
                    geometry: roi,
                    scale: 100,
                    maxPixels: 1e13,
                    tileScale: 16
                  });
  var threshold = OTSU(ee.Dictionary(his).get("VV"));
  var water = image.lte(threshold);
  water = water.set("date", image.get("date"));
  water = water.set("system:time_start", image.get("system:time_start"));
  water = water.updateMask(water).rename("water");
  var mask = water.connectedPixelCount(10, true);
  water = water.updateMask(mask.gte(10));
  var area = ee.Image.pixelArea().multiply(water);
  var dict = area.reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: roi,
    scale: 50,
    maxPixels: 1e13, 
    tileScale: 16
  });
  return water.set("area", ee.Number(dict.get("area")).divide(1000000));
});
print("waterImgCol", waterImgCol);
var dataList = waterImgCol.reduceColumns(ee.Reducer.toList(2), ["date", "area"])
                          .get("list");
dataList.evaluate(function(dList) {
  var yValues = [];
  var xValues = [];
  for (var i=0; i<dList.length; i++) {
    var data = dList[i];
    xValues.push(data[0]+"");
    yValues.push(data[1]);
  }
  
  var chart = ui.Chart.array.values(ee.List(yValues), 0, ee.List(xValues))
                .setSeriesNames(["水域面积"])
                .setOptions({
                  title: "水域面积变化趋势", 
                  hAxis: {title: "Date"},
                  vAxis: {title: "水域面积 (km^2)"},
                  legend: null,
                  lineWidth:1,
                  pointSize:2
                });
  print(chart);
});

var params = {
  framesPerSecond: 1,
  region: roi,
  min:-25,
  max:5,
  dimensions: 720,
};
print(ui.Thumbnail(joinCol, params));

var params = {
  framesPerSecond: 1,
  region: roi,
  palette: "blue",
  dimensions: 720,
};
 print(ui.Thumbnail(waterImgCol, params));

---

总结

本章内容基于GEE平台，利用哨兵卫星光学数据和雷达数据分别精确提取河南特大暴雨区域的水体动态变化和受灾分布情况。

我是沐风，致力于打造丝滑的知识分享，集勤劳睿智帅气于一体的沐风。

期待下期精彩案例分享，欢迎交流学习。

温馨提示：本栏目将于明年六月份转为付费专栏，所有盈利将捐献山区希望小学工程。

标签：水灾,get,ee,image,water,start,GEE,var,Scholars
来源： https://blog.csdn.net/willpeng0258/article/details/121851502