R: Faça o download de um DEM grande, altere a projeção e ajuste em menor escala

11

Este é um processo que leva apenas alguns segundos no software GIS. Minha tentativa de fazê-lo no R usa uma grande quantidade de memória e falha. Existe algo errado no meu código ou isso é algo que R não pode fazer? Eu li que R pode trabalhar dentro do Grass, posso usar uma função Grass de dentro do 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)
J. Win.
fonte
Apenas curioso, qual é o pacote que você está usando?
DJQ
@celenius: este pacote é chamadoraster
J. Win.

Respostas:

9

Do meu olhar para a fonte, rasterparece adivinhar se o conjunto de dados se encaixa na memória e, se for o caso, execute a operação na memória, caso contrário, no disco. Você pode forçá-lo a executar o cálculo configurando explicitamente chunksize(células a serem processadas por vez) e maxmemory(número máximo de células a serem lidas na memória):

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

Como alternativa, você pode executar a transformação com o GDAL diretamente:

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

Essa provavelmente será a opção mais rápida e não requer a configuração explícita de um ambiente GIS.

scw
fonte
Não foi o que aconteceu, mas o que aconteceu: setOptions(chunksize = 1e+04, maxmemory = 1e+06)Tempo de oito minutos, muito menos que o necessário para instalar e usar um GIS real.
J. Win.
@J. Winchester: atualizei minha resposta para incluir suas configurações, pois essa é a melhor abordagem. O autor do pacote provavelmente estaria interessado em saber quando e por que falha e, esperançosamente, atualizar os padrões para refletir isso.
scw 23/02
1
é uma boa idéia para adicionar compressão (lossless) e azulejos (o padrão é 256x256) para o GeoTIFF de gdalwarp se o seu alvo pode lidar com isso: -co COMPRESS = LZW -co TELHADO = YES
mdsumner
Eu avisei Robert Hijmans sobre o caso. Em um DEM menor, as configurações padrão são quase ideais, então isso é um mistério até agora.
J. Win.
ótimo! Isso também me permitiu exportar um RasterLayer para um (3GB) netcdf com writeRaster.
David LeBauer
3

Você também pode usar o pacote spgrass6 para a integração entre R e grama. O autor é Roger Bivand (o autor de sp)

Este pacote possui muitas funções para executar a grama completamente dentro de R (ou o inverso) e trocar dados entre R e grama

para obter mais informações: http://cran.r-project.org/web/packages/spgrass6/index.html

dickoa
fonte
1

OI,

Esse é um processo que leva apenas alguns segundos no software GIS. Minha tentativa de fazer isso em R usa uma grande quantidade de memória e falha.

Você respondeu às suas perguntas, faça isso em GRASS ou GDAL e deixe R para outras tarefas.

Pablo
fonte
1
Para completude: você deve procurar o pacote spgrass para executar a grama em R. #
319 johanvdw
1
E uma terceira opção é usar saga gis. Existe um módulo (RSAGA) que conecta a saga e a R. #
3111 johanvdw
Essa função R foi projetada para usar o GDAL, mas parece não estar sendo usada corretamente neste caso. Minha pergunta é "Como melhor posso realizar esta tarefa com R", não "Qual software GIS está disponível para executar essa tarefa".
J. Win.