NDVI数据集提取植被覆盖度FVC

最新推荐文章于 2024-09-20 16:06:30 发布

彭博锐

最新推荐文章于 2024-09-20 16:06:30 发布

阅读量1.9k

点赞数 29

文章标签： python R ENVI QGIS 笔记学习

本文链接：https://blog.csdn.net/pengbr/article/details/140236029

版权

植被覆盖度FVC

植被覆盖度（Foliage Vegetation Cover，FVC）是指植被冠层覆盖地表的面积比例，通常用来描述一个区域内植被的茂密程度或生长状况。它是生态学、环境科学以及地理信息系统等领域的重要指标，对于理解地表能量平衡、水文循环、碳循环等生态过程至关重要。

植被覆盖度可以通过地面调查、遥感和地理信息系统等技术手段来测量。在遥感应用中，通常使用植被指数（如归一化植被指数NDVI）来估算植被覆盖度。这些指数反映了植被的健康状况和生物量，是监测植被变化、评估生态系统健康和进行环境变化研究的重要工具。

FVC的计算公式（参考文献Cui et al.,2023）：
NDVImin为累计像元值5%的值，NDVImax为累计像元值95%的值
FVC的取值范围为[0,1]
在这里插入图片描述
表格参考自知乎用户@24K知识星球。

NDVI数据集

本次NDVI数据集选择的是GIMMS-3G+数据集，下载网址：GIMMS-3G+
我下载了2022年下半年的数据，数据格式是nc文件

首先用R语言将nc文件拆分成按月划分的tif文件

library(raster)
library(ncdf4)
ndvi<-brick("ndvi3g_geo_v1_2_2022_0712.nc4",varname="ndvi")
writeRaster(ndvi,filename = ndvi@data@names,bylayer=T,format="GTiff",overwrite=TRUE)

这样就获得了按月划分的tif文件。比如7月的文件为X2022.07.tif
由于比例系数是0.0001，所以每个波段的数值需要乘以0.0001
利用R语言将像元数值缩小10000倍。

library(raster)
library(sp)

# 打开.tif文件
input_file <- "X2022.07.tif"
output_file <- "X2022.07.scaled.tif"
raster_data <- raster(input_file)

# 获取波段的最小值和最大值
min_value <- minValue(raster_data)
max_value <- maxValue(raster_data)

# 设置新的最小值和最大值范围
new_min <- min_value/10000
new_max <- max_value/10000

# 调整波段范围
adjusted_raster <- ((raster_data - min_value) / (max_value - min_value)) * (new_max - new_min) + new_min

# 保存修改后的.tif文件
writeRaster(adjusted_raster, output_file, format="GTiff")

利用QGIS查看下归一化后的NDVI文件如下：
在这里插入图片描述
海洋部分的无数据nodata，冰覆盖的数值为-0.3

由于冰覆盖的数值较多，我不想让该数值影响到NDVImin和NDVImax值的计算，故写了一个python脚本将-0.3的值变成nodata

import rasterio

# 打开NDVI TIF文件
with rasterio.open('X2022.07.scaled.tif') as src:
    # 读取NDVI数据
    ndvi = src.read(1)

    # 判断是否有Nodata
    nodata = src.nodata
    if nodata is None:
        print('没有空数据')
    else:
        print(nodata)

    # 将NDVI值为-0.3的像素赋值为NoData
    ndvi[(ndvi == -0.3)] = nodata

    # 复制源文件的元数据
    profile = src.profile

    # 更新NoData值
    profile['nodata'] = nodata

    # 保存新的TIF文件
    with rasterio.open('X2022.07.scaled2.tif', 'w', **profile) as dst:
        dst.write_band(1, ndvi)

去除-0.3后的tif文件QGIS展示：
在这里插入图片描述
然后再利用python计算NDVImin和NDVImax

import rasterio
import numpy as np

# 打开NDVI栅格数据
with rasterio.open('X2022.07.scaled2.tif') as src:
    ndvi = src.read(1)  # 假设NDVI是在第一个波段

# 判断是否有Nodata
nodata = src.nodata
if nodata is None:
    print('没有空数据')
else:
    print(nodata)

# 排除NoData值。NoData的值需要排除，不然会影响到NDVImin和NDVImax的计算结果。
ndvi_valid = ndvi[(ndvi != nodata)]
print(ndvi_valid)

# 计算5%和95%作为最小值和最大值
ndvi_min = np.percentile(ndvi_valid, 5)
ndvi_max = np.percentile(ndvi_valid, 95)

print(f"NDVI最小值: {ndvi_min}")
print(f"NDVI最大值: {ndvi_max}")

终端显示：

-3.3999999521443642e+38
[-0.0356 -0.0511 -0.0456 ...  0.4683  0.5181  0.565 ]
NDVI最小值: 0.0869000032544136
NDVI最大值: 0.8664000034332275

利用ENVI软件的快速统计功能同样也可以获取累计百分数为5%和95%的值，结果和python运算的结果差不多。
在这里插入图片描述
接下来，就可以利用QGIS的栅格计算器功能获取FVC了。
QGIS中的栅格计算器在Raster（栅格）选项栏中。
根据FVC的计算公式输入栅格计算表达式：

 ( "X2022.07.scaled2@1"  - 0.0869000032544136 )  /  ( 0.8664000034332275 - 0.0869000032544136)

在这里插入图片描述
这样，就获得了FVC的图如下。可以看到取值不在0-1之间，所以需要将小于0的值都赋值为0，大于1的部分赋值为1，其余部分为其本身的值。

再次利用栅格计算器，将FVC的值调整到0-1之间。
栅格计算表达式：

("X202207FVC@1" < 0) * 0 + ("X202207FVC@1" > 1) * 1 + ("X202207FVC@1" >= 0 AND "X202207FVC@1" <= 1) * "X202207FVC@1"

在这里插入图片描述
现在就将FVC的值调整到0-1了。我还对图层的样式做了调整，突出显示了FVC<=0.2的地区

本文参考文献

1、Cui, B., Gui, D., Liu, Q., Abd‐Elmabod, S. K., Liu, Y., & Lu, B. (2024). Distribution and growth drivers of oases at a global scale. Earth’s Future, 12, e2023EF004086. https://doi.org/10.1029/2023EF004086
2、Pinzon, J.E., E.W. Pak, C.J. Tucker, U.S. Bhatt, G.V. Frost, and M.J. Macander. 2023. Global Vegetation Greenness (NDVI) from AVHRR GIMMS-3G+, 1981-2022. ORNL DAAC, Oak Ridge, Tennessee, USA. https://doi.org/10.3334/ORNLDAAC/2187