PIE-Engine实践篇｜巧用map算子提高效率

　　PIE-Engine 遥感云计算平台上集成了海量遥感数据、强大的算力和丰富的算子，使得大规模尺度长时间序列的遥感数据分析可以在线快速完成，其中一大利器便是map算子。本期我们来详细介绍map算子的作用，它和传统for循环的区别，以及如何巧用map算子优化代码，提高计算效率。

　　一、 map算子的作用

　　map算子是我们在云计算平台上处理影像时常用的算子，可以对集合(List、FeatureCollection、ImageCollection)中的每一个元素进行操作，操作完成返回的结果仍是一个集合对象。以FeatureCollection调用map算子为例，其基本形式如下：

　　var featureCol = pie.FeatureCollection('NGCC/CHINA_PROVINCE_BOUNDARY');

　　var featureColNew = featureCol.map(function (feature) {

　　var geometry = feature.geometry();

　　var featureNew = pie.Feature(geometry.centroid());

　　return featureNew;

　　})

　　上述代码取出featureCol 中的每一个feature，然后求取各feature的几何中心，得到一个新的矢量集合-featureColNew。

　　二、map算子和for循环的区别

　　当然，有些map操作用传统的for/while循环也可以实现，但二者在计算机制上有本质区别，for循环是在前端控制变量的循环，且是顺序执行代码，效率比较慢，一般在PIE-Engine里很少用。而map则是在后台服务器上直接作用于集合中的元素进行循环操作，执行效率很快，相当于“并行执行任务”， 但在map里一般不再用 print、Map、Export 以及Reducer等方法。下面分别对比在集合中使用map算子和for循环的不同之处。

　　针对List中的每个元素进行循环计算

　　//对List中的每个元素加1，采用map算子

　　var list1 = pie.List([1,2,3,2]);

　　var list2 = list1.map(function (value) {

　　return pie.Number(value).add(1);

　　});

　　print("list2",list2);

　　输出结果：

　　//对List中的每个元素加1，采用for循环

　　var list1 = pie.List([1,2,3,2]);

　　var list2 = [];

　　for(var i = 0;i

　　list2.push(pie.Number(list1.get(i)).add(1).getInfo());

　　}

　　print("list2",list2);

　　输出结果：

　　List集合调用map算子后便可直接作用于其中的元素，而采用for循环则还需要增加诸多限制，才能实现相同的功能。

　　针对FeatureCollection中的每个Feature进行循环计算

　　利用中国省级行政区划矢量数据计算每个省的面积，分别采用map算子和for循环进行计算。

　　(1)取FeatureCollection中的每个元素并计算其面积，采用map算子：

　　向上滑动阅览

　　var featureCol = pie.FeatureCollection().load('NGCC/CHINA_PROVINCE_BOUNDARY');

　　var featureColNew = featureCol.map(function (feature) {

　　var name = feature.get("name");

　　var geometry = feature.geometry();

　　var area = geometry.area().divide(1000000);

　　feature = feature.set("name",name);

　　feature = feature.set("area",area);

　　return feature;

　　})

　　var result = featureColNew.reduceColumns(pie.Reducer.toList(2),["name","area"]);

　　print("result",result);

　　Map.addLayer(featureCol,,"featureCol");

　　Map.setCenter(118,39.7,3);

　　代码链接：

　　https://engine.piesat.cn/engine-share/shareCode.html?id=cbe5f4f5efca4b5faa8344b0ab325c04

　　采用map算子获取中国各省的面积

　　(2)取FeatureCollection中的每个元素并计算其面积，采用for循环：

　　向上滑动阅览

　　var featureCol = pie.FeatureCollection().load('NGCC/CHINA_PROVINCE_BOUNDARY');

　　print(featureCol)

　　for(i = 0; i

　　var name = featureCol.getAt(i).get("name");

　　var geometry = featureCol.getAt(i).geometry();

　　var area = geometry.area().divide(1000000);

　　print(name,area);

　　}

　　Map.addLayer(featureCol,,"featureCollection");

　　Map.setCenter(118,39.7,3);

　　代码链接：

