Hi
I got a function to convert ndvi raster stack to CSVs. Each stack is divided into 100 subset, which is convert to csv. The code works for small raster stack, for large ones, I cannot load them into the memory, then it takes massive time to do the task.
is there a better way to do it?
The code is below:
require(gdal)
require(raster)
exportCSV <- function () {
? tif <- list.files(pattern='NDVI.tif$')
? wd <- getwd()
? ModisTile<-? substr(wd, nchar(wd)-5, nchar(wd))
? nImages <- length(tif)
? cat(paste("Stacking images ...", "\n"))
? s <- stack(tif)
? cat(paste("Loading values to RAM memory ...", "\n"))
#this step is skipped in case of large stacks, then it takes very long time
? s <- readAll(s)
#create the subsets bounding coordinates
? borders <- extent(s)
? Xmin <- borders at xmin
? Xmax <- borders at xmax
? Ymin <- borders at ymin
? Ymax <- borders at ymax
? xIncreament <-(Xmax-Xmin)/10
? yIncreament <-(Ymax-Ymin)/10
? cat(paste("Subsetting and writing NDVI values ...", "\n"))
? for (i in 1:10) {
??? for (j in 1:10) {
????? clip_xmin <- Xmin + xIncreament*(i-1)
????? clip_xmax <- Xmin + xIncreament*i
????? clip_ymin <- Ymin + yIncreament*(j-1)
????? clip_ymax <- Ymin + yIncreament*j
????? c_xmin <- format(round(clip_xmin,6), nsmall=6)
????? c_xmax <- format(round(clip_xmax,6), nsmall=6)
????? c_ymin <- format(round(clip_ymin,6), nsmall=6)
????? c_ymax <- format(round(clip_ymax,6), nsmall=6)????
?????
????? subset <- extent(c(clip_xmin, clip_xmax, clip_ymin, clip_ymax))
????? c <- crop(s, subset)
????? p <- as.data.frame(rasterToPoints(c))
????? csvName <- paste0(ModisTile, "_Xmin_",c_xmin, "_Xmax_",c_xmax, "_Ymin_",c_ymin, "_Ymax_",c_ymax,".csv")
????? cat(paste("Writing Subset... MOIDS Tile:", ModisTile,", X", i, "Y", j, "\n"))
????? write.table(p, csvName, row.names=F, sep=";", dec=".")
??? }
? }
}
Best,
Mohammad?PhD Candidate Institute of Crop Science and Resource Protection - Crop Science Research Group
Katzenburgweg 5 - 53115 Bonn - Germany
Tel.: +49 (0) 228 73 3258?????? Fax: +49 (0) 228 73 2870
abdelrazek at uni-bonn.de??????? http://www.lap.uni-bonn.de
Converting large RasterStack to CSVs fast
2 messages · Mohammad Abdel-Razek, Éder Comunello
Hello, Mohammad!
You may have some improvement in performance avoiding "for statements" and
using a "vectorized" code. You could try something like the code below.
If you can test with your data, i would appreciate if you inform the
results.
### <code r>
require(rgdal); require(raster)
getwd()
### download some data to test
getData("worldclim", var = "tmin", res = 10) ### tmin
fn <- dir("wc10", patt=".bil$", full=T)
fn <- fn[order(nchar(fn), fn)]; fn
# [1] "wc10/tmin1.bil" "wc10/tmin2.bil" "wc10/tmin3.bil" ...
### read images
s <- stack(fn) ### dimensions : 900, 2160, 1944000, 12 (nrow, ncol,
ncell, nlayers)
fromDisk(s)
### extents of subsets
bor <- extent(s); bor
res <- 45 ### subsets resolution
X <- unique(c(seq(bor at xmin, bor at xmax, by=res), bor at xmax)); X
Y <- unique(c(seq(bor at ymin, bor at ymax, by=res), bor at ymax)); Y
ext <- cbind(expand.grid(Xmin=X[-length(X)], Ymin=Y[-length(Y)]),
expand.grid(Xmax=X[-1], Ymax=Y[-1]))[,c(1,3,2,4)]
head(ext); nrow(ext)
plot(s, 1)
system.time(
sapply(1:nrow(ext), function(i) {
mask <- ext[i,]
subset <- with(mask, extent(c(Xmin, Xmax, Ymin, Ymax)))
plot(subset, add=T)
text(rowMeans(mask[,1:2]), rowMeans(mask[,3:4]), lab=i)
c <- crop(s, subset)
write.table(as.data.frame(rasterToPoints(c)), paste0("p",i,".txt"), )
}))
# user system elapsed
# 213.79 7.00 224.94
txt <- dir(patt="^p[0-9]+.txt$")
txt <- txt[order(nchar(txt), txt)]; txt
# [1] "p1.txt" "p2.txt" "p3.txt" "p4.txt" ...
### </code>
Cheers,
?der Comunello <c <comunello.eder at gmail.com>omunello.eder at gmail.com>
Dourados, MS - [22 16.5'S, 54 49'W]
?der Comunello <c <comunello.eder at gmail.com>omunello.eder at gmail.com>
Dourados, MS - [22 16.5'S, 54 49'W]
2015-05-15 5:43 GMT-04:00 Mohammad Abdel-Razek via R-sig-Geo <
r-sig-geo at r-project.org>:
Hi
I got a function to convert ndvi raster stack to CSVs. Each stack is
divided into 100 subset, which is convert to csv. The code works for small
raster stack, for large ones, I cannot load them into the memory, then it
takes massive time to do the task.
is there a better way to do it?
The code is below:
require(gdal)
require(raster)
exportCSV <- function () {
tif <- list.files(pattern='NDVI.tif$')
wd <- getwd()
ModisTile<- substr(wd, nchar(wd)-5, nchar(wd))
nImages <- length(tif)
cat(paste("Stacking images ...", "\n"))
s <- stack(tif)
cat(paste("Loading values to RAM memory ...", "\n"))
#this step is skipped in case of large stacks, then it takes very long time
s <- readAll(s)
#create the subsets bounding coordinates
borders <- extent(s)
Xmin <- borders at xmin
Xmax <- borders at xmax
Ymin <- borders at ymin
Ymax <- borders at ymax
xIncreament <-(Xmax-Xmin)/10
yIncreament <-(Ymax-Ymin)/10
cat(paste("Subsetting and writing NDVI values ...", "\n"))
for (i in 1:10) {
for (j in 1:10) {
clip_xmin <- Xmin + xIncreament*(i-1)
clip_xmax <- Xmin + xIncreament*i
clip_ymin <- Ymin + yIncreament*(j-1)
clip_ymax <- Ymin + yIncreament*j
c_xmin <- format(round(clip_xmin,6), nsmall=6)
c_xmax <- format(round(clip_xmax,6), nsmall=6)
c_ymin <- format(round(clip_ymin,6), nsmall=6)
c_ymax <- format(round(clip_ymax,6), nsmall=6)
subset <- extent(c(clip_xmin, clip_xmax, clip_ymin, clip_ymax))
c <- crop(s, subset)
p <- as.data.frame(rasterToPoints(c))
csvName <- paste0(ModisTile, "_Xmin_",c_xmin, "_Xmax_",c_xmax,
"_Ymin_",c_ymin, "_Ymax_",c_ymax,".csv")
cat(paste("Writing Subset... MOIDS Tile:", ModisTile,", X", i, "Y",
j, "\n"))
write.table(p, csvName, row.names=F, sep=";", dec=".")
}
}
}
Best,
Mohammad PhD Candidate Institute of Crop Science and Resource Protection
- Crop Science Research Group
Katzenburgweg 5 - 53115 Bonn - Germany
Tel.: +49 (0) 228 73 3258 Fax: +49 (0) 228 73 2870
abdelrazek at uni-bonn.de http://www.lap.uni-bonn.de
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo