今天学习的是如何利用GEE自动提取河宽。
如果大家手里有更多关于GEE有趣的应用，欢迎分享给我！
本文涵盖的研究思路和代码均来自文章RivWidthCloud: An Automated Google Earth Engine Algorithm for River Width Extraction From Remotely Sensed Imagery！欢迎大家去阅读学习！
目录
3.2.1 准备包含L8,5,7的影像集合（注意L7的影像只到03年）
一.代码调用
文章给的代码有JS和Python两个版本，这里只讲述JS版本，代码均可见于https://github.com/seanyx/RivWidthCloudPaper
样本示例
var fns = require('users/eeProject/RivWidthCloudPaper:rwc_landsat.js'); //引用大神的库
var imageId = "LC08_L1TP_022034_20130422_20170310_01_T1"; //以一张L8影像为例
var aoi = ee.Geometry.Polygon(
[[[-88.47207053763748, 37.46382559855354],
[-88.47207053763748, 37.375480838211885],
[-88.2592104302156, 37.375480838211885],
[-88.2592104302156, 37.46382559855354]]], null, false);
var rwc = fns.rwGenSR('Jones2019', 4000, 333, 500, aoi); //定义方法参数
var widths = rwc(imageId); //提取图像中河流宽度
// widths为河流每个点所在河段的宽度，为点集，大致表现如下图所示
widths = widths.map(function(f) {return(f.setGeometry(null))});
//删除点集的空间属性并导出为csv
Export.table.toDrive({
collection: widths,
description: imageId,
folder: "",
fileNamePrefix: imageId,
fileFormat: "CSV"});大致效果如下图所示（注意以上代码只是导出为csv文件，想查看如下图效果需要加上Map.addLayer(widths)
关于如何设置参数
rwGenSR(WATER_METHOD, MAXDISTANCE, FILL_SIZE, MAXDISTANCE_BRANCH_REMOVAL, AOI)- WATER_METHOD {String} 水体提取方法，这里有两个方法可选，'Jones2019' 和 'Zou2018'，该参数默认为'Zou2018'，注：两种方法文章来源都可在作者原文引文中找到，这里不再赘述。
- MAXDISTANCE {integer，单位：米}，即河流河岸到河流中心线的最大值，默认为4000m
- FILL_SIZE {integer,单位：像素}，河流中有像素数小于该阈值的异物（岛屿、堤坝）在计算河流中心线前会被去除，用于避免不必要的分段，默认为333像素
- MAXDISTANCE_BRANCH_REMOVAL {integer，单位：像素}，修建小于这一长度的河流分支，默认为500像素
-
AOI {ee.Geometry.Polygon} 研究区范围
关于输出的结果
| 列名 | 描述 | 单位 |
| 纬度 | 点所在纬度 | |
| 经度 | 点所在经度 | |
| 宽度 | 点所在河段宽度 | 米 |
| 流向的正交方向 | 点所在河流横断面的正交角度 | 弧度 |
| 海拔 | 基于MERIT DEM的整个河面的平均高程 | 米 |
| 图像id | 输入的Landsat图像的ID | 无单位 |
| 投影 | 图像投影 | 无单位 |
| 山影干扰程度 | 范围：0-1 （0表示毫无影响） | 无单位 |
| 冰雪干扰程度 | 范围：0-1 （0表示毫无影响） | 无单位 |
| 云干扰程度 | 范围：0-1 （0表示毫无影响） | 无单位 |
| 云层阴影 | 范围：0-1 （0表示毫无影响） | 无单位 |
| 河流宽度的横断面长度不足而导致宽度测量不准确的程度 | 范围：0-1 （0表示毫无影响） | 无单位 |
| 离图像边缘太近导致宽度测量不准确的程度 | 范围：0-1 （0表示毫无影响） | 无单位 |
二.方法原理
方法主要可分成四个部分，1.提取河流 2.提取河流中心线 3.提取每点处河流流向 4.按照流向正交方向计算河流宽度（原文里有更详细的描述，这里不再赘述）
三.代码详解
细心的小伙伴肯定发现了，该代码需要输入Landsat的图像ID，对于使用哨兵影像或我们不想使用文中的 提取水方法，我们该如何灵活的调用代码来达成所需。为此，我们需要对代码进行深刻的剖析。时常梳理大佬们的代码，能加强我们对GEE的灵活运用能力。
预警，这一部分会比较啰嗦。
3.1 首先我们要先了解作者各个部分代码的整体内容
- functions_landsat.js landsat预处理代码
- functions_waterClassification_Jones2019.js 提取水方法代码1
- functions_waterClassification_Zou2018.js 提取水代码2
-
functions_river.js 河流提取代码
-
functions_centerline_width.js 河流宽度提取代码
代码讲解不是按照代码运行顺序而是按照逻辑顺序，请注意！
3.2 影像预处理
3.2.1 准备包含L8,5,7的影像集合（注意L7的影像只到03年）
exports.merge_collections_std_bandnames_collection1tier1_sr = function() {
// merge landsat 5, 7, 8 collection 1 tier 1 imageCollections and standardize band names
// refer and define standard band names
var bn8 = ['B1', 'B2', 'B3', 'B4', 'B6', 'pixel_qa', 'B5', 'B7'];
var bn7 = ['B1', 'B1', 'B2', 'B3', 'B5', 'pixel_qa', 'B4', 'B7'];
var bn5 = ['B1', 'B1', 'B2', 'B3', 'B5', 'pixel_qa', 'B4', 'B7'];
var bns = ['uBlue', 'Blue', 'Green', 'Red', 'Swir1', 'BQA', 'Nir', 'Swir2'];
// create a merged collection from landsat 5, 7, and 8
var ls5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR").select(bn5, bns);
var ls7 = (ee.ImageCollection("LANDSAT/LE07/C01/T1_SR")
.filterDate('1999-04-15', '2003-05-30') // exclude LE7 images that affected by failure of Scan Line Corrector (SLC)
.select(bn7, bns));
var ls8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR").select(bn8, bns);
var merged = ls5.merge(ls7).merge(ls8);
return(merged);
};
3.2.2 按照给出的ID从影像集合中检索影像
exports.id2Img = function(id) {
return(ee.Image(exports.merge_collections_std_bandnames_collection1tier1_sr()
.filterMetadata('LANDSAT_ID', 'equals', id)
.first()));
};3.2.3 评估影像质量（云、冰、影。。）
以下类似代码在landsat预处理中经常遇得到，但是还是想细致的讲一下
var Unpack = function(qualityBand, startingBit, bitWidth) {
return(qualityBand
.rightShift(startingBit)
.bitwiseAnd(ee.Number(2).pow(bitWidth).subtract(1).int()));
}; //解压，原始landsat的BQA波段是将二进制压缩为十进制，这里是转换回来，并找到二进制上第（startingBit）位的数字，如果其值是（bitWidth）就为1，否则为0，所有像素按照此生成一张图像
var UnpackAllSR = function(bitBand) {
//按照官方的质量参数运用上述函数分别将云、冰雪等影像因素找出来，并合成一个多波段影像。举例：将BQA波段十进制转化为二进制后，像素值第五位数字是1的话该像素就为云，配合Unpack函数生成值为云的一张图像，最后和其他冰雪等图像合成一张多波段影像
var bitInfoSR = {
'Cloud': [5, 1],
'CloudShadow': [3, 1],
'SnowIce': [4, 1],
'Water': [2, 1]
};
var unpackedImage = ee.Image();
for (var key in bitInfoSR) {
unpackedImage = ee.Image.cat(
unpackedImage, Unpack(bitBand, bitInfoSR[key][0], bitInfoSR[key][1])
.rename(key));
}
return(unpackedImage.select(Object.keys(bitInfoSR)));
};
exports.AddFmaskSR = function(image) {
var temp = UnpackAllSR(image.select(['BQA']));
var fmask = (temp.select(['Water']).rename(['fmask'])
// 生成一张包含所有云、冰雪影响因素的图像fmask,按照水为1，云为4，冰雪为3等赋值
.where(temp.select(['SnowIce']), ee.Image(3))
.where(temp.select(['CloudShadow']), ee.Image(2))
.where(temp.select(['Cloud']), ee.Image(4)))
.mask(temp.select(['Cloud']).gte(0)); //保持和fmask一样的掩膜
return(image.addBands(fmask)); //最后得到影响因素的影像
};3.2.4 评估山影影响
这里作者用的dem数据为MERIT，GEE里也有其他的计算山影代码，这里按下不表
exports.CalcHillShadowSR = function(image) {
var dem = ee.Image("users/eeProject/MERIT").clip(image.geometry().buffer(9000).bounds());
var SOLAR_AZIMUTH_ANGLE = ee.Number(image.get('SOLAR_AZIMUTH_ANGLE'));
var SOLAR_ZENITH_ANGLE = ee.Number(image.get('SOLAR_ZENITH_ANGLE'));
return(ee.Terrain.hillShadow(dem, SOLAR_AZIMUTH_ANGLE, SOLAR_ZENITH_ANGLE, 100, true).rename(['hillshadow']));
};3.2.5 用于分类水的函数代码
这里是作者调用写有方法的其他js文件进行提取水，这里就不再继续解读提取水的方法了，感兴趣的同学可以去看看相关的文章和代码（大家也可以尝试此方法自己提取水）例如：
waterJones2019.ClassifyWaterJones2019(你的影像) （记得调用库）
exports.ClassifyWater = function(imgIn, method) {
if (method == 'Jones2019') {
var waterJones2019 = require('users/eeProject/RivWidthCloudPaper:functions_Landsat578/functions_waterClassification_Jones2019.js');
return(waterJones2019.ClassifyWaterJones2019(imgIn));
} else if (method == 'Zou2018') {
var waterZou2018 = require('users/eeProject/RivWidthCloudPaper:functions_Landsat578/functions_waterClassification_Zou2018.js');
return(waterZou2018.ClassifyWaterZou2018(imgIn));
}
};3.2.6 分类提取水并加入之前的质量波段
exports.CalculateWaterAddFlagsSR = function(imgIn, waterMethod) {
waterMethod = typeof waterMethod !== 'undefined' ? waterMethod : 'Jones2019';
//加入质量波段
var fmask = exports.AddFmaskSR(imgIn).select(['fmask']);
//将质量波段云、冰雪单独提取出来，用于后期判断其不同的影响
var fmaskUnpacked = fmask.eq(4).rename('flag_cloud')
.addBands(fmask.eq(2).rename('flag_cldShadow'))
.addBands(fmask.eq(3).rename('flag_snowIce'))
.addBands(fmask.eq(1).rename('flag_water'));
//对水进行分类
var water = exports.ClassifyWater(imgIn, waterMethod)
.where(fmask.gte(2), ee.Image.constant(0));
var hillshadow = exports.CalcHillShadowSR(imgIn).not().rename(['flag_hillshadow']);
var imgOut = ee.Image(water.addBands(fmask).addBands(hillshadow).addBands(fmaskUnpacked)
.setMulti({
'image_id': imgIn.get('LANDSAT_ID'),
'timestamp': imgIn.get('system:time_start'),
'scale': imgIn.projection().nominalScale(),
'crs': imgIn.projection().crs()
}));
// 设置影像的属性，包括ID，时间，分辨率，投影系
return(imgOut);
};3.3 提取河流
这一步的输入参数有四个，包括上一步得到的水提取图像、GREL、你定义或者使用默认参数的MAXDISTANCE和FILL_SIZE。
关于GREL参数，作者使用的是GRWL矢量产品V01.01数据集，即Global River Widths from Landsat，是一个全球河流数据库，包含5800万个河流中心线位置和宽度。（数据集链接：https://zenodo.org/record/1297434#.W8JkshNKh24）
var grwl = ee.FeatureCollection("users/eeProject/grwl");之前我以为作者是按照侵蚀、膨胀等形态学操作来提取河流，现在看来我想错了。
3.3.1 先依照已有数据得到河流大致中心线
筛选水体的范围内（bound) 的GREL数据（clDataset）
exports.GetCenterline = function(clDataset, bound) {
var cl = clDataset.filterBounds(bound);
return(cl);
};3.3.2 提取河流通道
image为水体的二值图像，centerline为上步得到的大致中心线数据，maxDistance为河岸到河流最大距离。
这一步主要是计算距离河流中心线的累计成本，关于累计成本计算方式可参加ArcGIS教程：成本距离 (空间分析)
exports.ExtractChannel = function(image, centerline, maxDistance) {
//下局image.not()为成本函数，这里水是0，陆地是1
var connectedToCl = image.not().cumulativeCost({
//源图像，先是将centerline转变为像素值为1的栅格图像，然后只保留水体内的水体中心线栅格信息
source: ee.Image().toByte().paint(centerline, 1).and(image),
//距离超过多少米不予考虑
maxDistance: maxDistance,
//false计算平面距离，true考虑大地曲面
geodeticDistance: false
}).eq(0);
//上文计算了中心线以外到河流中心线的累计成本，由于陆地的总成本将超过0（陆地成本函数为1），河流的总成本为0，因此最后选了eq(0)来提取河流
var channel = image.updateMask(connectedToCl).unmask(0).updateMask(image.gte(0)).rename(['channelMask']);
//将河流以外的值取值为0，并和输入的Image保持相同的掩膜范围
return channel;
};3.3.3 删除细小岛屿、堤坝等
这一步就是按照你之前给定的面积将面积小于该阈值的非水体变成水体，我之前关于GEE面向对象的部分有涉及。
exports.RemoveIsland = function(channel, FILL_SIZE) {
var fill = channel.not().selfMask().connectedPixelCount(FILL_SIZE).lt(FILL_SIZE);
var river = channel.where(fill, ee.Image(1)).rename(['riverMask']);
return(river);
};3.3.4 以上三部的总体汇总
相信在我讲解完所有函数后，大家就能明白下面各个步骤的含义了
exports.ExtractRiver = function(imgIn, clData, maxDist, minIslandRemoval) {
var waterMask = imgIn.select('waterMask');
var bound = waterMask.geometry();
var cl = exports.GetCenterline(clData, bound);
var channelMask = exports.ExtractChannel(waterMask, cl, maxDist);
var riverMask = exports.RemoveIsland(channelMask, minIslandRemoval);
return(imgIn.addBands(channelMask).addBands(riverMask));
};3.4 提取河流中心线
3.4.1 计算河流中像素到河岸的距离
输入参数为统计区域邻域大小neighborhoodSize和空间分辨率scale
exports.CalcDistanceMap = function(img, neighborhoodSize, scale) {
//设置focal_max的迭代次数为1，2，在二值图像中即做两次宽度为1，2的膨胀处理，两者相减即得到了河岸的边界
var imgD2 = img.focal_max(1.5, 'circle', 'pixels', 2);
var imgD1 = img.focal_max(1.5, 'circle', 'pixels', 1);
var outline = imgD2.subtract(imgD1);
//计算河流像素到边界的距离
var dpixel = outline.fastDistanceTransform(neighborhoodSize).sqrt();
//单位由像素换算成米
var dmeters = dpixel.multiply(scale);
var DM = dmeters.mask(dpixel.lte(neighborhoodSize).and(imgD2));
//掩膜得到河流部分，输出的图像值为该河流像素到河岸的距离
return(DM);
};3.4.2 计算到河岸距离的梯度图像
由于河流中心线到河岸的距离要比周围一圈水体像素大，因此需对距离图像求梯度方法来提取中心线（卷积方法，这个较基础在此不表）。函数里设置了针对不同河宽的求梯度方法，这里默认为9*9像素范围的方法2。
exports.CalcGradientMap = function(image, gradMethod, scale) {
var dx, dy, g, k_dx, k_dy;
if (gradMethod == 1) {
var grad = image.gradient();
dx = grad.select(['x']);
dy = grad.select(['y']);
g = dx.multiply(dx).add(dy.multiply(dy)).sqrt();
}
if (gradMethod == 2) {
k_dx = ee.Kernel.fixed(3, 3,
[[ 1/8, 0, -1/8],
[ 2/8, 0, -2/8],
[ 1/8, 0, -1/8]]);
k_dy = ee.Kernel.fixed(3, 3,
[[ -1/8, -2/8, -1/8],
[ 0, 0, 0],
[ 1/8, 2/8, 1/8]]);
dx = image.convolve(k_dx);
dy = image.convolve(k_dy);
g = dx.multiply(dx).add(dy.multiply(dy)).divide(scale.multiply(scale)).sqrt();
}
if (gradMethod == 3) {
k_dx = ee.Kernel.fixed(3, 1,
[[-0.5, 0, 0.5]]);
k_dy = ee.Kernel.fixed(1, 3,
[[0.5],
[0],
[-0.5]]);
dx = image.convolve(k_dx);
dy = image.convolve(k_dy);
g = dx.multiply(dx).add(dy.multiply(dy)).divide(scale.multiply(scale));
}
return(g);
};3.4.3 提取河流骨架
提取骨架涉及到更为复杂的形态学操作，这里给出两篇参考资料
Hit-and-Miss Transform 和 Donchyts, 2016
var HitOrMiss = function(image, se1, se2) {
var e1 = image.reduceNeighborhood(ee.Reducer.min(), se1);
var e2 = image.not().reduceNeighborhood(ee.Reducer.min(), se2);
return(e1.and(e2));
};
var SplitKernel = function(kernel, value) { //kernl里值为value的返回1否则返回0
var result = [];
for(var r = 0; r < kernel.length; r++) {
var row = [];
for(var c = 0; c < kernel.length; c++) {
row.push(kernel[r][c] == value ? 1 : 0);
}
result.push(row);
}
return(result);
};
exports.Skeletonize = function(image, iterations, method) { //本文默认方法1
var se1w = [[2, 2, 2],
[0, 1, 0],
[1, 1, 1]];
if(method == 2) {
se1w = [[2, 2, 2],
[0, 1, 0],
[0, 1, 0]];
}
var se11 = ee.Kernel.fixed(3, 3, SplitKernel(se1w, 1));
var se12 = ee.Kernel.fixed(3, 3, SplitKernel(se1w, 2));
var se2w = [[2, 2, 0],
[2, 1, 1],
[0, 1, 0]];
if(method == 2) {
se2w = [[2, 2, 0],
[2, 1, 1],
[0, 1, 1]];
}
var se21 = ee.Kernel.fixed(3, 3, SplitKernel(se2w, 1));
var se22 = ee.Kernel.fixed(3, 3, SplitKernel(se2w, 2));
var result = image;
for(var i = 0; i < iterations; i++) {
for(var j=0; j<4; j++) { // 对kernel进行四轮旋转
result = result.subtract(HitOrMiss(result, se11, se12));
se11 = se11.rotate(1);
se12 = se12.rotate(1);
result = result.subtract(HitOrMiss(result, se21, se22));
se21 = se21.rotate(1);
se22 = se22.rotate(1);
}
}
return(result.rename(['clRaw']));
};
调用上述方法，得到骨架图像
exports.CalcOnePixelWidthCenterline = function(img, GM, hGrad) {
var imgD2 = img.focal_max(1.5, 'circle', 'pixels', 2);
var cl = ee.Image(GM).mask(imgD2).lte(hGrad).and(img);
// 迭代skeletonization方法2次
var cl1px = exports.Skeletonize(cl, 2, 1);
return(cl1px);
};3.4.4 优化得到的中心线（骨架）
关于优化方法的函数，方法也是来自于上步的形态学方法。
下面是计算中心线端点和拐点两个函数，当一条线点过多时，它们就会成为冗余点，不利于我们提取1像素宽的中心线
var ExtractEndpoints = function(CL1px) { //计算中心线的端点
var se1w = [[0, 0, 0],
[2, 1, 2],
[2, 2, 2]];
var se11 = ee.Kernel.fixed(3, 3, SplitKernel(se1w, 1));
var se12 = ee.Kernel.fixed(3, 3, SplitKernel(se1w, 2));
var result = CL1px;
for(var i=0; i<4; i++) {
result = result.subtract(HitOrMiss(result, se11, se12));
se11 = se11.rotate(1);
se12 = se12.rotate(1);
}
var endpoints = CL1px.subtract(result);
return(endpoints);
};
var ExtractCorners = function(CL1px) { //计算中心线的拐点
var se1w = [[2, 2, 0],
[2, 1, 1],
[0, 1, 0]];
var se11 = ee.Kernel.fixed(3, 3, SplitKernel(se1w, 1));
var se12 = ee.Kernel.fixed(3, 3, SplitKernel(se1w, 2));
var result = CL1px;
for(var i=0; i<4; i++) {
result = result.subtract(HitOrMiss(result, se11, se12));
se11 = se11.rotate(1);
se12 = se12.rotate(1);
}
var cornerPoints = CL1px.subtract(result);
return(cornerPoints);
};下面为调用上述函数对中心线进行精简
exports.CleanCenterline = function(cl1px, maxBranchLengthToRemove, rmCorners) {
//计算点密度
var nearbyPoints = cl1px.mask(cl1px).reduceNeighborhood({
reducer: ee.Reducer.count(),
kernel: ee.Kernel.circle(1.5),
skipMasked: true});
//点密度小的像素点，端点
var endsByNeighbors = nearbyPoints.lte(2);
//点密度大的像素，一般为支流相交点
var joints = nearbyPoints.gte(4);
//进行必要掩膜后，用累计成本函数计算到支流端点距离（支线长度），较短长度的支流被保留下来，其余被掩膜掉，这里阈值是输入的参数maxBranchLengthToRemove
var costMap = cl1px.mask(cl1px).updateMask(joints.not()).cumulativeCost({
source: endsByNeighbors.mask(endsByNeighbors),
maxDistance: maxBranchLengthToRemove,
geodeticDistance: false});
var branchMask = costMap.gte(0).unmask(0);
var cl1Cleaned = cl1px.updateMask(branchMask.not()); //去除小于阈值的支流
var ends = ExtractEndpoints(cl1Cleaned); //去除冗余的端点（这里冗余的端点指什么我还真不清楚。。）
cl1Cleaned = cl1Cleaned.updateMask(ends.not());
if (rmCorners) {
var corners = ExtractCorners(cl1Cleaned);
cl1Cleaned = cl1Cleaned.updateMask(corners.not()); //去除多于的拐点（拐弯处形成直角的多余像素）
}
return(cl1Cleaned);
};3.4.5 汇总以上函数提取河流中心线
这里CleanCenterline函数迭代了两次
exports.CalculateCenterline = function(imgIn) {
var crs = imgIn.get('crs');
var scale = ee.Number(imgIn.get('scale'));
var riverMask = imgIn.select('riverMask');
var distM = exports.CalcDistanceMap(riverMask, 256, scale);
var gradM = exports.CalcGradientMap(distM, 2, scale);
var cl1 = exports.CalcOnePixelWidthCenterline(riverMask, gradM, 0.9);
var cl1Cleaned1 = exports.CleanCenterline(cl1, 300, true);
var cl1px = exports.CleanCenterline(cl1Cleaned1, 300, false);
var imgOut = imgIn
.addBands(cl1px.rename('cleanedCL'))
.addBands(cl1.rename('rawCL'))
.addBands(gradM.rename('gradientMap'))
.addBands(distM.rename('distanceMap'));
return(imgOut);
};
3.5 提取河流宽度
3.5.1 计算河流横截面角度
这里是建立了9*9的卷积框架，外围像素值为角度，通过外围像素和河流中位线相交的两个角度值来推算出河流的流向的正交角度（横截面角度）。最后得到中位线正交角度图像。
exports.CalculateAngle = function(clCleaned) {
//定义卷积模板
var w3 = (ee.Kernel.fixed(9, 9, [
[135.0, 126.9, 116.6, 104.0, 90.0, 76.0, 63.4, 53.1, 45.0],
[143.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 36.9],
[153.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 26.6],
[166.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 14.0],
[180.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1e-5],
[194.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 346.0],
[206.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 333.4],
[216.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 323.1],
[225.0, 233.1, 243.4, 256.0, 270.0, 284.0, 296.6, 306.9, 315.0]]));
var combinedReducer = ee.Reducer.sum().combine(ee.Reducer.count(), null, true);
//统计外围像素和河流中位数相交的数量和求和
var clAngle = (clCleaned.mask(clCleaned)
.rename(['clCleaned'])
.reduceNeighborhood({
reducer: combinedReducer,
kernel: w3,
inputWeight: 'kernel',
skipMasked: true}));
//相交数量大于2的像素被掩膜掉
var clAngleNorm = (clAngle
.select('clCleaned_sum')
.divide(clAngle.select('clCleaned_count'))
.mask(clAngle.select('clCleaned_count').gt(2).not()));
//当只有一个像素相交时，角度加90变成正交值
clAngleNorm = (clAngleNorm
.where(clAngle.select('clCleaned_count').eq(1), clAngleNorm.add(ee.Image(90))));
return(clAngleNorm.rename(['orthDegree']));
}; 
3.5.2 计算河流宽度
先将中位线角度图像变成具有角度属性的点集，为下面遍历每个点做准备（这里点已经从波段中汲取了到岸边距离的属性
var clPoints = (clAngleNorm.rename(['angle'])
.addBands(ee.Image.pixelCoordinates(crs))
.addBands(ee.Image.pixelLonLat().rename(['lon', 'lat']))
.addBands(DM.rename(['toBankDistance']))
.sample({
region: bound,
scale: scale,
projection: null,
factor: 1,
dropNulls: true
}));根据角度画出一条线段，和河岸的两个交点即为测量河流宽度的两个端点。（因为统一定义线段的长度，太长太短都不好，于是画出等于到岸边距离3倍的一条的线段，即上图红色Wm
var GetXsectionEnds = function(f) {
var xc = ee.Number(f.get('x'));
var yc = ee.Number(f.get('y'));
var orthRad = ee.Number(f.get('angle')).divide(180).multiply(Math.PI); //换算成弧度
var halfWidth = ee.Number(f.get('toBankDistance')).multiply(1.5);
//这里将到岸边的距离扩大到1.5倍
var cosRad = halfWidth.multiply(orthRad.cos());
var sinRad = halfWidth.multiply(orthRad.sin());
var p1 = ee.Geometry.Point([xc.add(cosRad), yc.add(sinRad)], crs);
var p2 = ee.Geometry.Point([xc.subtract(cosRad), yc.subtract(sinRad)], crs);
var xlEnds = (ee.Feature(ee.Geometry.MultiPoint([p1, p2]), {
'xc': xc,
'yc': yc,
'longitude': f.get('lon'),
'latitude': f.get('lat'),
'orthogonalDirection': orthRad,
'MLength': halfWidth.multiply(2), //距离岸边距离的三倍
'p1': p1,
'p2': p2,
'crs': crs,
'image_id': sceneID
}));
return(xlEnds);
};
var xsectionsEnds = clPoints.map(GetXsectionEnds); 遍历计算线段Wm通过和之前提取河流的图像相乘，得到河流内的线段，并统计线段内的值来检验线段端点是否在水中，或者超出给定边界。最后统计根据河流二值图像统计线段Wm的均值，利用均值乘Wm总长度得到最后的宽度。
下面代码中输入参数endInfo为去除小岛屿等杂物后的河流二值图像，endStat为河流二值图像。
var xsectionsEnds = clPoints.map(GetXsectionEnds);
var endStat = (endInfo.reduceRegions({
collection: xsectionsEnds,
reducer: ee.Reducer.anyNonZero().combine(ee.Reducer.count(), null, true), //# test endpoints type 1. if in water or 2. if extends over the image bound
scale: scale,
crs: crs}));
// ## calculate the width of the river and other flags along the xsection lines
var xsections1 = endStat.map(SwitchGeometry);
var combinedReducer = ee.Reducer.mean();
var xsections = (segmentInfo.reduceRegions({
collection: xsections1,
reducer: combinedReducer,
scale: scale,
crs: crs}));
最后结合上述函数导出结果
exports.prepExport = function(f) {
f = f.set({
'width': ee.Number(f.get('MLength')).multiply(f.get('channelMask')),
'endsInWater': ee.Number(f.get('any')).eq(1), //提取的河宽线段是否在水中，是为1
'endsOverEdge': ee.Number(f.get('count')).lt(2) //提取的河宽线段是否会超出给定边界，是为1
});
var fOut = ee.Feature(ee.Geometry.Point([f.get('longitude'), f.get('latitude')]), {})
.copyProperties(f, null, ['any', 'count', 'MLength', 'xc', 'yc', 'channelMask']);
return(fOut);
};
exports.CalculateWidth = function(imgIn) { //收集之前各种属性导出
var crs = imgIn.get('crs');
var scale = imgIn.get('scale');
var imgId = imgIn.get('image_id');
var bound = imgIn.select('riverMask').geometry();
var angle = imgIn.select('orthDegree');
var dem = ee.Image("users/eeProject/MERIT");
var infoEnds = imgIn.select('riverMask');
var infoExport = imgIn.select('channelMask')
.addBands(imgIn.select('^flag.*'))
.addBands(dem.rename('flag_elevation'));
var dm = imgIn.select('distanceMap');
var widths = exports.GetWidth(angle, infoExport, infoEnds, dm, crs, bound, scale, imgId)
.map(exports.prepExport);
return(widths);
};3.6 整合成最终函数
var AssignDefault = function(x, dv) { //当x没定义时 ，一律选择默认值dv
return(typeof x !== 'undefined' ? x : dv);
};
exports.rwGenSR = function(aoi, WATER_METHOD, MAXDISTANCE, FILL_SIZE, MAXDISTANCE_BRANCH_REMOVAL) {
// 分配默认值
WATER_METHOD = AssignDefault(WATER_METHOD, 'Jones2019');
MAXDISTANCE = AssignDefault(MAXDISTANCE, 4000);
FILL_SIZE = AssignDefault(FILL_SIZE, 333);
MAXDISTANCE_BRANCH_REMOVAL = AssignDefault(MAXDISTANCE_BRANCH_REMOVAL, 500);
aoi = AssignDefault(aoi, null);
var grwl = ee.FeatureCollection("users/eeProject/grwl");
var lsFun = require('users/eeProject/RivWidthCloudPaper:functions_Landsat578/functions_landsat.js');
var riverFun = require('users/eeProject/RivWidthCloudPaper:functions_river.js');
var clWidthFun = require('users/eeProject/RivWidthCloudPaper:functions_centerline_width.js');
var tempFUN = function(image) {
aoi = ee.Algorithms.If(aoi, aoi, image.geometry());
image = image.clip(aoi);
// 提取水
var imgOut = lsFun.CalculateWaterAddFlagsSR(image, WATER_METHOD);
// 提取河流
imgOut = riverFun.ExtractRiver(imgOut, grwl, MAXDISTANCE, FILL_SIZE);
// 计算中心线
imgOut = clWidthFun.CalculateCenterline(imgOut);
// 计算正交方向
imgOut = clWidthFun.CalculateOrthAngle(imgOut);
// 计算宽度
var widthOut = clWidthFun.CalculateWidth(imgOut);
return(widthOut);
};
return(tempFUN);
};到此结束，以后有空的话会整理整理地表温度反演的代码~
本次整理比较仓促，如有不足欢迎私信或评论指点！！！
2120
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
