使用GEE进行岸线自动提取

最新推荐文章于 2023-07-15 19:13:19 发布

GEE水生态空间

最新推荐文章于 2023-07-15 19:13:19 发布

阅读量3.7k

点赞数 12

文章标签： google earth

本文链接：https://blog.csdn.net/qq_32497109/article/details/129880050

版权

使用GEE进行岸线自动提取

2023 年 1 月 18 日

在本文中，我将概述一种在 Google Earth Engine 中从卫星图像中提取海岸线的方法。此方法是可扩展的，并自动将海岸线提取为矢量折线。完整的代码链接在帖子末尾可用。

该方法包括以下步骤

创建无云合成图像
提取所有水体
移除内陆水域和小岛屿
将栅格转换为矢量
简化和提取海岸线

我们将详细介绍每个步骤并查看实现结果所需的 Google Earth Engine API 代码。

创建无云合成图像

为了能够正确检测海岸线，我们需要无云图像。推荐的方法是创建合成图像。您可以选择合适的日期范围（每月/每年），应用云遮罩去除多云像素，然后创建中值合成以获得无云图像。

提示：多云地区

如果您在一个非常多云的地区工作，中值合成可能仍然有云。在这种情况下，尝试创建一个百分位合成（即 25 个百分位合成），这会给您带来更好的结果。在EEFA Book 的 F4.1 章中了解如何创建百分位数组合。

您可以使用任何多光谱传感器来提取海岸线。至少有一个 NIR 波段，最好是一个 SWIR 波段，以可靠地检测水是有用的。Landsat 和 Sentinel-2 是此任务中使用最广泛的数据集。

在此示例中，我们从 Landsat Level 2（表面反射率）收集开始，过滤选定年份在选定区域收集的图像，应用云遮罩，然后创建中值合成。

// Define parameters

var startDate = ee.Date('2020-01-01');

var endDate = startDate.advance(1, 'year');

var collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')

.filter(ee.Filter.date(startDate, endDate))

.filter(ee.Filter.bounds(geometry));

// Apply cloud mask

var maskL8sr = function(image) {

var qaMask = image.select('QA_PIXEL')

.bitwiseAnd(parseInt('11111', 2)).eq(0);

var saturationMask = image.select('QA_RADSAT').eq(0);

// Apply the scaling factors to the appropriate bands.

var opticalBands = image.select('SR_B.')

.multiply(0.0000275).add(-0.2);

var thermalBands = image.select('ST_B.*')

.multiply(0.00341802).add(149.0);

// Replace the original bands with the scaled ones

// and apply the masks.

return image.addBands(opticalBands, null, true)

.addBands(thermalBands, null, true)

.updateMask(qaMask)

.updateMask(saturationMask);

};

var collection = collection.map(maskL8sr);

// Select and Rename Bands

var collection = collection.select(

['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'],

['blue', 'green', 'red', 'nir', 'swir1', 'swir2']

);

// Create a median composite and clip

var composite = collection.median().clip(geometry);

var rgbVis = {

bands: ['red', 'green', 'blue'],

min: 0.0,

max: 0.3

};

Map.addLayer(composite, rgbVis, 'Composite Image');

警告：处理潮汐高度

一年中所有图像的中值合成将代表所有观察结果的平均海岸线位置——包括涨潮和退潮。这种复合材料不适合比较多年来海岸线的变化。要比较随时间的变化——您必须确保图像处于潮汐的同一阶段。您可以阅读这篇论文，了解一种屏蔽高潮和低潮像素的方法。

另一种方法是使用潮汐估计并过滤在同一潮汐阶段捕获的图像。来源将取决于感兴趣的国家和地区。

– HYCOM surface_elevation数据集可能对推导潮汐信息有用

如果您知道任何其他数据集，请在评论中告诉我。

提取所有水体

有几种方法可以从卫星图像中检测和提取地表水。我们想要一种无需手动识别阈值的自动化方法——因此我们将使用 Otsu 动态阈值。在这里，我们计算自动水提取指数 (AWEI) 并使用 Otsu 阈值找到给定图像中水像素的最佳阈值。您也可以尝试使用其他指数，例如 NDWI 或 MNDWI，它们可能在您感兴趣的地区效果更好。如果这种方法没有给出令人满意的结果，我建议你看看无监督聚类方法，它结合了不同的指标，并使用机器学习方法自动识别水像素。

// Code to compute AWEI and detect water

// using Otsu thresholding

var detectWater = function(image) {

var awei = image.expression(

'4 * (GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2)', {

'GREEN': image.select('green'),

'NIR': image.select('nir'),

'SWIR1': image.select('swir1'),

'SWIR2': image.select('swir2'),

}).rename('awei');

// Otsu Thresholding

var thresholding = require(

'users/gena/packages:thresholding');

var scale = 100;

var bounds = geometry;

var cannyThreshold = 0.7;

var cannySigma = 1;

var minValue = -0.2;

var th = thresholding.computeThresholdUsingOtsu(

awei, scale, bounds,

cannyThreshold, cannySigma, minValue);

// Create a Land-Water Image using Otsu Threshold

// You can replace th with a manual threshold if

// Otsu results are not satisfactory

var water = awei.gt(th).rename('water');

return water;

};

