如何使用下面提供的代码创建 R 循环？

编程问答 2022-05-04

问题描述

拜托，我需要帮助来创建一个循环，该循环将在包含 483 个 R 文件的 hdflist 上执行以下代码中所示的计算。我添加了一个包含两个 .hdf 文件和 shapefile 的链接以供试用。该代码对于单个 .hdf 文件似乎工作得很好，但我仍在为循环而苦苦挣扎。谢谢

download files from here https://beardatashare.bham.ac.uk/getlink/fi2gNzWbuv5H8Gp7Qg2aemdM/

# import .hdf file into R using get_subdatasets to access the subsets in the file`
sub <- get_subdatasets("MOD13Q1.A2020353.h18v08.006.2021003223721.hdf")

# convert red and NIR subsets and save them as raster`
gdalwarp(sub[4],'red_c.tif')
gdalwarp(sub[5],'NIR_c.tif')

# import red and NIR raster back into R`
# scale the rater while at it`

r_r=raster('red_c.tif') * 0.0001 
r_N=raster('NIR_c.tif') * 0.0001 

# calculate sigma using (0.5*(NIR+red))`
sigma <- (0.5*(r_N+r_r))

# calculate knr using exp((-(NIR-red)^2)/(2*sigma^2))`
knr <- exp((-(r_N-r_r)^2)/(2*sigma^2))

# calculate kndvi using (1 - knr) / (1 + knr)`
kndvi <- (1 - knr) / (1 + knr)

# import shapefile into R`
shp=readOGR(".","National_Parks")
options(stringsAsFactors = FALSE)

#change crs of shapefile to crs of one of the rasters`
shp2 <- spTransform(shp,crs(kndvi))

# use extent to crop/clip raster`
## set extent`
e <- extent(910000,980000,530000,650000)

## clip using crop function`

crop_kndvi <- crop(kndvi,e)

# mask raster using the shapefile`
kndvi_mask <- mask(crop_kndvi,shp2)

然后将 kndvi_mask 保存为 483 个文件的光栅

final kndvi raster for just one hdf file analyzed

解决方法

您可以将代码包装在一个函数中，然后在 hdf 路径上lapply。这样，如果您的循环太慢，则很容易将其并行化。 你可以试试这个：

