Friday 21 September 2018

raster - R: Download a large DEM, change projection, and adjust to smaller scale


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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...