[R-sig-Geo] Converting large RasterStack to CSVs fast

[É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
>

	[[alternative HTML version deleted]]