　　https://engine.piesat.cn/engine-share/shareCode.html?id=4198d136975841269e367aec1ae47e70

　　运行结果：

　　▲for循环输出各省面积

　　上述代码采用map算子，2秒即可出结果，而for循环则需要40秒，等待时间较长，从而体现map算子的优势。

　　针对ImageCollection中的每个Image进行循环计算

　　(1)取ImageCollection中的每个元素计算NDVI，采用map算子：

　　向上滑动阅览

　　var featureCol = pie.FeatureCollection().load('NGCC/CHINA_PROVINCE_BOUNDARY');

　　var beijing = featureCol.filter(pie.Filter.eq("name", "北京市"));

　　var bjGeo = beijing.getAt(0).geometry();

　　var visParams = ;

　　Map.addLayer(beijing, visParams, "北京市");

　　Map.centerObject(beijing, 6);

　　var imageCol = pie.ImageCollection("LC08/01/T1")

　　.filterDate("2019-12-01", "2019-12-31")

　　.filterBounds(bjGeo);

　　print("imageCol", imageCol);

　　var imageColNDVI = imageCol.map(function (image) {

　　// NDVI计算

　　var img_Nir = image.select("B5");

　　var img_Red = image.select("B4");

　　var img_NDVI = img_Nir.subtract(img_Red).divide(img_Nir.add(img_Red)).rename("NDVI");

　　return image.addBands(img_NDVI);

　　});

　　print("imageColNDVI", imageColNDVI);

　　//NDVI绘制样式

　　var visParamNDVI = {

　　min: -0.2,

　　max: 0.8,

　　palette: 'CA7A41, CE7E45, DF923D, F1B555, FCD163, 99B718, '+

　　'74A901, 66A000, 529400,3E8601, 207401, 056201, 004C00,'+

　　'023B01, 012E01, 011D01, 011301'

　　};

　　var imageNDVI = imageColNDVI.select("NDVI").mosaic().clip(bjGeo);

　　Map.addLayer(imageNDVI, visParamNDVI, "Layer_NDVI");

　　代码链接：

　　https://engine.piesat.cn/engine-share/shareCode.html?id=4ec895ce6e3b411daab877713ca81f45

　　输出结果：

　　(2)取ImageCollection中的每个元素计算NDVI，采用for循环：

　　向上滑动阅览

　　var featureCol = pie.FeatureCollection().load('NGCC/CHINA_PROVINCE_BOUNDARY');

　　var beijing = featureCol.filter(pie.Filter.eq("name", "北京市"));

　　var bjGeo = beijing.getAt(0).geometry();

　　var visParams = { color: "ff0000ff", fillColor: "00000000" };

　　Map.addLayer(beijing, visParams, "北京市");

　　Map.centerObject(beijing, 6);

　　var imageCol = pie.ImageCollection("LC08/01/T1")

　　.filterDate("2019-8-01", "2019-8-30")

　　.filterBounds(bjGeo);

　　print("imageCol", imageCol);

　　var newCol=[];

　　for (i = 0; i

　　var image = imageCol.getAt(i);

　　var img_Nir = image.select("B5");

　　var img_Red = image.select("B4");

　　var img_NDVI = img_Nir.subtract(img_Red).divide(img_Nir.add(img_Red))

　　.rename("NDVI");

　　image = image.addBands(img_NDVI);

　　newCol.push(image);

　　}

　　var newImageCol=pie.ImageCollection().fromImages(newCol);

　　print("newImageCol", newImageCol);

　　//NDVI绘制样式

　　var visParamNDVI = {

　　min: -0.2,

　　max: 0.8,

　　palette: 'CA7A41, CE7E45, DF923D, F1B555, FCD163, 99B718, ' +

　　'74A901, 66A000, 529400,3E8601, 207401, 056201, 004C00,' +

　　'023B01, 012E01, 011D01, 011301'

　　};

　　var imageNDVI = newImageCol.select("NDVI").mosaic().clip(bjGeo);

　　Map.addLayer(imageNDVI, visParamNDVI, "Layer_NDVI");

　　代码链接：

