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

本文链接：https://blog.csdn.net/qq_31988139/article/details/127029537

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

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

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")

名称	类型	无效值	空间分辨率	时间分辨率	描述信息
B1	Byte	0	30m	年	时间范围不连续，分别为2000年，2010年，2020年共三年的地表覆盖产品

土地分类的分类状况：

名称	内容	代码
耕地	用于种植农作物的土地，包括水田、灌溉旱地、雨养旱地、菜地、牧草种植地、大棚用地、以种植农作物为主间有果树及其他经济乔木的土地，以及茶园、咖啡园等灌木类经济作物种植地。	10
林地	乔木覆盖且树冠盖度超过30%的土地，包括落叶阔叶林、常绿阔叶林、落叶针叶林、常绿针叶林、混交林，以及树冠盖度为10-30%的疏林地。	20
草地	天然草本植被覆盖，且盖度大于10%的土地，包括草原、草甸、稀树草原、荒漠草原，以及城市人工草地等。	30
灌木地	灌木覆盖且灌丛覆盖度高于30%的土地，包括山地灌丛、落叶和常绿灌丛，以及荒漠地区覆盖度高于10%的荒漠灌丛。	40
湿地	位于陆地和水域的交界带，有浅层积水或土壤过湿的土地，多生长有沼生或湿生植物。包括内陆沼泽、湖泊沼泽、河流洪泛湿地、森林/灌木湿地、泥炭沼泽、红树林、盐沼等。	50
水体	陆地范围液态水覆盖的区域，包括江河、湖泊、水库、坑塘等。	60
苔原	寒带及高山环境下由地衣、苔藓、多年生耐寒草本和灌木植被覆盖的土地，包括灌丛苔原、禾本苔原、湿苔原、高寒苔原、裸地苔原等。	70
人造地表	由人工建造活动形成的地表，包括城镇等各类居民地、工矿、交通设施等，不包括建设用地内部连片绿地和水体。	80
裸地	植被覆盖度低于10%的自然覆盖土地，包括荒漠、沙地、砾石地、裸岩、盐碱地等。	90
冰川和永久积雪	由永久积雪、冰川和冰盖覆盖的土地，包括高山地区永久积雪、冰川，以及极地冰盖等。	100

date

string

影像时间

代码：


//加载显示石河子市矢量边界数据
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最小值