代码之家 › 专栏 › 技术社区 › Neal Barsch

在R中仅读取光栅的一个裁剪或范围

geotiff tiff raster geospatial r

Neal Barsch · 技术社区 · 7 年前

我正在处理每日TIF文件(数千份每日文件)中的大量数据。我正在分析一个shapefile区域中光栅的平均值,该区域可能重复了数千个形状层。

我目前的.tif文件适用于整个国家,而实际上我只需要每个shapefile层的每个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))

然后,通常的Velox(或光栅)提取,重复。

2 回复 | 直到 7 年前

Val Srinivas 7 年前

如果您对数据感兴趣而不是光栅对象(或VRT/TIFF),那么这可能是一种方法:

如果您从文件“加载”光栅,并不一定意味着它将作为文档加载到内存中。 raster 解释:

在许多情况下,例如,当从文件创建RasterLayer时,它(最初)在(RAM)内存中不包含任何单元(像素)值,它只有描述RasterLayer的参数。可以使用getValues、extract和相关函数访问单元格值。可以使用setvalues和replacement指定新值。

因此,您可以先“索引”光栅,然后使用 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

loki 7 年前

在这种情况下,我会创建一个 virtual raster file (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的文件大小不应超过某些KB。