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)
Answer
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.
No comments:
Post a Comment