　　https://engine.piesat.cn/engine-share/shareCode.html?id=ef81bdf97c8d45809c3101f5f8dee4c4

　　上述代码采用map算子和for循环的计算时间基本相同， 但在使用for循环时要重新构造一个存放新的image集合的数组，并在循环结束后重新构造一个由这些image组成的ImageCollection，在步骤上较为复杂。

　　因此，在遥感云计算平台中，涉及到对集合中的每个元素进行操作时，推荐使用map算子，快速且便利。但前述提到，在map算子中不宜再用 print、Map、Export 以及Reducer等方法，如果需要对集合中的每个元素进行统计操作，则可用for循环代替，以下代码展示了使用for循环统计北京各城区提取的建筑区面积：

　　向上滑动阅览

　　var chn = pie.FeatureCollection('NGCC/CHINA_PROVINCE_BOUNDARY');

　　var beijing = chn.filter(pie.Filter.eq("name", "北京市"));

　　var bjGeo = beijing.getAt(0).geometry();

　　var chn2 = pie.FeatureCollection('NGCC/CHINA_COUNTY_BOUNDARY');

　　var beijing2 = chn2.filter(pie.Filter.eq("cname", "北京市"));

　　var bjGeo2 = beijing2.getAt(0).geometry();

　　Map.centerObject(beijing, 7)

　　var visParams = { color: "ff0000ff", fillColor: "00000000" };

　　Map.addLayer(beijing, visParams, "北京市边界");

　　Map.addLayer(beijing2, visParams, "北京市各区");

