前言

農(nóng)業(yè)是國(guó)民經(jīng)濟(jì)的基礎(chǔ)，是人類社會(huì)的衣食之源、生存之本。遙感技術(shù)由于具有探測(cè)范圍廣、信息獲取快、宏觀性強(qiáng)等特點(diǎn)，已廣泛應(yīng)用于農(nóng)作物分類識(shí)別和作物面積估算的研究中。傳統(tǒng)的利用遙感手段進(jìn)行農(nóng)作物分類識(shí)別的方法，大多基于中低分辨率遙感影像（如Landsat、MODIS等）開展，但大量的混合像元會(huì)限制識(shí)別精度的提高。隨著高分辨率的遙感影像不斷免費(fèi)開放，快速精準(zhǔn)地獲取作物信息為政府部門合理制定糧食政策、精準(zhǔn)發(fā)放農(nóng)作物補(bǔ)貼、適時(shí)調(diào)配農(nóng)業(yè)機(jī)械等多方面具有重要意義。

新疆作為中國(guó)最大的棉花種植基地，其棉花產(chǎn)量能夠左右國(guó)際棉價(jià)，在我國(guó)棉花產(chǎn)業(yè)中占有舉足輕重的作用，而棉花作為新疆的主要經(jīng)濟(jì)作物，對(duì)新疆農(nóng)業(yè)經(jīng)濟(jì)和社會(huì)發(fā)展也至關(guān)重要。本應(yīng)用案例主要是基于PIE-Engine平臺(tái)利用Sentinel-2A L2A級(jí)數(shù)據(jù)快速識(shí)別、提取新疆維吾爾族自治區(qū)石河子市的棉花耕地，并反演棉花的NDVI（植被指數(shù)）、NDWI（水體指數(shù)）以及SPAD（葉綠素），以期獲取高精度、高分辨率的區(qū)域棉花耕地分布及其生長(zhǎng)情況，服務(wù)一帶一路核心區(qū)的農(nóng)業(yè)信息化建設(shè)。

PIE-Engine平臺(tái)介紹

航天宏圖致力于加速我國(guó)遙感技術(shù)的發(fā)展進(jìn)程，依托行業(yè)多年技術(shù)積累，獨(dú)立自主研發(fā)了安全可控的開放式遙感云計(jì)算平臺(tái)：PIE-Engine（Pixel Information Expert Engine，像素專家引擎），實(shí)現(xiàn)了遙感數(shù)據(jù)按需獲取、運(yùn)算以及專題信息聚焦服務(wù)，以滿足對(duì)地觀測(cè)數(shù)據(jù)獲取能力飛速增長(zhǎng)帶來的信息高效化處理和服務(wù)需求。目前平臺(tái)數(shù)據(jù)總量已經(jīng)超過1.5PB，存儲(chǔ)了國(guó)內(nèi)與國(guó)外的近80種遙感數(shù)據(jù)集，超過250萬景影像數(shù)據(jù)，涵蓋了光學(xué)、微波、高光譜、高程、人口、氣象、夜光等多種數(shù)據(jù)集，國(guó)內(nèi)數(shù)據(jù)包括高分、風(fēng)云、海洋等系列，國(guó)外數(shù)據(jù)包括Landsat、MODIS、Sentinel以及Himawari等。

棉花提取方法

首先，利用GlobaLand 30全球地表覆蓋分類數(shù)據(jù)獲取石河子市的所有耕地區(qū)域。另外，通過分析石河子市棉花耕地在影像中的敏感特征（光譜與指數(shù)），棉花在紅光波段的反射率相比其它波段表現(xiàn)出了更高的敏感性，并且該區(qū)域處于生長(zhǎng)季的棉花NDVI指數(shù)普遍大于0.5。由此，我們可以根據(jù)這兩個(gè)敏感特征在影像中的耕地區(qū)域內(nèi)快速識(shí)別與提取棉花耕地，并計(jì)算棉花的NDVI、NDWI以及葉綠素。計(jì)算公式如下：

式中：B3為影像的綠光波段反射率，B4為影像的紅光波段反射率，B6為紅邊2波段反射率，B7為紅邊3波段反射率，B8為近紅外波段的反射率。

提取結(jié)果

示例代碼

代碼鏈接：

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

向下滑動(dòng)閱覽

1.//加載顯示石河子市矢量邊界數(shù)據(jù)??

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數(shù)據(jù)集并篩選耕地??

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.//計(jì)算影像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.//根據(jù)敏感特征提取棉花并顯示??

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.//計(jì)算棉花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.//計(jì)算棉花耕地的面積，單位：公頃??

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.//計(jì)算棉花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.//計(jì)算棉花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);?

小編為各位正在關(guān)注和學(xué)習(xí)PIE-Engine Studio的用戶強(qiáng)烈安利一期教學(xué)視頻，感興趣的小伙伴，快快點(diǎn)擊查看呀！今后，我們會(huì)在相關(guān)的公眾號(hào)文章下面，附送一期視頻教程，請(qǐng)大家持續(xù)關(guān)注哦~