r - 如何使用下面提供的代码创建 R 循环?
问题描述
请,我需要帮助来创建一个循环,该循环将在包含 R 中的 483 个文件的 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 个文件的栅格
解决方案
您可以将代码包装在一个函数中,然后lapply
通过 hdf 路径。这样,如果您的循环太慢,将很容易并行化它。你可以试试这个:
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], paste0(name,'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"
在lapply
我的电脑上,该功能在 45 秒内运行。例如,您可以通过从包中lapply
替换它来轻松实现并行化。仅仅 2 个文件是不值得的,但是如果你有数百个文件,你可以大大加快这个过程: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()
使用sfLapply
此功能需要 20 秒才能运行。这是一个很好的改进
推荐阅读
- c# - 具有 azure AD 身份验证的长寿网页,超时
- c# - 在 C# win 应用程序中创建电子邮件通知功能
- sql - SQL - 按不同列的值为空的列对两行进行分组
- add-on - 面临创建插件、Visual Studio、G1ANT 的问题
- implicit - Gekko 和 CoolProp
- php - Phpmyadmin 显示 mbstring 丢失
- c# - 我们什么时候应该在 C# 中使用默认接口方法?
- c# - 为什么简单的 await Task.Delay(1) 可以实现并行执行?
- postgresql - Postgres 查询优化器在按条件添加另一个订单后生成错误计划
- visual-studio - std::shared_ptr 的读取访问冲突