var water = detectWater(composite);

var waterVis = {min:0, max:1, palette: ['white', 'blue']};

Map.addLayer(water, waterVis, 'All Water');

移除内陆水域和小岛屿

我们的主要兴趣是陆水边界，因此我们移除了所有内陆水体和海洋中的小岛。这是通过删除所有小于所选尺寸的水像素簇来完成的。

// Export resolution in meters is at which the coastline will

// be vectorized.

// Higher resolution (such as 10) will give smooother results

var exportResolutionMeters = 30;

// We need a threshold in meters to remove inland waterbodies

// and coastal islands. If you have a large waterbody close

// to the coastline, set this higher.

var waterbodySizeMeters = 2000;

var islandSizeMeters = 1000;

// This function takes a binary Land-Water image and

// removes inland water and small islands

function removeInlandWaterAndIslands(waterImage) {

// reduceConnectedComponents expects an interger image

waterImage = waterImage.int();

// Define neighborhood based on user parameters

var connectedPixelsLand = ee.Number(waterbodySizeMeters)

.divide(exportResolutionMeters).int();

var connectedPixelsWater = ee.Number(islandSizeMeters)

.divide(exportResolutionMeters).int();

// Remove inland water

var landFilled = waterImage.addBands(waterImage)

.reduceConnectedComponents(

ee.Reducer.median(), 'water', connectedPixelsLand)

.unmask(99).eq(99).and(waterImage.neq(0));

// Remove small islands

var waterFilled = landFilled.addBands(landFilled)

.reduceConnectedComponents(

ee.Reducer.median(), 'water_1', connectedPixelsWater)

.unmask(99).eq(99).and(landFilled.neq(1));

// Land-Water Boundary

return waterFilled;

}

var landWaterBoundary = removeInlandWaterAndIslands(water);

var landWaterBoundaryVis = {

min:0,

max:1,

palette: ['blue', 'white']

};

Map.addLayer(landWaterBoundary, landWaterBoundaryVis,

'Land-Water Boundary (Raster)');

将栅格转换为矢量

接下来，我们将二值水陆图像矢量化以获得陆地区域的多边形轮廓。这也可以作为许多类型分析的最终输出。因此，您可以将其导出为 shapefile 并进一步使用。

// Convert the coastline image to vector

var vectors = ee.Image(landWaterBoundary).selfMask()

.reduceToVectors({

geometry: geometry,

scale: exportResolutionMeters,

eightConnected: true,

maxPixels: 1e10,

tileScale: 16

});

Map.addLayer(vectors, {color: 'blue'},

'Land-Water Boundary (Vector)');

简化和提取海岸线

许多应用程序需要生成的几何图形为线而不是多边形。因此，我们采用生成的多边形，对其进行简化并提取坐标。然后我们使用坐标来构造ee.Geometry.MultiLineString()特征。

1个

2个

3个

4个

5个

6个

8个

// This function takes vectorized polygons and

// extracts a polyline

var simplifyAndExtractCoastline = function(vectors){

// Simplify vectors

var processedVectors = vectors.map(function(f) {

var coords = f.geometry()

.simplify({maxError: exportResolutionMeters})

.coordinates();

// Buffer the geometry by a pixel to avoid rasterizing

// the boundary polygon

var bufferDistance = ee.Number(

exportResolutionMeters).multiply(-1);

return f

.setGeometry(

ee.Geometry.MultiLineString(coords)

.intersection(geometry.buffer(bufferDistance)));

});

return processedVectors;

};

var coastlineVector = simplifyAndExtractCoastline(

vectors);

Map.addLayer(coastlineVector, {color: 'red'},

'Coastline (Vector)');

导出结果

GEE 的代码编辑器环境仅限运行 5 分钟的交互式计算。如果您尝试从一个大区域中提取海岸线，您可能会遇到图块或内存错误。这仅仅意味着您无法以交互方式可视化结果。要在这种情况下获得结果，您必须导出结果。导出可以运行更长时间并且可以访问更大的资源池。您可以将结果导出到 Google 云端硬盘，甚至可以导出到您自己的资产存储中。

// Exports

// Exports can run for longer time and have more resources

// If any of your layers time-out or give tile errors,

// Run the exports instead

Export.table.toDrive({

collection: coastlineVector,

description: 'Extracted_Coastline_Vector',

folder: 'earthengine',

fileNamePrefix: 'coastline',

fileFormat: 'SHP'})

Export.image.toDrive({

image: landWaterBoundary,