PIE-engine 教程 ——新疆石河子市棉花种植面积提取（阈值法）案例分析

我们之前有有一个案例用来提取大蒜的种植面积，这里我们首先来解释一下阈值法

阈值法，关于阈值法的介绍很多人还不知道，现在让我们一起来看看吧！

1、阈值分割法是一种基于区域的图像分割技术，原理是把图像象素点分为若干类。

2、图像阈值化分割是一种传统的最常用的图像分割方法，因其实现简单、计算量小、性能较稳定而成为图像分割中最基本和应用最广泛的分割技术。

3、它特别适用于目标和背景占据不同灰度级范围的图像。

4、它不仅可以极大的压缩数据量，而且也大大简化了分析和处理步骤，因此在很多情况下，是进行图像分析、特征提取与模式识别之前的必要的图像预处理过程。

5、图像阈值化的目的是要按照灰度级，对像素集合进行一个划分，得到的每个子集形成一个与现实景物相对应的区域，各个区域内部具有一致的属性，而相邻区域不具有这种一致属性。

6、这样的划分可以通过从灰度级出发选取一个或多个阈值来实现。

以上阈值法的介绍来自：阈值法（关于阈值法的介绍）_华夏报道网 (lhyzedu.com)

本文同样是利用NDVI的计算（使用的sentinel2影像用来计算），同时我们利用全球30米分辨率的土地分类数据进行，我们首先看一下使用的结果：

全球地表覆盖数据(GlobeLand30)是中国研制的30米空间分辨率全球地表覆盖数据，2014年已发布2000和2010版，自然资源部于2017年对该数据进行了更新，形成了2020版。
GlobeLand30中包括10个种类，分别为耕地、林地、草地、灌木地、湿地、水体、苔原、人造地表、裸地、冰川和永久积雪。GlobeLand30研制所用的分类影像主要为30米多光谱影像，包括Landsat系列的TM5、ETM+、OLI多光谱影像和中国环境减灾卫星HJ-1多光谱影像，2020版数据还使用了16米的GF1多光谱影像。在影像无云（少云）前提下，优先选择生产基准年或更新年度±2年内植被生长季的多光谱影像，影像获取困难地区则放宽获取时间，从而保证影像的全球覆盖度。

GlobeLand30 V2010数据精度评价由同济大学牵头完成。从全球853幅数据中抽取80个图幅，布设超过15万个检验样本，得出GlobeLand30 V2010数据的总体精度为83.50%，Kappa系数0.78。GlobeLand30 V2020数据精度评价由中国科学院空天信息创新研究院牵头完成。基于景观形状指数抽样模型进行全套数据布点，共布设样本超过23万个。得出GlobeLand30 V2020数据的总体精度为85.72%，Kappa系数0.82。

数据集ID: 

NGCC/GLOBELAND30

时间范围: 2000年-2020年

范围: 全国

来源: NGCC

复制代码段: 

var images = pie.ImageCollection("NGCC/GLOBELAND30")

土地分类的分类状况：

代码：

 
//加载显示石河子市矢量边界数据
var shz = pie.FeatureCollection("NGCC/CHINA_COUNTY_BOUNDARY")
    .filter(pie.Filter.eq("name", "石河子市"))
    .first()
    .geometry();
Map.centerObject(shz, 9);
Map.addLayer(shz, { color: "ff0000ff", fillColor: "00000000", width: 1 }, "石河子市");
 
//加载显示全球地表覆盖GlobeLand30数据集并筛选耕地
var lc = pie.ImageCollection('NGCC/GLOBELAND30')
    .filterDate("2019", "2021")
    .select("B1")
    .first()
    .eq(10)
    .clip(shz);
Map.addLayer(lc, { uniqueValue: ["0", "1"], palette: ["f5fffa", "9acd32"] }, "耕地", false);
 
//加载显示2020年7月石河子市Sentinel-2 L2A合成影像
var img = pie.Image("user/15321576968/SHZ_S2_L2A");
Map.addLayer(img.select(["B4", "B3", "B2"]), { min: 0, max: 3000 }, "石河子市S2_L2A合成影像", false);
 
//计算影像NDVI、NDWI以及SPAD
var green = img.select("B3");
var red = img.select("B4");
var re2 = img.select("B6");
var re3 = img.select("B7");
var nir = img.select("B8");
var NDVI = (nir.subtract(red)).divide(nir.add(red)).rename("NDVI");
var NDWI = (green.subtract(nir)).divide(green.add(nir)).rename("NDWI");
var SPAD = (re3.divide(re2)).multiply(25.34).subtract(28.06).rename("SPAD");
 
//根据敏感特征提取棉花并显示
var cot_RGB = img.select(["B4", "B3", "B2"])
    .updateMask(lc.eq(1).and(nir.gt(4500)).and(NDVI.gt(0.5)));
Map.addLayer(cot_RGB, { min: 0, max: 3000 }, "棉花_RGB")
 
