Только чтение урожая или экстента растра в R - PullRequest
Новая
       108

Только чтение урожая или экстента растра в R

0 голосов
/ 13 июня 2018

Я работаю с огромным количеством данных в ежедневных tif-файлах (тысячи ежедневных файлов).Я анализирую среднее значение растра в области шейп-файла, которое повторяется на тысячах слоев формы.

Мои текущие TIF-файлы предназначены для всей страны, хотя на самом деле мне нужна только небольшая площадь каждого TIF-файла для каждого слоя шейп-файла.Каждый слой формы имеет различный набор дней для извлечения, т. Е. Фрейм данных, подобный этому: 

df <- data.frame(shplayer=c("shp_layerbuffer1","shp_layerbuffer2", "shp_layerbuffer3"), start=c("2000_02_18", "2004_03_19", "2010_08_20"), end=c("2010_01_09", "2005_03_15", "2012_09_04"))

Есть ли в R способ обрезать .tif (или любой файл растрового типа) ДОчитаете файл?Таким образом, я мог прочитать только обрезанную область каждого из файлов tif

Я предусмотрел что-то вроде (повторяя по всему набору дат): 

library(sf)
library(raster)
shp_layerbuffer1 <- st_read("myshpfolder", layer="shp_layerbuffer1", quiet=T)

###EXAMPLE BUT doesn't work to crop the raster as it comes in
tempraster <- raster::raster("mytif_2000_02_18.tif", ext=extent(shp_layerbuffer1))

Затем обычная скорость (илирастр) извлечь, повторить.

Ответы [ 2 ]

0 голосов
/ 14 июня 2018

Если вас интересуют данные, а не растровый объект (или vrt / tiff), то это может быть подходом:

Если вы «загружаете» растр из файла, он необязательно означает, что он загружается в память, как объясняется в документации raster:

Во многих случаях, например, когда RasterLayer создается из файла, он (изначально) не содержит никаких ячеек(пиксельные) значения в (RAM) памяти, он имеет только те параметры, которые описывают RasterLayer.Вы можете получить доступ к значениям ячеек с помощью getValues, extract и связанных функций.Вы можете назначать новые значения с помощью setValues ​​и с заменой.

Таким образом, вы можете сначала «проиндексировать» растр, а затем использовать getValuesBlock для считывания значений, попадающих в интересующую вас область. 

library(raster)

f <- system.file('external/test.grd',package = 'raster')

# index
r <- raster(f)

# check if in memory

inMemory(r)
#[1] FALSE # output

# this would be an extent from your overlapping shapefile
e <- extent(r,58,68,40,50)

# get cells from extent; either use cells as index directly or convert to rowcol

rowcol <- rowColFromCell(r,cellsFromExtent(r,e))

v <- getValuesBlock(r,row=rowcol[1,1],nrows=(rowcol[nrow(rowcol),1] - rowcol[1,1]),
                    col=rowcol[1,2],ncols=(rowcol[nrow(rowcol),2] - rowcol[1,2]))

v
# [1] 295.392 225.710 258.595 310.336 324.666 322.702 307.217 283.611 263.987 156.609 322.978 297.565 301.797 315.971
# [15] 323.605 326.920 317.894 297.138 270.027 225.769 337.503 323.388 314.900 308.877 314.556 343.555 344.035 315.400
# [29] 291.188 270.876 337.866 324.632 307.220 278.294 264.379 348.519 356.456 322.450 301.790 285.815 331.383 318.950
# [43] 299.246 262.026 230.869 294.012 320.274 312.777 300.513 285.317 321.075 311.447 296.952 278.519 270.279 283.797
# [57] 296.415 298.861 293.150 277.573 306.692 300.772 289.376 273.141 258.457 258.638 272.966 283.977 284.621 271.690
# [71] 286.749 286.617 275.618 247.307 198.884 193.865 240.465 268.687 277.303 271.431 260.336 271.458 263.977 231.071
# [85] 161.407 154.181 222.417 258.672 271.711 272.642 235.458 258.810 257.763 239.553 209.409 205.905 234.162 255.668
# [99] 266.260 270.532

РЕДАКТИРОВАТЬ: 

# check if raster was loaded after getValuesBlock
inMemory(r)
#[1] FALSE
0 голосов
/ 13 июня 2018

В таком случае я бы создал виртуальный растровый файл (VRT) с пакетом gdalUitls.Это самый быстрый способ создания подмножества (или даже виртуальной мозаики) из одного или нескольких наборов растровых данных.Для этого VRT вы можете установить желаемый экстент , используя sf::extent.

Минимальный (не воспроизводимый) пример: 

library(gdalUitls)
library(raster)
library(sf)

shp_layerbuffer1 <- st_read("myshpfolder", layer="shp_layerbuffer1", quiet=T)

gdalbuildvrt(gdalfile = "yourraster.tif", 
             output.vrt = "tmp.vrt", 
             te = st_bbox(shp_layerbuffer1))

tempraster <- raster("tmp.vrt")

VRT в основном создает виртуальный холст для нового "растрового" файла.Затем он ссылается на соответствующие строки и столбцы из (одного или нескольких) базовых растровых файлов.Следовательно, вам не нужно создавать целый новый растр, а просто ссылаться на существующие файлы.Размер файла вашего VRT не должен превышать несколько килобайт.

...