GEE：栅格计算（GoogleEarthEngine平台）

以下列出几个ArcGIS中常用到的栅格计算器公式
同样，这些栅格计算公式也可以在GEE平台上实现。

需求	公式
将空值赋予某一特定的数值（如100）	CON(ISNULL(“raster”),100,“raster”)
将某个特定的栅格值（100）赋值为0	Con(“raster”==-100,0,“raster”)
将某一特定的栅格值(如100)赋值为空值	setnull(“raster”==100,“raster”)
将栅格数值在某一条件范围内（如大于100）赋值为0	Con(“raster” >100,0,“raster”)
使目标影像中的多个栅格数值赋予无效值	SetNull(("raster == 100)
栅格计算器中的平方和开方运算	（1）Power（raster，2），Power（raster，0.5）
	（2）Squre()，SquareRoot()

---

where()、clamp()和unitScale()

作用1：提取NDVI、DEM以及各个波段在某一区间的值。
作用2：去除影像中的噪声。
作用3：将某一区间的值重新赋值。

updateMask（）——掩膜掉不符合条件的像素

如果要掩膜掉不符合条件的像素请看这篇文章https://blog.csdn.net/qq_35591253/article/details/120649323

---

例一：https://code.earthengine.google.com/7f0920f9719f807107d72ee08dcebace?noload=true

//原始的NDVI栅格图像。（此时的NDVI值在区间[-1,+1]）

//使用.where(NDVI.lt(0.3),0.3)之后的效果
//将NDVI小于0.3的值全部赋值为0.3
.where(NDVI.lt(0.3),0.3)

//使用.clamp(0,1)之后，的NDVI栅格图像。(此时的NDVI值在区间[0,1])
.clamp(0,1)

//将小于0.2的NDVI值赋值为1，大于0.39的赋值为3，在区间[0.2,0.39]的赋值为2
var NDVI = NDVI.where(greenest.lt(0.2),1)
               .where(greenest.gte(0.2).and(greenest.lte(0.39)),2)
               .where(greenest.gt(0.39),3)

//.unitScale(A,B)
//将栅格数据的值类似于归一化操作，只不过将值转换到A和B之间
.unitScale(0,0.5)

---

最大最小值归一化
将值域为【最小值，最大值】映射为【0，1】

//归一化
function Normalization(image){
  var image_d = image.subtract(Min);//MIn，Max分别是最小值和最大值，需要提前计算好填进来
  return image_d.divide(Max-Min).copyProperties(image, image.propertyNames())
}

等比例映射法
将值域为 A 区间等比例映射到 B 区间
举例：将【-0.2，1.0】的值域映射到【0，255】

//等比例映射法
var a = 255-0;
var b = 1-(-0.2)
varscale=a/b;

//将-0.2-1.0的值域映射到0-255
var newNDVI = (NDVINoAbnormalValue. subtract(-0.2)) . multip1y(scale).select( 'NDVI');

---

有时候我们需要选择NDVI值大于0的像素；又比如，计算的指数值太大，可视化的时候其余像素会因为这个最大的误差而错误的显示，甚至被拉伸至0，导致影像不可用。为了将大于这些阈值的误差像素剔除，我们必须使用栅格计算，来处理遥感影像。
那么怎么在GEE平台上做栅格计算呢？

---

我们可以用.where()和.clamp()
具体语法可以参考https://zhuanlan.zhihu.com/p/76343298
这里抽取他的例子。

var DEM = ee.Image("USGS/SRTMGL1_003");

var Land_Cover = ee.Image("ESA/GLOBCOVER_L4_200901_200912_V2_3").select('landcover');

var High_Land = Land_Cover.where(DEM.lt(4000),0)

Map.setCenter(88.96, 33.77,4)

Map.addLayer(DEM.lt(4000))

Map.addLayer(High_Land)

这里的.where()操作就是将DEM小于4000的值全部设置为0。
同ArcGIS里面的栅格计算。

---

var DEM = ee.Image("USGS/SRTMGL1_003");

var DEM_Clamp = DEM.clamp(450,795)

Map.setCenter(106.5662, 29.5589,10)

Map.addLayer(DEM_Clamp)

这里的.clamp()操作就是将DEM的值在450到795区间的抽取出来，形成新的图层。
注意： 如果一个像素值大于795，就给它赋值795，小于450的，都给它赋值450。
其中，450为最小值，795为最大值。

---

var L8_One = ee.Image(ee.ImageCollection("LANDSAT/LC08/C01/T1")
               .filterBounds(ee.Geometry.Point(106.5776, 29.5577)).first());
               
var L8_Unitscale = L8_One.unitScale(0,32767)

print(L8_One,L8_Unitscale)

Map.addLayer(L8_One,		{"bands":["B5","B4","B3"],"max":30000})

Map.addLayer(L8_Unitscale,	{"bands":["B5","B4","B3"],"max":0.92})

这个类似归一化，只不过归一化是将栅格值转化为-1到+1之间，.unitScale()操作是将栅格值转换到0和+1之间。

---

例子代码：

//定义addNDVI函数，该函数的输入是image，然后对image执行NDVI计算，并将结果命名为ndvi
//并将生成的ndvi添加到原有的image中
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

//选择2019年的L8影像
var lc82019 = L8.filterDate('2019-08-01', '2019-8-31');

//将筛选出来的2019年的影像均执行NDVI计算并添加该结果至原影像
var withNDVI = lc82019.map(addNDVI);

//利用qualityMosaic()函数执行NDVI最大值筛选，并用gz的边界进行裁剪
var greenest = withNDVI.qualityMosaic('NDVI').clip(roi);

//对筛选出来的像素，执行NDVI计算（由于筛选出来的是包含所有波段
//所以可以直接输出全波段影像而不仅仅是NDVI）
var ndvi = greenest.normalizedDifference(['B5', 'B4']).rename('bestNDVI');

//添加生成的最优NDVI到图层窗口
Map.addLayer(ndvi,{palette: ['black', 'green', 'red']},'ndvi')
print("over")

//例一
var ndvi1 = ndvi.clamp(0,1)

// //例二
// var ndvi1 = ndvi.where(ndvi.lt(0.3),0.3)
Map.addLayer(ndvi1,{palette: ['black', 'green', 'red']},'ndvi')

---

附加：
剔除异常值
计算和打印剔除异常值以后的最大值和最小值

//剔除异常值
var ndvi1 = NDVI.where(NDVI.lt(-0.2),-0.2);
var NDVINoAbnormalValue = ndvi1.where(ndvi1.gt(1.0),1.0);

//计算和打印剔除异常值以后的最大值和最小值
var NDVINoAbnormalValueMin = (NDVINoAbnormalValue.select('NDVI').reduceRegion({
    reducer: ee.Reducer.min(),
    geometry: geometry,
    scale: 30, 
    crs: "EPSG:4326",
    bestEffort: true,
    tileScale: 5
}));
var NDVINoAbnormalValueMax = (NDVINoAbnormalValue.select('NDVI').reduceRegion({
    reducer: ee.Reducer.max(),
    geometry: geometry,
    scale: 30, 
    crs: "EPSG:4326",
    bestEffort: true,
    tileScale: 5
}));
print('剔除异常值以后的最大值和最小值:', NDVINoAbnormalValueMin, NDVINoAbnormalValueMax)