[R-sig-Geo] Converting large RasterStack to CSVs fast
Mohammad Abdel-Razek
kimofos at yahoo.com
Fri May 22 13:51:29 CEST 2015
Hello Comunello!
Thanks for the reply. Your vectorized code is faster than the one I have, by a margin of 5%, which is an improvement. I think the limiting step is clipping. Is there any paralleled version of clip?
Mohammad
From: "r-sig-geo-request at r-project.org" <r-sig-geo-request at r-project.org>
To: r-sig-geo at r-project.org
Sent: Saturday, May 16, 2015 12:00 PM
Subject: R-sig-Geo Digest, Vol 141, Issue 16
Send R-sig-Geo mailing list submissions to
r-sig-geo at r-project.org
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
or, via email, send a message with subject or body 'help' to
r-sig-geo-request at r-project.org
You can reach the person managing the list at
r-sig-geo-owner at r-project.org
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-Geo digest..."
Today's Topics:
1. Re: Issue with ogrInfo (?der Comunello)
2. Re: Converting large RasterStack to CSVs fast (?der Comunello)
3. create a shapefile (Gustavo Dalposso)
4. Re: create a shapefile (Felinto COSTA)
----------------------------------------------------------------------
Message: 1
Date: Fri, 15 May 2015 07:41:57 -0400
From: ?der Comunello <comunello.eder at gmail.com>
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] Issue with ogrInfo
Message-ID:
<CABmC8g=ZAFGTu=RJ6BaZjD76gstYwrRJ=rRPqPCuE6WDKVkaqQ at mail.gmail.com>
Content-Type: text/plain; charset="UTF-8"
Hello,
I think the problem is that your link is recovering only html code from
github site. It's necessary to modify it to get the "real" JSON file.
### <code r>
sapply(c("rgeos", "maptools", "rgdal"), require, char=T)
url0 <- "
https://github.com/kjhealy/uk-elections/blob/master/maps/topo_wpc.json"
download.file(url0, dest=basename(url0), mode="wb")
# downloaded 26 KB
uk.map <- readOGR(basename(url0), "wpc")
# Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv
= use_iconv, :
# Cannot open file
url1 <- "
https://raw.githubusercontent.com/kjhealy/uk-elections/master/maps/topo_wpc.json
"
download.file(url1, dest=basename(url0), mode="wb")
# downloaded 2.3 MB
uk.map <- readOGR(basename(url1), "wpc")
# OGR data source with driver: GeoJSON
# Source: "topo_wpc.json", layer: "wpc"
# with 632 features
# It has 2 fields
plot(uk.map)
### </code>
?der Comunello <c <comunello.eder at gmail.com>omunello.eder at gmail.com>
Dourados, MS - [22 16.5'S, 54 49'W]
[[alternative HTML version deleted]]
------------------------------
Message: 2
Date: Fri, 15 May 2015 12:35:44 -0400
From: ?der Comunello <comunello.eder at gmail.com>
To: "r-sig-geo at r-project.org" <r-sig-geo at r-project.org>
Subject: Re: [R-sig-Geo] Converting large RasterStack to CSVs fast
Message-ID:
<CABmC8gmd8aV6WBRpXRGnwbmRht8Z_MX3WWd9=N999RzoaZVoRA at mail.gmail.com>
Content-Type: text/plain; charset="UTF-8"
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]]
------------------------------
Message: 3
Date: Fri, 15 May 2015 19:09:07 +0000
From: Gustavo Dalposso <gustavodalposso at hotmail.com>
To: "r-sig-geo at r-project.org" <r-sig-geo at r-project.org>
Subject: [R-sig-Geo] create a shapefile
Message-ID: <SNT152-W5068688FE73BB915BA33D2BBC70 at phx.gbl>
Content-Type: text/plain; charset="iso-8859-1"
Hello friends of R
I have a map in BMP format (a figure).
The map consists of 14 areas.
I have the location of two coordinates. (x1, y1) and (x2, y2)
See the attached image.
Question: Is it possible to create a shapefile? Which package I use?
Atte.
Gustavo Henrique Dalposso
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20150515/df1e9d08/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: map.JPG
Type: image/jpeg
Size: 11778 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20150515/df1e9d08/attachment-0001.jpe>
------------------------------
Message: 4
Date: Fri, 15 May 2015 16:44:13 -0300
From: Felinto COSTA <incorpld at onda.com.br>
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] create a shapefile
Message-ID: <55564C8D.1040504 at onda.com.br>
Content-Type: text/plain; charset="UTF-8"
Gustavo.
Se vc. tiver esse mapa em formato CAD (dxf ou dwg) pode usar o Cad2Shape
(http://www.guthcad.com.au/cad2shape.htm).
Ele permite algumas utiliza?s em sua vers?"free".
Se n?tiver, tente passar para os formatos de CAD usando, talvez, o
Img2Cad (http://www.img2cad.com/).
J?sei apenas o primeiro.
Felinto COSTA
On 15/5/2015 16:09, Gustavo Dalposso wrote:
> Hello friends of R
>
>
> I have a map in BMP format (a figure).
> The map consists of 14 areas.
> I have the location of two coordinates. (x1, y1) and (x2, y2)
> See the attached image.
>
> Question: Is it possible to create a shapefile? Which package I use?
>
> Atte.
> Gustavo Henrique Dalposso
>
>
> _______________________________________________
> 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]]
------------------------------
Subject: Digest Footer
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
------------------------------
End of R-sig-Geo Digest, Vol 141, Issue 16
******************************************
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list