走天涯徐小洋

V1

2023/03/15阅读:18主题:前端之巅同款

NDVI值域恢复

如何将2000-2020年中国30米年最大NDVI数据集0-255的值域恢复正常

有同学反映说不知道如何把0-255值域的NDVI数据恢复为正常的NDVI值域,其实这个用R很容易实现的,在这里给出代码。

数据来源

数据来源于董金玮等在国家生态科学数据中心共享的2000-2020年中国30米年最大NDVI数据集。

数据集下载网址:http://www.nesdc.org.cn/sdo/detail?id=60f68d757e28174f0e7d8d49

在网页上有数据说明文档下载按钮,下载数据说明文档可以查看数据使用说明
在网页上有数据说明文档下载按钮,下载数据说明文档可以查看数据使用说明

数据引用方式

董金玮;周岩;尤南山. 2000-2020年中国30米年最大NDVI数据集[DB/OL]. 国家生态科学数据中心, 2021. doi:10.12199/nesdc.ecodb.rs.2021.012. cstr:15732.11.nesdc.ecodb.rs.2021.012

使用R语言恢复NDVI正常值域

在说明文档中给出了恢复NDVI值域的公式:

说明文档中关于恢复正常值域的说明
说明文档中关于恢复正常值域的说明

R语言代码

先把数据读取进来,在这里读取了一块样例数据,根据预览图和统计结果可以发现,目前数据值域为0-255

  • plot()绘制了一下草图
  • summary()查看读取的栅格值域
library(terra)
ltndvi = rast("./NDVImax2015/NDVImax2015.tif")
plot(ltndvi)
summary(ltndvi)
预览读取的数据
预览读取的数据
NDVI的最大最小值、一分位数、中值、均值、三分位数、最大值、空值个数
NDVI的最大最小值、一分位数、中值、均值、三分位数、最大值、空值个数

单核运行栅格四则运算

按照说明文档中的公式,再除以一万,即可将NDVI值域转换为[-1,1],对于R语言的栅格计算,直接写一个四则运算公式即可完成。不过R里面默认的四则运算是单核计算,计算速度较慢。下面是代码和计算结果,同样也进行了绘图和统计计算结果。

t1=proc.time()
ltndvi2 = (ltndvi*12000/255-2000)/10000
t2=proc.time()
t=t2-t1
print(paste0('执行时间:',t[3][[1]],'秒'))
plot(ltndvi2)
summary(ltndvi2)
已转换为正常值域
已转换为正常值域
计算结果,单核运行15秒
计算结果,单核运行15秒

并行四则运算

先把公式改成一个自定义函数,然后进行栅格计算,并行,由于下面的代码并行开了28个核,数据要被分成28份,然后并行计算,但是实际数据量并不大,这样就导致在分拆数据过程花费了大量时间,实际并行计算并没有起到加速作用。

由于这只是用很小的研究区进行了一个实验,在大数据量的时候并行还是很有价值的。

fun_ndvi =  function(x){
  a1 = (x*12000/255-2000)/10000
  return(a1)
}

t1=proc.time()
ltndvi3 = app(ltndvi, fun_ndvi, cores=28)
t2=proc.time()
t=t2-t1
print(paste0('执行时间:',t[3][[1]],'秒'))
plot(ltndvi3)
summary(ltndvi3)
可能是开的核太多,数据量太小,反而速度降低了很多
可能是开的核太多,数据量太小,反而速度降低了很多

这样就完成了NDVI值域的转换。

更多阅读

  1. 对比R语言Raster包和Terra包栅格计算
  2. R语言栅格时间序列回归分析——整体和逐像元计算,并行计算,快到飞起!

分类:

数学

标签:

数学编程

作者介绍

走天涯徐小洋
V1