基于PIE-Engine的新疆地区棉花种植面积提取
基于PIE-Engine的新疆地区棉花种植面积提取
前言
农业是国民经济的基础，是人类社会的衣食之源、生存之本。遥感技术由于具有探测范围广、信息获取快、宏观性强等特点，已广泛应用于农作物分类识别和作物面积估算的研究中。传统的利用遥感手段进行农作物分类识别的方法，大多基于中低分辨率遥感影像（如Landsat、MODIS等）开展，但大量的混合像元会限制识别精度的提高。随着高分辨率的遥感影像不断免费开放，快速精准地获取作物信息为政府部门合理制定粮食政策、精准发放农作物补贴、适时调配农业机械等多方面具有重要意义。
新疆作为中国最大的棉花种植基地，其棉花产量能够左右国际棉价，在我国棉花产业中占有举足轻重的作用，而棉花作为新疆的主要经济作物，对新疆农业经济和社会发展也至关重要。本应用案例主要是基于PIE-Engine平台利用Sentinel-2A L2A级数据快速识别、提取新疆维吾尔自治区石河子市的棉花耕地，并反演棉花的NDVI（植被指数）、NDWI（水体指数）以及SPAD（叶绿素），以期获取高精度、高分辨率的区域棉花耕地分布及其生长情况，服务一带一路核心区的农业信息化建设。
石河子市概况图
PIE-Engine平台介绍
航天宏图致力于加速我国遥感技术的发展进程，依托行业多年技术积累，独立自主研发了安全可控的开放式遥感云计算平台：PIE-Engine（Pixel Information Expert Engine，像素专家引擎），实现了遥感数据按需获取、运算以及专题信息聚焦服务，以满足对地观测数据获取能力飞速增长带来的信息高效化处理和服务需求。目前平台数据总量已经超过1.5PB，存储了国内与国外的近80种遥感数据集，超过250万景影像数据，涵盖了光学、微波、高光谱、高程、人口、气象、夜光等多种数据集，国内数据包括高分、风云、海洋等系列，国外数据包括Landsat、MODIS、Sentinel以及Himawari等。
PIE-Engine Studio代码在线编辑器
棉花提取方法
首先，利用GlobaLand 30全球地表覆盖分类数据获取石河子市的所有耕地区域。另外，通过分析石河子市棉花耕地在影像中的敏感特征（光谱与指数），棉花在红光波段的反射率相比其它波段表现出了更高的敏感性，并且该区域处于生长季的棉花NDVI指数普遍大于0.5。由此，我们可以根据这两个敏感特征在影像中的耕地区域内快速识别与提取棉花耕地，并计算棉花的NDVI、NDWI以及叶绿素。计算公式如下：
式中：B3为影像的绿光波段反射率，B4为影像的红光波段反射率，B6为红边2波段反射率，B7为红边3波段反射率，B8为近红外波段的反射率。
提取结果
棉花提取结果
棉花NDVI
棉花NDWI
棉花叶绿素
示例代码
代码链接：
https://engine.piesat.cn/engine-share/shareCode.html?id=4481acfe5c3d4ee79f2b739f056bc730
代码运行结果
向下滑动阅览
1.//加载显示石河子市矢量边界数据
2.var shz = pie.FeatureCollection("NGCC/CHINA_COUNTY_BOUNDARY")
3. .filter(pie.Filter.eq("name", "石河子市"))
4. .first()
5. .geometry();
6.Map.centerObject(shz, 9);
7.Map.addLayer(shz, {color: "ff0000ff", fillColor: "00000000", width: 1}, "石河子市");
8.
9.//加载显示全球地表覆盖GlobeLand30数据集并筛选耕地
10.var lc = pie.ImageCollection('NGCC/GLOBELAND30')
11. .filterDate("2019", "2021")
12. .select("B1
13. .first()
14. .eq(10)
15. .clip(shz);
16.Map.addLayer(lc, {uniqueValue: ["0", "1"], palette: ["f5fffa", "9acd32"]}, "耕地", false);
17.
18.//加载显示2020年7月石河子市Sentinel-2 L2A合成影像
19.var img = pie.Image("user/3408/SHZ_S2_L2A");
20.Map.addLayer(img.select(["B4", "B3", "B2"]), {min:0, max: 3000}, "石河子市S2_L2A合成影像",false);
21.
22.//计算影像NDVI、NDWI以及SPAD
23.var green = img.select("B3");
24.var red = img.select("B4");
25.var re2 = img.s
26.var re3 = img.select("B7");
27.var nir = img.select("B8");
28.var NDVI = (nir.subtract(red)).divide(nir.add(red)).rename("NDVI");
29.var NDWI = (green.subtract(nir)).divide(green.add(nir)).rename("NDWI");
30.var SPAD = (re3.divide(re2)).multiply(25.34).subtract(28.06).rename("SPAD");
31.
32.//根据敏感特征提取棉花并显示
33.var cot_RGB = img.select(["B4", "B3", "B2"])
34. .updateMask(lc.eq(1).and(nir.gt(4500)).and(NDVI.gt(0.5)));
35.Map.addLayer(cot_RGB,
36.
37.//掩膜获取棉花NDVI
38.var cot = NDVI.updateMask(lc.eq(1).and(nir.gt(4500)).and(NDVI.gt(0.5)));
39.
40.//按不同阈值显示NDVI
41.var cot_NDVI = cot.where(cot.gt(0.50).and(cot.lte(0.75)), 1)
42. .where(cot.gt(0.75).and(cot.lte(0.80)), 2)
43. .where(cot.gt(0.80).and(cot.lte(0.85)), 3)
44. .where(cot.gt(0.85).and(cot.lte(0.90)), 4)
45. .where(cot.gt(0.90).and(cot.lte(1.00)), 5);
46.Map.addLayer(cot_N
47.
48.//计算棉花NDVI最大值、最小值以及均值
49.var NDVI_max = cot.reduceRegion(pie.Reducer.max(), shz, 30).get("NDVI").getInfo();
50.print("NDVI最大值", NDVI_max);
51.var NDVI_min = cot.reduceRegion(pie.Reducer.min(), shz, 30).get("NDVI").getInfo();
52.print("NDVI最小值", NDVI_min);
53.var NDVI_mean = cot.reduceRegion(pie.Reducer.mean(), shz, 30).get("NDVI").getInfo();
54.print("NDVI平均值", NDVI_mean);
55.
56.//计算棉花耕地的面积，单位：公顷
57.var area = cot.pixelArea()
58.
59. .reduceRegion(pie.Reducer.sum(), shz, 30)
60. .getInfo()
61. .constant/1000000;
62.print("棉花种植面积", area);
63.
64.//掩膜获取棉花NDWI
65.var wat = NDWI.updateMask(lc.eq(1).and(nir.gt(4500)).and(NDVI.gt(0.5)));
66.
67.//按不同阈值显示NDWI
68.var water = wat.where(wat.gt(-1.00).and(wat.lte(-0.90)), 1)
69. .where(wat.gt(-0.90).and(wat.lte(-0.85)), 2)
70. .where(wat.gt(-0.85).and(wat.lte(-0.80)), 3)
71. .where(wat.gt(-0.80).and(wat.lte(-0.70)), 4)
72. .where(wat.gt(-0.70), 5);
73.Map.addLayer(water, {min: 1, max: 5, palette: ["FF0000", "FFC800", "B6FF8F", "33C2FF", "0000FF"]}, "棉花-NDWI");
74.
75.//计算棉花NDWI最大值、最小值以及均值
76.var NDWI_max = wat.reduceRegion(pie.Reducer.max(), shz, 30).get("NDWI").getInfo();
77.print("NDWI最大值", NDWI_max);
78.var NDWI_min = wat.reduceRegion(pie.Reducer.min(), shz, 30).get("NDWI").getInfo();
79.print("NDWI
80.var NDWI_mean = wat.reduceRegion(pie.Reducer.mean(), shz, 30).get("NDWI").getInfo();
81.print("NDWI平均值", NDWI_mean);
82.
83.//掩膜获取棉花SPAD
84.var crp = SPAD.updateMask(lc.eq(1).and(nir.gt(4500)).and(NDVI.gt(0.5)));
85.
86.//按不同阈值显示SPAD
87.var chlorophyll = crp.where(crp.gt(0.0).and(crp.lte(1.5)), 1)
88. .where(crp.gt(1.5).and(crp.lte(2.5)), 2)
89. .where(crp.gt(2.5).and(crp.lte(3.5)), 3)
90. .wh
91. .where(crp.gt(4.5), 5);
92.Map.addLayer(crp, {min:1, max:5, palette: ["FFFF80", "71EB2F", "55FF00", "216E9E", "0C1078"]}, "叶绿素");
93.
94.//计算棉花SPAD最大值、最小值以及均值
95.var SPAD_max = crp.reduceRegion(pie.Reducer.max(), shz, 30).get("SPAD").getInfo();
96.print("叶绿素最大值", SPAD_max);
97.var SPAD_min = crp.reduceRegion(pie.Reducer.min(), shz, 30).get("SPAD").getInfo();
98.print("叶绿素最小值", SPAD_min);
99.var SPAD_mean = crp.reduceRegion(pie.Reducer.mean
100.print("叶绿素平均值", SPAD_mean);
小编为各位正在关注和学习PIE-Engine Studio的用户强烈安利一期教学视频，感兴趣的小伙伴，快快点击查看呀！今后，我们会在相关的公众号文章下面，附送一期视频教程，请大家持续关注哦~
