如何用不同程度的栅格创建栅格砖？

Question

我是 R 的新手，所以这个问题非常基础，但我一直在努力解决，找不到有效的解决方案。我想从同一区域的一些陆地卫星图像创建一个光栅砖。 它们以HDF-EOS 格式下载，我使用Modis Reprojection Tool 将它们转换为.tif。

生成的栅格具有相同的投影，但它们的范围、分辨率和原点不同。

我尝试了几种方法，总结如下：

1. 手动定义子集范围并对所有栅格进行子集化。然后尝试用子集栅格制作砖块

2. 重新采样栅格，为它们提供相同的列数和行数。理想情况下，这将确保光栅单元对齐并且可以放入光栅砖中。此选项创建了一个砖块，其中栅格没有值，它们是空的。

我想知道我应该遵循什么概念来纠正程度。创建一个空栅格，然后用导入的 landsat 图像的值填充它是否正确（且有效）？你能看出我在哪里犯错了吗？ 如果相关，我正在使用 Mac OSX 版本 10.9.1，并使用 rgdal 版本 0.8-14

任何帮助将不胜感激！

谢谢

我在这里添加我一直在使用的代码：

# .tif files have been creating using the Modis Reprojection Tool. Input
# files used for this Tool was LANDSAT HDF-EOS imagery.

library(raster)
library(rgdal)

setwd()=getwd()

# Download the files from dropbox:
dl_from_dropbox <- function(x, key) {
  require(RCurl)
  bin <- getBinaryURL(paste0("https://dl.dropboxusercontent.com/s/", key, "/", x),
                      ssl.verifypeer = FALSE)
  con <- file(x, open = "wb")
  writeBin(bin, con)
  close(con)
  message(noquote(paste(x, "read into", getwd())))
}
dl_from_dropbox("lndsr.LT52210611985245CUB00-vi.NDVI.tif", "qb1bap9rghwivwy")
dl_from_dropbox("lndsr.LT52210611985309CUB00-vi.NDVI.tif", "sbhcffotirwnnc6")
dl_from_dropbox("lndsr.LT52210611987283CUB00-vi.NDVI.tif", "2zrkoo00ngigfzm")



# Create three rasters
tif1 <- "lndsr.LT52210611985245CUB00-vi.NDVI.tif"
tif2 <- "lndsr.LT52210611985309CUB00-vi.NDVI.tif"
tif3 <- "lndsr.LT52210611987283CUB00-vi.NDVI.tif"

r1 <- raster(tif1, values=TRUE)
r2 <- raster(tif2, band=1, values=TRUE)
r3 <- raster(tif3, band=1, values=TRUE)

### Display their properties
# projection is identical for the three rasters
projection(r1)
# "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
projection(r2)
projection(r3)

# Extents are different
extent(r1)
# class       : Extent 
# xmin        : -45.85728 
# xmax        : -43.76855 
# ymin        : -2.388705 
# ymax        : -0.5181549
extent(r2)
# class       : Extent 
# xmin        : -45.87077 
# xmax        : -43.78204 
# ymin        : -2.388727 
# ymax        : -0.5208711 
extent(r3)
# class       : Extent 
# xmin        : -45.81952 
# xmax        : -43.7173 
# ymin        : -2.405129 
# ymax        : -0.5154312

# origin differs for all
origin(r1)
# 5.644590e-05 -8.588605e-05
origin(r2)
# 0.0001122091 -0.0001045107
origin(r3)
# 6.949976e-05 -5.895945e-05

# resolution differs for r2
res(r1)
# 0.0002696872 0.0002696872
res(r2)
# 0.0002696875 0.0002696875
res(r3)
# 0.0002696872 0.0002696872


## Try different approaches to create a raster brick

# a- define a subset extent, and subset all the rasters
plot(r1, main="layer1 NDVI")
de <- drawExtent(show=TRUE, col="red")
de
# class       : Extent 
# xmin        : -45.36159 
# xmax        : -45.30108 
# ymin        : -2.002435 
# ymax        : -1.949501
e <- extent(-45.36159,-45.30108,-2.002435,-1.949501)
# Crop each raster with this extent
r1c <- crop(r1,e)
r2c <- crop(r2,e)
r3c <- crop(r3,e)
# Make raster brick
rb_a <- brick(r1c,r2c,r3c)
# Error in compareRaster(x) : different extent

# b- Resample each raster
s <- raster(nrow=6926, ncol=7735)  # smallest nrow and ncol among r1,r2 and r3
r1_res <- resample(r1,s, method="ngb")
r2_res <- resample(r2,s, method="ngb")
r3_res <- resample(r3,s, method="ngb")
# Resampling gives for the three rasters the following message:
# Warning message:
#   In .local(x, y, ...) :
#   you are resampling y a raster with a much larger cell size, 
#   perhaps you should use "aggregate" first

# Make raster brick
rb_c <- brick(r1, r2, r3)
# Error in compareRaster(x) : different extent

Answer 1

这里有一些东西可以帮助你。因为我没有你的 .tif 文件，所以只是一些提示。您是否检查了您的栅格 s 上的范围？这是世界的大小，只有那些列，它的单元格非常大。因此，您必须在重新采样之前向栅格添加范围。根据您的信息，我做了这样的事情：

# create an extent that includes all your data
e<-extent(-46, -43, -2, -0.6)

# create a raster with that extent, and the number of rows and colums to achive a
# similar resolution as you had before, you might have to do some math here....
# as crs, use the same crs as in your rasters before, from the crs slot
s<-raster(e, nrows=7000, ncols=7800, crs=r1@crs)

# use this raster to reproject your original raster (since your using the same crs,
# resample should work fine
r1<-resample(r1, s, method="ngb")

节日快乐， 本

PS 一种更好的方式来找到您想要的范围和分辨率：

# dummy extent from your rasters, instead use lapply(raster list, extent)
a<-extent(-45.85728, -43.76855, -2.388705, -0.5181549)
b<-extent(-45.87077, -43.78204, -2.388727, -0.5208711) 
c<-extent(-45.81952 ,-43.7173 , -2.405129 ,-0.5154312)
extent_list<-list(a, b, c)

# make a matrix out of it, each column represents a raster, rows the values
extent_list<-lapply(extent_list, as.matrix)
matrix_extent<-matrix(unlist(extent_list), ncol=length(extent_list)
rownames(matrix_extent)<-c("xmin", "ymin", "xmax", "ymax")

# create an extent with the extrem values of your extent
best_extent<-extent(min(matrix_extent[1,]), max(matrix_extent[3,]),
min(matrix_extent[2,]), max(matrix_extent[4,]))

# the range of your extent in degrees
ranges<-apply(as.matrix(best_extent), 1, diff)
# the resolution of your raster (pick one) or add a desired resolution
reso<-res(r1)
# deviding the range by your desired resolution gives you the number of rows and columns
nrow_ncol<-ranges/reso

# create your raster with the following
s<-raster(best_extent, nrows=nrow_ncol[2], ncols=nrow_ncol[1], crs=r1@crs)