顾及物候信息提取冬小麦面积（基于AI earth）

ytkz2022-10-172022-10-18

农业遥感

农业遥感的一个基础方向是农作物分类，农作物分类的基础在于要了解农作物的生长周期，打个比方就是：

1.植物越绿，NDVI值越高

2.植物不绿了，NDVI值变低

3.植物的生长期、成熟期、收割期，即植物是什么时候绿，什么时候黄。

NDVI = （近红外-红）/（近红外+红）

上个式子重点观察分子，也就是近红外-红，用近红外波段和红波段的值来判断植物是否绿油油。

背后的原理基于人眼对颜色的判别。首先说一个更基础的问题：正常的植物为什么是绿色的？那是因为植物会吸收蓝光和红光，对绿光不敏感，所以反射绿光，并且我们眼睛只对可见光敏感，看不见近红外波段之外的光。总之，正常植物大量反射近红外光，小量反射绿光，大量吸收蓝光和红光。

冬小麦的物候及分类策略

冬小麦，和其他植物不一样的地方在于，它是10月播种，这时慢慢变绿但总体不绿；过冬是保持期；生长期是次年的3月到4月，这个阶段慢慢变绿，而且是很绿；收获期是次年的5月到6月，这阶段植被黄了，甚至是秃了。

总结起来就是：

1.如果该像素是冬小麦，那么10月的影像中它的NDVI很低，但不完全低

2.如果该像素是冬小麦，那么3、4月的影像中它的NDVI逐渐变高、保持高值

3.如果该像素是冬小麦，那么5、6月的影像中它的NDVI逐渐变低

4.水体、建筑物的NDVI值全年为负值

简单的数学公式分类策略

把以上4点总结为数学公式，如果该像素为冬小麦，则该像素满足

$\begin{cases} NDVI_{3月}>0.33 \\ NDVI_{10月}<0.5\\ \end{cases}$

根据冬小麦的物候再去区分常绿森林，如果该像素为冬小麦，该像素也满足

$\begin{cases} NDVI_{10月}<0.5 \\ NDVI_{3月}>0.48 \\NDVI_{5月}>-0.12\\ NDVI_{5月}<0.17 \\ \end{cases}$

或许该像素满足

$\begin{cases}NDVI_{10月}<0.5 \\ NDVI_{3月}<0.48\\NDVI_{5月}>-0.12\\NDVI_{5月}<0.5\times NDVI_{3月} \end{cases}$

基于AI earth实现冬小麦提取

# 初始化
import aie

aie.Authenticate()
aie.Initialize()

# 定义矢量范围
feature_collection = aie.FeatureCollection('China_City') \
                        .filter(aie.Filter.eq('city', '许昌市'))
region = feature_collection.geometry()

# 分类规则
# 指定检索数据集，可设置检索的空间和时间范围，以及属性过滤条件（如云量过滤等）
def getl8_ndvi(startdate,enddate):
    dataset = aie.ImageCollection('LANDSAT_LC08_C02_T1_L2') \
                 .filterBounds(region) \
                 .filterDate(startdate, enddate) 
                 # .filter(aie.Filter.lte('eo:cloud_cover', 20.0))
    ndvi =  dataset.map(get_ndvi)  
    return ndvi

def get_ndvi(image):
    ndvi = image.normalizedDifference(['SR_B5', 'SR_B4'])
    return ndvi  

# 黄淮海地区冬小麦典型物候期。播种期 10-11月，旺长期3-4月，成熟期5-6月
NDVI_median = getl8_ndvi('2020-10-11', '2020-11-10').median().clip(region) 
NDVI_max = getl8_ndvi('2021-03-20', '2021-04-20').max().clip(region) 
NDVI_min = getl8_ndvi('2021-05-20', '2021-06-30').min().clip(region) 

# 冬小麦提取规则集
mask1 = NDVI_max.gt(aie.Image.constant(0.33))     # 过滤水体、不透水面、裸地等非植被区域
mask2 = NDVI_median.lt(aie.Image.constant(0.50))  # 过滤森林、春播作物

mask3 = NDVI_max.lte(aie.Image.constant(0.48))    # < 0.48存在
mask4 = NDVI_max.gt(aie.Image.constant(0.48))     # 区分植被区域和非植被区域，> 0.48 为植被区域


mask5 = NDVI_min.gt(aie.Image.constant(-0.12))    # 筛选冬小麦区域，成熟/收获期冬小麦植被指数下降
mask6 = NDVI_min.lt(aie.Image.constant(0.17))   
mask7 = NDVI_min.lt(NDVI_max.add(aie.Image.constant(0.5)))

wheat1 = mask1.And(mask2).And(mask3).And(mask5).And(mask6)
wheat2 = mask1.And(mask2).And(mask4).And(mask5).And(mask7)
wheat  = wheat1.add(wheat2).where((wheat1.add(wheat2)).gt(aie.Image.constant(0)),aie.Image.constant(1))

在AI earth 上可视化

map = aie.Map(
    center=region.getCenter(),
    height=800,
    zoom=7
)

vis_params = {
    'color': '#00FF00'
}

map.addLayer(
    region,
    vis_params,
    'region',
    bounds=region.getBounds()
)

mask_vis  = {
    'min': 0,
    'max': 1,
    'palette': ['#ffffff', '#008000']    # 0:白色, 1:绿色
}

ndvi_vis  = {
    'min': -0.2,
    'max': 0.6,
    'palette': ['#d7191c', '#fdae61', '#ffffc0', '#a6d96a', '#1a9641']
}

map.addLayer(NDVI_median,ndvi_vis, 'NDVI_median', bounds=region.getBounds())
map.addLayer(NDVI_max,ndvi_vis, 'NDVI_max', bounds=region.getBounds())
map.addLayer(NDVI_min,ndvi_vis, 'NDVI_min', bounds=region.getBounds())
map.addLayer(wheat,mask_vis, 'wheat', bounds=region.getBounds())    # 绿色区域为小麦

map