Gimms NDVI3g v1 數(shù)據(jù)產(chǎn)品介紹

關(guān)于此數(shù)據(jù)集處理方式網(wǎng)上找了不少資料，一般都是用matlab或者是R來(lái)處理，matlab處理方法直接百度都能搜到很多，這里就不逐一介紹了，主要講一講如何用R進(jìn)行處理，得到整月最大值合成，或者是半月影像。

老板事先安排了數(shù)據(jù)下載工作，于是就在ecocast用wget下載好了數(shù)據(jù)，本以為完事了，老板又安排對(duì)數(shù)據(jù)進(jìn)行提取，得到半月影像數(shù)據(jù)，嘗試了gdal，讀取雖然沒(méi)報(bào)錯(cuò)，但是值不對(duì)，再加上電腦重裝系統(tǒng)沒(méi)有matlab，但是之前有用過(guò)R，于是嘗試使用R來(lái)處理數(shù)據(jù)。

在使用R之前也嘗試過(guò)GDAL來(lái)讀取，提取半月影像數(shù)據(jù)，使用gdalinfo看到數(shù)據(jù)一些信息，一共倆subdataset ndvi包含12個(gè)波段，percentile也是12個(gè)波段，percentile相當(dāng)于一個(gè)包含精度信息的數(shù)據(jù)集，里面不同區(qū)間的像元值代表了ndvi數(shù)據(jù)是如何得到的。使用gdal_translate提取的ndvi數(shù)據(jù)有些問(wèn)題，之前覺(jué)得可能ndvi數(shù)據(jù)與percentile有關(guān)聯(lián)，不能分開(kāi)讀取，老板認(rèn)為可能是與gdal或者netcdf版本有關(guān)系，這里也沒(méi)深究。在使用R處理數(shù)據(jù)前，首先介紹一個(gè)R包，gimms，這個(gè)是類(lèi)似于官方的一個(gè)R包，專(zhuān)門(mén)用來(lái)處理Gimms ndvi3g（此數(shù)據(jù)集分為v0和v1，分別對(duì)應(yīng)envi格式和netcdf4格式）數(shù)據(jù)，可以對(duì)數(shù)據(jù)進(jìn)行下載，提取，最大值合成。

在R中進(jìn)行包管理還是很簡(jiǎn)單的，直接install.packages("gimms","rgdal") ，rgdal是gimms包的依賴(lài)包。

Gimms NDVI3g v1 數(shù)據(jù)產(chǎn)品處理

接下來(lái)直接上代碼（本人一直用java，R只是接觸過(guò)幾次，所以代碼結(jié)構(gòu)比較混亂），由于提前下載好了數(shù)據(jù)，所以這里就直接對(duì)硬盤(pán)上數(shù)據(jù)進(jìn)行處理。

# downloadGimms() 可以使用這個(gè)函數(shù)來(lái)下載

#

library("gimms","rgdal")

#這個(gè)包還有些函數(shù)，可以對(duì)數(shù)據(jù)進(jìn)行精度控制，這里暫時(shí)用不到。

for (value in 1982:2015) {

    fn1 <- paste0("你的輸入路徑/ndvi3g_geo_v1_",value,"_","0106",".nc4")

    ndvi3g <- rasterizeGimms(x = fn1)

    sprintf("正在讀取...%s_%s",value,"0106")

    for (variable in 1:12) {

      name <- paste0("你的輸出路徑/ndvi3g_geo_v1_" , value ,"_","0106","-", variable, ".tif")

      sprintf("正在處理...%s",name)

      writeRaster(x = ndvi3g[[variable]],filename = name)

    }

}

for (value in 1981:2015) {

  fn2 <- paste0("你的輸入路徑/ndvi3g_geo_v1_",value,"_","0712",".nc4")

  ndvi3g <- rasterizeGimms(x = fn2)

  sprintf("正在讀取...%s_%s",value,"0712")

for (variable in 1:12) {

  name <- paste0("你的輸出路徑/ndvi3g_geo_v1_" ,value,"_","0712","-", variable, ".tif")

  sprintf("正在處理...%s",name)

  writeRaster(x = ndvi3g[[variable]],filename = name)

}

    }

最大值合成 月ndvi數(shù)據(jù)處理

gimms_files_date <- downloadGimms(x = as.Date("2000-01-01"),y = as.Date("2000-06-06"),dsn="E:/ndvi")

mvc <- monthlyComposite(ndvi3g_geo_v1_2000_0106,indices = monthlyIndices(gimms_files_date))

ndvi3g <- rasterizeGimms(x = gimms_files_date)

mvc <- monthlyComposite(ndvi3g,indices = monthlyIndices(gimms_files_date))

writeRaster(x = mvc,filename = "E:\ndvi\ndvi.tif")

image.png