Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

This is a process that takes just a few seconds in GIS software. My attempt to do it in R uses a large amount of memory then fails. Is there something wrong in my code, or is this just something R cannot do? I have read R can work inside Grass, can I use a Grass function from inside R?

library(raster)

# I have many environmental rasters in this format
new_r <- raster(ncol=615, nrow=626, xmn=-156.2, xmx=-154.8, ymn=18.89, ymx=20.30)
res(new_r) <- 0.00225
projection(new_r) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0"

R> new_r ### not too big with a few hundred cells per side
class       : RasterLayer 
dimensions  : 627, 622, 1  (nrow, ncol, nlayers)
ncell       : 389994 
resolution  : 0.00225, 0.00225  (x, y)
projection  : +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs +towgs84=0,0,0 
extent      : -156.2, -154.8, 18.89, 20.3  (xmin, xmax, ymin, ymax)
values      : none

# I get the DEM at much higher resolution (zipfile is 182Mb)
zipurl <- "ftp://soest.hawaii.edu/coastal/webftp/Hawaii/dem/Hawaii_DEM.zip"
DEMzip <- download.file(zipurl, destfile = "DEMzip")
unzip("DEMzip", exdir = "HIDEM")
HIDEM <- raster("HIDEM/hawaii_dem")

R> HIDEM ### 10m resolution, file is way too big
class       : RasterLayer 
dimensions  : 15067, 13136, 1  (nrow, ncol, nlayers)
ncell       : 197920112 
resolution  : 10, 10  (x, y)
projection  : +proj=utm +zone=5 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 
extent      : 179066, 310426, 2093087, 2243757  (xmin, xmax, ymin, ymax)
values      : HIDEM/hawaii_dem 
min value   : 0 
max value   : 4200 

# the following line fails (after a long time)
new_HIDEM <- projectRaster(HIDEM, new_r)
share|improve this question
    
Just curious, what is the package you are using? –  djq Feb 23 '11 at 1:23
    
@celenius: this package is called raster –  J. Winchester Feb 23 '11 at 2:19

3 Answers 3

up vote 7 down vote accepted

From my look at the source, raster looks to guess if the dataset fits into memory, and if so, perform the operation in memory, otherwise on disk. You can force it to perform the calculation by explicitly setting chunksize (cells to process at a time) and maxmemory (maximum number of cells to read into memory):

setOptions(chunksize = 1e+04, maxmemory = 1e+06)

Alternatively, you could perform the transformation with GDAL directly:

gdalwarp -t_srs '+proj=utm +zone=5 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0' HIDEM/hawaii_dem hawaii_dem_utm.tif

This will likely be the fastest option, and doesn't require setting up a GIS environment explicitly.

share|improve this answer
    
That didn't do it, but this did: setOptions(chunksize = 1e+04, maxmemory = 1e+06) Time eight minutes, much less than it would take to install and use a real GIS. –  J. Winchester Feb 23 '11 at 0:05
    
@J. Winchester: I've updated my response to include your settings as that's the better approach. The package author, would likely be interested to hear when and why it crashes, and hopefully update the defaults to reflect this. –  scw Feb 23 '11 at 0:19
    
it's a good idea to add (lossless) compression and tiling (defaults to 256x256) to the GeoTIFF from gdalwarp if your target can handle it: -co COMPRESS=LZW -co TILED=YES –  mdsumner Feb 23 '11 at 13:10
    
I let Robert Hijmans know about the case. On a smaller DEM, the default settings are near-optimal, so this is a mystery so far. –  J. Winchester Feb 23 '11 at 19:36
    
great! This also allowed me to export a RasterLayer to a (3GB) netcdf with writeRaster. –  David Nov 15 '13 at 22:21

You can also use the spgrass6 package for the integration between R and grass. The author is Roger Bivand (the author of sp)

This package have many function to completly run grass inside R (or the reverse) and exchange data between R and grass

for more information : http://cran.r-project.org/web/packages/spgrass6/index.html

share|improve this answer

HI,

his is a process that takes just a few seconds in GIS software. My attempt to do it in R uses a >large amount of memory then fails.

You answered your questions, do that in GRASS or GDAL and leave R for other tasks.

share|improve this answer
1  
For completeness: you should look to the spgrass package to run grass in R. –  johanvdw Feb 22 '11 at 21:56
1  
And a third option is using saga gis. There is a module (RSAGA) which connects saga and R. –  johanvdw Feb 22 '11 at 22:02
    
This R function is designed to use GDAL, but seems not to be using it well in this case. My question is "How can I best accomplish this task with R", not "What GIS software is available that can do this task." –  J. Winchester Feb 23 '11 at 3:12

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.