　　function rmCloud(image) {

　　var qa = image.select("BQA");

　　var cloudMask = qa.bitwiseAnd(1

　　return image.updateMask(cloudMask);

　　}

　　var l8 = pie.ImageCollection("LC08/01/T1")

　　.filterDate("2020-7-1", "2020-8-30")

　　.filterBounds(bjGeo)

　　.select(["B2", "B3", "B4", "B5", "B6", "BQA"])

　　.map(rmCloud);

　　//计算用到的指数

　　function imgCalculate(image) {

　　var green = image.select("B3");

　　var red = image.select("B4");

　　var nir = image.select("B5");

　　var swir1 = image.select("B6");

　　var ndvi = (nir.subtract(red)).divide(nir.add(red)).rename("NDVI");

　　var mndwi = green.subtract(swir1)

　　.divide(green.add(swir1))

　　.rename("MNDWI");

　　return image.addBands(ndvi).addBands(mndwi);

　　}

　　var imgCol = l8.map(imgCalculate)

　　.mean()

　　.clip(bjGeo);

　　Map.addLayer(imgCol, { min: 0, max: 3000, bands: ["B4", "B3", "B2"] }, "imgCol");

　　var ndvi = imgCol.select("NDVI");

　　var nonVeg = imgCol.updateMask(ndvi.lte(0.6));

　　Map.addLayer(nonVeg, { min: 0, max: 3000, bands: ["B4", "B3", "B2"] }, "nonVeg");

　　var mndwi = imgCol.select("MNDWI");

　　var nonVegWater = nonVeg.updateMask(mndwi.lte(0.3));

　　Map.addLayer(nonVegWater, { min: 0, max: 3000, bands: ["B4", "B3", "B2"] }, "nonVegWater");

　　var nonVegWater2 = nonVegWater.select("B2").gt(0);

　　Map.addLayer(nonVegWater2, { min: 0, max: 1, palette: "ff0000" }, "nonVegWater2");

　　var areaImage = nonVegWater2.pixelArea().multiply(nonVegWater2.gt(0));

　　var area = areaImage.reduceRegion(pie.Reducer.sum(), bjGeo, 30);

　　print("北京市建筑用地面积(单位：平方千米): ", pie.Number(area.get("constant")).divide(1000000));

　　var newFeatures = [];

　　for (var i = 0; i

　　var temp = nonVegWater2.clip(beijing2.getAt(i).geometry());

　　var areaImg = temp.pixelArea().multiply(temp.gt(0));

　　var area = areaImg.reduceRegion(pie.Reducer.sum(), beijing2.getAt(i).geometry(), 30);

　　newFeatures.push(beijing2.getAt(i).set("area", area.get("constant")));

　　}

　　beijing2 = pie.FeatureCollection(newFeatures);

　　print("beijing2",beijing2);

　　var result = beijing2.reduceColumns(pie.Reducer.toList(2), ["name", "area"]);

　　print("result", result);

　　代码链接：

　　https://engine.piesat.cn/engine-share/shareCode.html?id=bf885222bf584f33aaa96c2a0093c03d

　　输出结果：

　　三、巧用map算子，提高计算效率

　　map算子带来便利的同时，也将大量的计算转移给了后台计算服务。通常我们在影像处理时会计算多个指数，单个指数的计算公式通常会写成函数的形式方便调用。多个map算子共同使用时，相当于多次循环，所涉及到的运算量也很大，怎样使用map算子才能优化计算逻辑，提高计算效率呢?

　　下面以去云后计算裸土指数BSI、改进的归一化差异水体指数MNDWI、增强型的裸土指数EBSI、建筑用地指数NDBI等4个指数的计算为例，来说明如何在代码级别进行优化来减少计算量从而提高计算效率。

　　方法1，将每个指数的计算单独写成函数，在筛选完影像后依次用map算子调用：

　　向上滑动阅览

　　var chn = pie.FeatureCollection('NGCC/CHINA_PROVINCE_BOUNDARY');

　　var beijing = chn.filter(pie.Filter.eq("name", "北京市"));

　　var bjGeo = beijing.getAt(0).geometry();

　　Map.centerObject(beijing, 7)

　　var visParams = { color: "ff0000ff", fillColor: "00000000" };

　　Map.addLayer(beijing, visParams, "北京市边界");

　　//去云

　　function rmCloud(image) {

　　var qa = image.select("BQA");

　　var cloudMask = qa.bitwiseAnd(1

　　return image.updateMask(cloudMask);

　　}

　　//裸土指数BSI: [(B06 + B04)-(B05 + B02)]/[(B06 + B04)+(B05 + B02)]

　　function BSI(image) {

　　var nir = image.select("B5");

　　var swir1 = image.select("B6");

　　var red = image.select("B4");

　　var blue = image.select("B2");

　　var bsi = (swir1.add(red).subtract(nir.add(blue)))

　　.divide(swir1.add(red).add(nir.add(blue)));

　　return image.addBands(bsi.rename("BSI"));

　　}

　　//改进的归一化差异水体指数MNDWI: (B03 - B06)/(B03 + B06)

　　function MNDWI(image) {

　　var mndwi = image.select("B3").subtract(image.select("B6"))

　　.divide(image.select("B3").add(image.select("B6")))

　　return image.addBands(mndwi.rename("MNDWI"));

　　}

　　//增强型的裸土指数EBSI：(BSI - MNDWI)/(BSI + MNDWI)

　　function EBSI(image) {

　　var mndwi = image.select("MNDWI");

　　var bsi = image.select("BSI");

　　var ebsi = (bsi.subtract(mndwi)).divide(bsi.add(mndwi));

　　return image.addBands(ebsi.rename("EBSI"));

　　}

　　//建筑用地指数NDBI: (B06 - B05)/(B06 + B05)

　　function NDBI(image) {

　　var ndbi = image.select("B6").subtract(image.select("B5"))

　　.divide(image.select("B6").add(image.select("B5")))

　　return image.addBands(ndbi.rename("NDBI"));

　　}

　　//加载Landsat 8影像数据集合

　　var imgColl8 = pie.ImageCollection("LC08/01/T1")

　　.filterDate("2020-4-1", "2020-4-30")

　　.filterBounds(bjGeo)

　　.select(["B2", "B3", "B4", "B5", "B6", "B7", "BQA"])

　　.map(rmCloud)

　　.map(BSI)

　　.map(MNDWI)

　　.map(EBSI)

　　.map(NDBI)

　　.mean()

　　.clip(bjGeo);

　　print("imgColl8",imgColl8);

　　Map.addLayer(imgColl8, { min: 0, max: 3000, bands: ["B4", "B3", "B2"] }, "imgCol1");

　　Map.addLayer(imgColl8.select("BSI"), { min: -0.2, max: 0.2 }, "BSI");

　　Map.addLayer(imgColl8.select("NDBI"), { min: -0.2, max: 0.2 }, "NDBI");

　　Map.addLayer(imgColl8.select("MNDWI"), , "MNDWI");

　　Map.addLayer(imgColl8.select("EBSI"), { min: 0, max: 1 }, "EBSI");

　　代码链接：

　　https://engine.piesat.cn/engine-share/shareCode.html?id=402cadc45985400c8650620b9d93bad9

　　上述代码共调用了5次map算子，相当于对筛选后的影像进行了5次遍历，除去云操作外，每一次都是对每景影像进行波段运算完后添加计算后的指数，将大量的算力消耗在了重复读取和写入上，输出计算结果用时85秒，完整显示用时135秒。

　　依次调用多个map计算每个指数

　　方法2，将去云及指数运算封装成一个函数，在筛选完影像后仅用map算子调用一次：

　　向上滑动阅览

　　var chn = pie.FeatureCollection('NGCC/CHINA_PROVINCE_BOUNDARY');

　　var beijing = chn.filter(pie.Filter.eq("name", "北京市"));

　　var bjGeo = beijing.getAt(0).geometry();

　　Map.centerObject(beijing, 7)

　　var visParams = { color: "ff0000ff", fillColor: "00000000" };

　　Map.addLayer(beijing, visParams, "北京市边界");

　　var imgColl8 = pie.ImageCollection("LC08/01/T1")

　　.filterDate("2020-4-1", "2020-4-30")

　　.filterBounds(bjGeo)

　　.select(["B2", "B3", "B4", "B5", "B6", "B7", "BQA"]);

　　//影像处理函数

　　function imgCalculate(image) {

　　var qa = image.select("BQA");

　　var cloudMask = qa.bitwiseAnd(1

　　image = image.updateMask(cloudMask);

　　var blue = image.select("B2");

　　var green = image.select("B3");

　　var red = image.select("B4");

　　var nir = image.select("B5");

　　var swir1 = image.select("B6");

　　var bsi = (swir1.add(red).subtract(nir.add(blue)))

　　.divide(swir1.add(red).add(nir.add(blue)))

　　.rename("BSI");

　　var mndwi = green.subtract(swir1)

　　.divide(green.add(swir1))

　　.rename("MNDWI");

　　var ebsi = (bsi.subtract(mndwi)).divide(bsi.add(mndwi))

　　.rename("EBSI");

　　var ndbi = swir1.subtract(nir)

　　.divide(swir1.add(nir))

　　.rename("NDBI");

　　return image.addBands(mndwi).addBands(bsi).addBands(ebsi)

　　.addBands(ndbi);

　　}

　　imgColl8 = imgColl8.map(imgCalculate).mean().clip(bjGeo);

　　print("imgColl8", imgColl8);

　　Map.addLayer(imgColl8, { min: 0, max: 3000, bands: ["B4", "B3", "B2"] }, "imgCol1");

　　Map.addLayer(imgColl8.select("BSI"), { min: -0.2, max: 0.2 }, "BSI");

　　Map.addLayer(imgColl8.select("NDBI"), { min: -0.2, max: 0.2 }, "NDBI");

　　Map.addLayer(imgColl8.select("MNDWI"), , "MNDWI");

　　Map.addLayer(imgColl8.select("EBSI"), { min: 0, max: 1 }, "EBSI");

　　代码链接：

　　https://engine.piesat.cn/engine-share/shareCode.html?id=ce08fdb1cca147acb6ae6da7629fb49f

　　将指数计算封装到一个函数中调用

　　上述代码仅调用了1次map算子，对筛选后的影像进行了去云操作以及4个指数波段的添加，仅需对每景影像读取写入一次，大大节省了算力，只用8秒即可输出结果并完整显示，相比于方法1，计算速度提升近10倍。

　　学会了这种思路后，平常我们在编写代码时，注意将计算过程优化，就可以提高计算效率，减少等待时间。