基于R语言的NDVI的Sen-MK趋势检验

绝地灬酷狼 2024-04-17 06:17 41阅读 0赞

本实验拟分析艾比湖地区2010年至2020年间的NDVI数据，数据从MODIS遥感影像中提取的NDVI值，在GEE遥感云平台上将影像数据下载下来。代码如下：

import ee
    import geemap 
    geemap.set_proxy(port=7890)# 设置全局网络代理
    Map = geemap.Map()
    
    # 指定艾比湖地区数据范围
    region = ee.Geometry.BBox(82.433,
                              44.367,
                              84.5,
                              45.267)
                              
    def get_mean_ndvi(year):
      
      y1,y2 = f'{
        year}-01-01', f'{
        year+1}-01-01'
      
      ndvi_collection = (ee.ImageCollection('MODIS/MCD43A4_006_NDVI') 
                         .filterDate(y1,y2) 
                         .filterBounds(region))
      ndvi = ndvi_collection.reduce(ee.Reducer.mean())
      
      geemap.ee_export_image_to_drive(
        ndvi, description=f'ndvi{
        year}', folder='image',scale=1000,region=region)
      
    for y in range(2010,2021):
      get_mean_ndvi(y)

影像会下载到Google云盘中，通过手动下载到本地，其根目录结构如下：  
![在这里插入图片描述][0fdaefaf6d7f4761b5868abeec8ef557.png]

图1 根目录结构

下载该10年间的数据后，打开RStdio并导入将趋势检验中将使用的R包。代码如下：

library(sp)
    library(raster)
    library(rgdal)
    library(trend)
    
    setwd('E:/CN/NDVI')
    fl <- list.files(pattern = '*tif$')
    firs <- raster(fl[1])
    
    for (i in 1:10) {
        
        r <- raster(fl[i])
        firs <- stack(firs, r)
    }
    
    fun <- function(y){
        
        if(length(na.omit(y)) <10) return(c(NA, NA, NA))   #删除数据不连续含有NA的像元
        av <- mean(y,na.rm=T)
        MK_estimate <- sens.slope(ts(na.omit(y), start = 2010, end = 2020, frequency = 1), conf.level = 0.95) #Sen斜率估计
        slope <- MK_estimate$estimate
        MK_test <- MK_estimate$p.value
    #    Zs <- MK_estimate$statistic
        return(c(av, slope, MK_test))
    }
    
    e <- calc(firs, fun)   #栅格计算
    #e_Zs <- subset(e,1)  #提取Z统计量
    e_mean <- subset(e,1) #提取均值图层
    e_slope <- subset(e,2)   #提取sen斜率
    e_MKtest <- subset(e,3)   #提取p值
    
    plot(e_mean)
    plot(e_slope)
    plot(e_MKtest)
    
    writeRaster(e_mean, "E:/CN/NDVI/e_Zs.tif", format="GTiff", overwrite=T)
    writeRaster(e_slope, "E:/CN/NDVI/e_slope.tif", format="GTiff", overwrite=T)
    writeRaster(e_MKtest, "E:/CN/NDVI/e_MKtest.tif", format="GTiff", overwrite=T)

![在这里插入图片描述][a7798cb751ec433192b63afa7bdaaf53.png]

图2 2010-2020年间艾比湖地区NDVI均值图层

![在这里插入图片描述][ca1c5b0c1ac44044b5f3fb11c5f28d23.png]

图3 R语言运行界面  
![在这里插入图片描述][b8cfc98b11a049bc89ff8baeaec51e7e.png]

图4 p值  
![在这里插入图片描述][03e79f81b8164450bddde8904737eb24.png]

图5 sen斜率  
![在这里插入图片描述][ce2140e76ec8421eb409911fe5206b62.png]

图6 Z统计量

R语言计算完slope和Z值后，根据这两个结果就可以进行NDVI趋势制图了。

一、变化趋势划分  
结合SNDVI和Z统计量划分NDVI变化趋势：  
1、slope  
\-0.0005~0.0005稳定区域  
大于或等于0.0005植被改善区域  
小于-0.0005为植被退化区域  
2、Z统计量

二、Slope划分  
置信水平0.05  
Z绝对值大于1.96显著  
Z绝对值小于等于1.96不显著  
Slope被划分为三级：  
SNDVI≤−0.0005 植被退化  
−0.0005≥SNDVI≥0.0005 植被生长稳定  
SNDVI≥0.0005 植被改善  
使用重分类（Reclassify）对slope进行划分  
由于slope.tif文件研究区范围外的值非空，所以在这里先裁剪了一下  
裁剪所用矢量和栅格数据坐标系需要一致，否则范围容易出错  
统一使用了WGS84地理坐标系作为空间参考  
使用Model builder构建地理处理流  
![在这里插入图片描述][0616681a82dc4dfab0fc46f92be56b39.png]

图7 重分类

三、Slope划分过程  
重分类结果：  
\-1退化  
0稳定  
1改善  
![在这里插入图片描述][ca0df52378eb46999833dbc145c366fa.png]

图8 重分类结果

四、Z值划分  
对Z值进行重分类，确定显著性  
|Zs|≤0.196 未通过95%置信度检验，不显著  
|Zs|≥0.196 通过95%置信度检验，显著  
![在这里插入图片描述][c5c6b91684e947528248215557743fbc.png]

图9 重分类  
五、Z值重分类  
重分类结果：  
1不显著  
2显著  
![在这里插入图片描述][b117d7e19fcf4c659fa4b7a647ffdbc7.png]

图10 重分类结果

六、变化趋势计算  
使用栅格计算器将Slope和Z值计算结果相乘，最后得到趋势变化划分  
\-2严重退化  
\-1轻微退化  
0稳定不变  
1轻微改善  
2明显改善  
![在这里插入图片描述][deacee97e0ce4c0882d7fad4bce1a001.png]

图11 栅格计算器相乘  
![在这里插入图片描述][f653446d31304509bfbaedd8e2cfcde8.png]

图12 arcgis计算NDVI趋势图

[0fdaefaf6d7f4761b5868abeec8ef557.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/e7cd702dafd342a1837f098ca57efd2f.png
[a7798cb751ec433192b63afa7bdaaf53.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/ec4eb10432e74e379b796531cfe48cfa.png
[ca1c5b0c1ac44044b5f3fb11c5f28d23.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/bd04e7f0f3e74aef8f6a385af69f4625.png
[b8cfc98b11a049bc89ff8baeaec51e7e.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/c2751619e0d946ce8705a3290348dc4f.png
[03e79f81b8164450bddde8904737eb24.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/333e6ba08fdf4ac99ffac4b19e49508e.png
[ce2140e76ec8421eb409911fe5206b62.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/c66ea3313a1b44c5a03e6c21dc237270.png
[0616681a82dc4dfab0fc46f92be56b39.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/169c98fc2fd442b9ac205442767482d3.png
[ca0df52378eb46999833dbc145c366fa.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/8ef08b14e1f840de978abfd6c0fdb959.png
[c5c6b91684e947528248215557743fbc.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/0ddd523a525c4f70b4a231ace8a75ad5.png
[b117d7e19fcf4c659fa4b7a647ffdbc7.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/685275abb35a4d0087c3b2b164277988.png
[deacee97e0ce4c0882d7fad4bce1a001.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/d9d07398597d40eca1eb820714002bd4.png
[f653446d31304509bfbaedd8e2cfcde8.png]: https://image.dandelioncloud.cn/pgy_files/images/2024/04/15/f2c0887396be4c5d8dede30bb73b6eea.png