//掩膜获取棉花NDVI
var cot = NDVI.updateMask(lc.eq(1).and(nir.gt(4500)).and(NDVI.gt(0.5)));
 
//按不同阈值显示NDVI
var cot_NDVI = cot.where(cot.gt(0.50).and(cot.lte(0.75)), 1)
    .where(cot.gt(0.75).and(cot.lte(0.80)), 2)
    .where(cot.gt(0.80).and(cot.lte(0.85)), 3)
    .where(cot.gt(0.85).and(cot.lte(0.90)), 4)
    .where(cot.gt(0.90).and(cot.lte(1.00)), 5);
Map.addLayer(cot_NDVI, { min: 1, max: 5, palette: ["E64C00", "E6E600", "55FF00", "33C2FF", "0000FF"] }, "棉花-NDVI");
 
//计算棉花NDVI最大值、最小值以及均值
var NDVI_max = cot.reduceRegion(pie.Reducer.max(), shz, 30).get("NDVI").getInfo();
print("NDVI最大值", NDVI_max);
var NDVI_min = cot.reduceRegion(pie.Reducer.min(), shz, 30).get("NDVI").getInfo();
print("NDVI最小值", NDVI_min);
var NDVI_mean = cot.reduceRegion(pie.Reducer.mean(), shz, 30).get("NDVI").getInfo();
print("NDVI平均值", NDVI_mean);
 
//计算棉花耕地的面积，单位：平方千米
var area = cot.pixelArea()
    .multiply(cot)
    .reduceRegion(pie.Reducer.sum(), shz, 30)
    .getInfo()
    .constant / 1000000;
print("棉花种植面积", area);
 
//掩膜获取棉花NDWI
var wat = NDWI.updateMask(lc.eq(1).and(nir.gt(4500)).and(NDVI.gt(0.5)));
 
//按不同阈值显示NDWI
var water = wat.where(wat.gt(-1.00).and(wat.lte(-0.90)), 1)
    .where(wat.gt(-0.90).and(wat.lte(-0.85)), 2)
    .where(wat.gt(-0.85).and(wat.lte(-0.80)), 3)
    .where(wat.gt(-0.80).and(wat.lte(-0.70)), 4)
    .where(wat.gt(-0.70), 5);
Map.addLayer(water, { min: 1, max: 5, palette: ["FF0000", "FFC800", "B6FF8F", "33C2FF", "0000FF"] }, "棉花-NDWI");
 
//计算棉花NDWI最大值、最小值以及均值
var NDWI_max = wat.reduceRegion(pie.Reducer.max(), shz, 30).get("NDWI").getInfo();
print("NDWI最大值", NDWI_max);
var NDWI_min = wat.reduceRegion(pie.Reducer.min(), shz, 30).get("NDWI").getInfo();
print("NDWI最小值", NDWI_min);
var NDWI_mean = wat.reduceRegion(pie.Reducer.mean(), shz, 30).get("NDWI").getInfo();
print("NDWI平均值", NDWI_mean);
 
//掩膜获取棉花SPAD
var crp = SPAD.updateMask(lc.eq(1).and(nir.gt(4500)).and(NDVI.gt(0.5)));
 
//按不同阈值显示SPAD
var chlorophyll = crp.where(crp.gt(0.0).and(crp.lte(1.5)), 1)
    .where(crp.gt(1.5).and(crp.lte(2.5)), 2)
    .where(crp.gt(2.5).and(crp.lte(3.5)), 3)
    .where(crp.gt(3.5).and(crp.lte(4.5)), 4)
    .where(crp.gt(4.5), 5);
Map.addLayer(crp, { min: 1, max: 5, palette: ["FFFF80", "71EB2F", "55FF00", "216E9E", "0C1078"] }, "叶绿素");
 
//计算棉花SPAD最大值、最小值以及均值
var SPAD_max = crp.reduceRegion(pie.Reducer.max(), shz, 30).get("SPAD").getInfo();
print("叶绿素最大值", SPAD_max);
var SPAD_min = crp.reduceRegion(pie.Reducer.min(), shz, 30).get("SPAD").getInfo();
print("叶绿素最小值", SPAD_min);
var SPAD_mean = crp.reduceRegion(pie.Reducer.mean(), shz, 30).get("SPAD").getInfo();
print("叶绿素平均值", SPAD_mean);

结果：

NDVI最大值

0.9621863183224476

NDVI最小值

0.5003236245954693

NDVI平均值

0.8635523634730466

棉花种植面积

56.16408868367092

NDWI最大值

-0.4891443167305236

NDWI最小值

-0.9020195960807839

NDWI平均值

-0.766581956295227

叶绿素最大值

17.378812949640288

叶绿素最小值

-4.081157757496737

叶绿素平均值

3.5509758941190572

 同样大家可以用画图工具进行自己所选定区域的分析，但是我们得先进行自己进行sentinel2的NDVI计算