[R-sig-Geo] gstat memory

[Dylan Beaudette]

On Friday 30 October 2009, Murray Richardson wrote:
> Hello R-SIG-GEO list,
>
> I know this has come up before but I am having an ongoing memory problem
> with the gstat package (gstat out of dynamic memory) that I can't seem
> to solve.
>
> I am using R to interpolate DEMs from LiDAR xyz point files and mosaic
> them together via RSAGA.  The script uses a loop to load in each xyz
> point file, and interpolate over a regular lattice of points using idw
> from gstat.  Although large, the computational requirements within each
> iteration of the loop should be well within my system's ability so it
> seems like it is a cumulative effect (note I can restart the process at
> the last loop that triggered the error and the iteration is
> successful).  I am removing temp objects and running gc() at the end of
> every loop.
>
> It proceeds normally for about 5-10 iterations and physical memory use
> on my system (VISTA 64 bit, 12GB RAM) gradually increases over time
> until I get the "gstat out of memory" error. Here is the relevant
> portion of the loop:
>

Hi,

If you aren't constrained to using R/SAGA, GRASS should be able to most of the 
interpolation in a single pass. Check out the modules v.surf.rst and 
v.surf.idw. I have head that people have routinely used v.surf.rst to 
interpolate point clouds 10-60Gb in size.

Cheers,
Dylan


>
>    path<-"E:/LidarData/1_Ground_First_Return/UTM17/"
>    filen<-paste(path, centre_tile$TILECODE, ".xyz", sep="")
>    xyzi<-read.table(filen, sep=",", header=F)
>    names(xyzi)<-c("x","y","z","i")
>
>    for(j in 1:length(tile_group_names)){
>        filen<-paste(path, tile_group_names[j], ".xyz", sep="")
>        tmp_xyzi<-read.table(filen, sep=",", header=F)
>        names(tmp_xyzi)<-c("x","y","z","i")
>        xyzi<-rbind(xyzi,tmp_xyzi)
>
>    }
>        coordinates(xyzi)=~x+y
>    # note grid_coords is just the regular lattice that is created from
> the current tile coordinates
>    coordinates(grid_coords)=~x+y
>
>    # do the interpolation
>    interp<-idw(z~1, xyzi, grid_coords, nmax=4,maxdist=2, idp=1.0)
>
>    names(interp)<-c("z","var")
>    slot(interp,"data")<-data.frame(slot(interp,"data"))
>    finalSPntsDF<-SpatialPointsDataFrame(interp, data.frame(interp$z),
>                proj4string = CRS(as.character(NA)), match.ID = TRUE)
>
>
>    finalSPDF<-SpatialPixelsDataFrame(finalSPntsDF, data.frame(interp$z),
> tolerance = sqrt(.Machine$double.eps),
>      proj4string = CRS(as.character(NA)))
>
> ...
> On a related note - I have tried using the 64 bit version of Revolution
> but unfortunately the gstat package has not been ported.  Has anyone
> contemplated or begun a 64 bit port of this package?
>
> Thanks
>
> Murray
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341