Читать только обрезку или экстент растра в R

Я работаю с огромным количеством данных в ежедневных файлах 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))

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


r geospatial raster tiff geotiff
person Neal Barsch    schedule 13.06.2018    source источник

Ответы (2)


arrow_upward
3
arrow_downward

Если вас интересуют данные, а не растровый объект (или 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
person Val    schedule 14.06.2018
comment
@loki Нет, не верю. Вызов raster не загружает данные, только метаданные. Только getValues или что-то подобное загрузит данные в память ... вот тут-то и пригодится подмножество. - person Val; 14.06.2018
comment
Не волнуйтесь ... Лично мне очень нравится vrt ... может быть, немного переборщить для этой цели - но на самом деле все сводится к личному вкусу, ИМО :) - person Val; 14.06.2018
comment
Вы абсолютно правы. В последнее время я научился любить vrts и довольно часто ими пользуюсь - person loki; 14.06.2018
comment
Я не уверен, что vrts быстрее, чем чтение прямых tif-файлов в velox, по крайней мере, для моих целей извлечения растров. Смотрите мои тесты здесь: stackoverflow.com/questions/50870080/ - person Neal Barsch; 15.06.2018
comment
@NealBarsch Ваш первоначальный вопрос заключался в том, как частично читать файлы (также известные как экстенты кадрирования) в память. Оба решения ответили на него, и если вам понравился ответ Локи, я не вижу причин, чтобы не принимать его. Аспект производительности, кажется, рассматривается в вашем последующем вопрос - person Val; 15.06.2018
comment
В основном я пытался спровоцировать еще несколько споров о том, что является самым быстрым :) - person Neal Barsch; 15.06.2018

arrow_upward
4
arrow_downward

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

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

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 не должен превышать нескольких килобайт.

person loki    schedule 13.06.2018
comment
Привет, Локи, большое спасибо за это! Я продолжаю получать Error in paste(parameter_values[[which(names(parameter_values) == X)]], : non-string argument to internal 'paste' какие-нибудь идеи, почему? - person Neal Barsch; 14.06.2018
comment
Так что проверьте это: если вы кадрируете один раз в растре, это происходит намного быстрее, но не при загрузке файла. Также спасибо за вашу работу над этим! - person Neal Barsch; 16.06.2018