R: Téléchargez un grand DEM, modifiez la projection et ajustez à une échelle plus petite

Il s'agit d'un processus qui ne prend que quelques secondes dans le logiciel SIG. Ma tentative de le faire dans R utilise une grande quantité de mémoire puis échoue. Y a-t-il quelque chose de mal dans mon code, ou est-ce juste quelque chose que R ne peut pas faire? J'ai lu que R peut fonctionner dans Grass, puis-je utiliser une fonction Grass depuis l'intérieur de 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)

raster grass r J. Win.
la source

Juste curieux, quel est le package que vous utilisez?

djq

@celenius: ce paquet s'appelleraster

J. Win.

Réponses:

De mon regard sur la source, rastercherche à deviner si l'ensemble de données s'inscrit dans la mémoire, et si oui, effectuez l'opération en mémoire, sinon sur le disque. Vous pouvez le forcer à effectuer le calcul en définissant explicitement chunksize(cellules à traiter à la fois) et maxmemory(nombre maximal de cellules à lire en mémoire):

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

Alternativement, vous pouvez effectuer la transformation directement avec GDAL:

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

Ce sera probablement l'option la plus rapide et ne nécessite pas la configuration explicite d'un environnement SIG.

scw
la source

Cela ne l'a pas fait, mais cela l'a fait: setOptions(chunksize = 1e+04, maxmemory = 1e+06)Temps de huit minutes, beaucoup moins que ce qu'il faudrait pour installer et utiliser un vrai SIG.

J. Win.

@J. Winchester: J'ai mis à jour ma réponse pour inclure vos paramètres car c'est la meilleure approche. L'auteur du package serait probablement intéressé de savoir quand et pourquoi il se bloque et, espérons-le, de mettre à jour les valeurs par défaut pour refléter cela.

scw

c'est une bonne idée d'ajouter la compression (sans perte) et la mosaïque (par défaut 256x256) au GeoTIFF de gdalwarp si votre cible peut le gérer: -co COMPRESS = LZW -co TILED = YES

mdsumner

J'ai informé Robert Hijmans de l'affaire. Sur un DEM plus petit, les paramètres par défaut sont presque optimaux, c'est donc un mystère jusqu'à présent.

J. Win.

génial! Cela m'a également permis d'exporter un RasterLayer vers un netcdf (3 Go) avec writeRaster.

David LeBauer

Vous pouvez également utiliser le package spgrass6 pour l'intégration entre R et grass. L'auteur est Roger Bivand (l'auteur de sp)

Ce paquet a de nombreuses fonctions pour exécuter complètement l'herbe à l'intérieur de R (ou l'inverse) et échanger des données entre R et l'herbe

pour plus d'informations: http://cran.r-project.org/web/packages/spgrass6/index.html

Dickoa
la source

SALUT,

C'est un processus qui ne prend que quelques secondes dans un logiciel SIG. Ma tentative de le faire dans R utilise une> grande quantité de mémoire puis échoue.

Vous avez répondu à vos questions, faites-le dans GRASS ou GDAL et laissez R pour d'autres tâches.

Pablo
la source

Pour être complet: vous devriez regarder le paquet spgrass pour faire pousser de l'herbe dans R.

johanvdw

Et une troisième option consiste à utiliser saga gis. Il existe un module (RSAGA) qui relie la saga et R.

johanvdw

Cette fonction R est conçue pour utiliser GDAL, mais ne semble pas bien l'utiliser dans ce cas. Ma question est "Comment puis-je accomplir au mieux cette tâche avec R", et non "Quel logiciel SIG est disponible pour effectuer cette tâche."

J. Win.