library(gdalUtils)
library(raster)
library(rgdal)
#set the directory where you have .hdf files. In my case I downloaded your data in "D:/download"
setwd("D:/download")
#function to save the masked index in your current working directory
#the final files name will depend on the name of the input hdf files
myfun <- function(path){
  name <- basename(tools::file_path_sans_ext(path))
  sub <- get_subdatasets(path)
  gdalwarp(sub[4],paste0(name,'_red_c.tif'))
  gdalwarp(sub[5],'NIR_c.tif'))
  r_r=raster(paste0(name,'_red_c.tif')) * 0.0001 
  r_N=raster(paste0(name,'NIR_c.tif')) * 0.0001 
  
  sigma <- (0.5*(r_N+r_r))
  knr <- exp((-(r_N-r_r)^2)/(2*sigma^2))
  kndvi <- (1 - knr) / (1 + knr)
  crop_kndvi <- crop(kndvi,e)
  kndvi_mask <- mask(crop_kndvi,shp2,filename=paste0(name,"_kndvi_mask.tif"))
}



    #list the hdf file in your current working directory. Thanks to setwd("D:/download") there is no need to specify the path argument of list.files().
   b#however for the for peace of mind:
    hdf <- list.files(path=getwd(),pattern = "hdf",full.names = T)
    #since your shop is always the same you could keep this part out of the function
    shp=readOGR(".","National_Parks")
    options(stringsAsFactors = FALSE)
    shp2 <- spTransform(shp,"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m
                        +no_defs ")
    e <- extent(910000,980000,530000,650000)
    #now run your function across the hdf files path
    lapply(hdf,myfun)

现在在您的工作目录中您可以找到所有已保存的 if

   list.files(pattern = "tif")
[1] "MOD13Q1.A2020337.h18v08.006.2020358165204_kndvi_mask.tif"
[2] "MOD13Q1.A2020337.h18v08.006.2020358165204_red_c.tif"     
[3] "MOD13Q1.A2020337.h18v08.006.2020358165204NIR_c.tif"      
[4] "MOD13Q1.A2020353.h18v08.006.2021003223721_kndvi_mask.tif"
[5] "MOD13Q1.A2020353.h18v08.006.2021003223721_red_c.tif"     
[6] "MOD13Q1.A2020353.h18v08.006.2021003223721NIR_c.tif"

在我的 PC 上使用 lapply 时，该函数在 45 秒内运行。 例如，您可以通过将 lapply 替换为 sfLapply 包中的 snowfall 来轻松并行化 library(snowfall) #open cluster with as many node as hdf file sfInit(parallel=TRUE,cpus=length(hdf)) # Load the required packages inside the cluster sfLibrary(raster) sfLibrary(rgdal) sfLibrary(gdalUtils) sfExportAll() system.time(sfLapply(hdf,myfun)) sfStop() 。对于仅 2 个文件来说这是不值得的，但如果您有数百个文件，则可以大大加快进程：

sfLapply

使用 df.set_index("Date",inplace=True) df=df.T df 此函数需要 20 秒才能运行。这是一个很好的改进

,

您可以通过以下方式使用 terra 做到这一点。 terra 是 raster 的替代品；它更快，更通用。例如，使用 terra，您可以跳过 gdalwarp 步骤。

您可以编写一个大 for-loop，但我更喜欢使用函数，然后在循环或 lapply 中调用它们。

此外，代替您的光栅代数方法，将 kndvi 计算包装到其自己的函数中并与 lapp 一起使用可能会更有效。我认为这是一种更好的方法，因为代码更清晰，并且允许您重复使用 kndvi 函数。

library(terra)
parks <- vect("National_Parks.shp")
parks <- project(parks,"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m")
e <- ext(910000,650000)

lapp 使用的 kndvi 函数

kndvi <- function(red,NIR) {
    red <- red * 0.0001 
    NIR <- NIR * 0.0001 
    sigma <- (0.5 * (NIR + red))
    knr <- exp((-(NIR-red)^2)/(2*sigma^2))
    (1 - knr) / (1 + knr)
}

主要功能。请注意，我在其他函数之前使用了 crop ；从而节省了大量不必要的处理。

fun <- function(f) {
    outf <- gsub(".hdf$","_processed.tif",f) 
    # if file.exists(outf) return(rast(outf))
    r <- rast(f)[[4:5]]
    # or r <- sds(f)[4:5]
    r <- crop(r,e)
    kn <- lapp(r,kndvi)
    name <- substr(basename(f),9,16)
    mask(kn,parks,filename=outf,overwrite=TRUE,names=name)
}

获取文件名并将函数与循环或 lapply 一起使用，如 Elia 所示。

ff <- list.files(pattern="hdf$",full=TRUE)

x <- list()
for (i in 1:length(ff)) {
    print(ff[i]); flush.console()
    x[[i]] <- fun(ff[i])
}
z <- rast(x)
z

#class       : SpatRaster 
#dimensions  : 518,302,2  (nrow,ncol,nlyr)
#resolution  : 231.6564,231.6564  (x,y)
#extent      : 909946.2,979906.4,530029.7,650027.7  (xmin,xmax,ymin,ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#sources     : MOD13Q1.A2020337.h18v08.006.2020358165204_processed.tif  
#              MOD13Q1.A2020353.h18v08.006.2021003223721_processed.tif  
#names       :     A2020337,A2020353 
#min values  : 0.0007564131,0.0028829363 
#max values  :    0.7608207,0.7303495

这在我的计算机上每个文件大约需要 1 秒。

或者作为您要求的 for-loop：

ff <- list.files(pattern="hdf$",full=TRUE)
for (f in ff) {
    print(f); flush.console()
    outf <- gsub(".hdf$",f) 
    r <- rast(f)[[4:5]]
    r <- crop(r,names=name)
}

outf <- list.files(pattern="_processed.tif$",full=TRUE)
x <- rast(outf)
, 
hdf_files <- list.files("foldername",pattern = ".hdf")

for(f in files) { ... }

对于保存代码，您可以使用字符串“f”为文件创建一个名称，这样它们就不会相互重叠保存。

hdfloopsrrraster