projectRaster 在应用于 R 中的列表对象时无法更改 crs

projectRaster fails to change crs when applied to a list object in R

我想将 6 个栅格堆叠在一个名为 allrasters 的列表中,但首先必须修复 crs 和范围不一致的问题。这是我的代码尝试将列表中的第二个栅格设置为列表中第三个栅格的 crs:

projectRaster(allrasters[[2]], crs=crs(allrasters[[3]]))

然而,当我 运行 这段代码并检查时,allrasters[[2]] 仍然是 proj.merc 并且没有任何改变...

栅格信息:

crs(allrasters[[2]])
CRS arguments:
 +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0
+x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext
+no_defs

crs(allrasters[[3]])
CRS arguments:
 +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5
+x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 

我假设你想要的是:

allrasters[[2]] <- projectRaster(allrasters[[2]], crs=crs(allrasters[[3]]))

也就是你忘了给projectRaster

的输出赋值

我认为你需要几个步骤:

  1. 您需要在同一投影中获取所有栅格
  2. 您需要找到所有栅格的完整范围,就好像它们镶嵌在一起一样
  3. 您需要对栅格重新采样,以便它们具有相同的范围、分辨率和投影
  4. 您将堆叠光栅。

这是我用一些虚假数据创建的示例,应该可以帮助您完成此操作:


##Loading necessary packages##
library(raster)
library(rgeos)
library(tmaptools)

#For reproducibility#
set.seed(52)

##Creating fake rasters with different extents and projections##
R1<-raster(nrow=100, ncol=100, xmn=44.52, xmx=45.1, ymn=-122.1, ymx=-121.2, crs=crs("+init=epsg:4267"))
R2<-raster(nrow=100, ncol=100, xmn=44.49, xmx=45.8, ymn=-122.0, ymx=-121.3, crs=crs("+init=epsg:4326"))
R3<-raster(nrow=100, ncol=100, xmn=44.48, xmx=45.1, ymn=-122.5, ymx=-121.5, crs=crs("+init=epsg:4979"))
R4<-raster(nrow=100, ncol=100, xmn=44.55, xmx=45.6, ymn=-122.2, ymx=-121.0, crs=crs("+init=epsg:4269"))

values(R1)<-rnorm(10000, 500, 10)
values(R2)<-rnorm(10000, 1000, 60)
values(R3)<-rnorm(10000, 300, 10)
values(R4)<-rnorm(10000, 2500, 70)

##Creating a list of the rasters##
tmp<-list(R1,R2,R3,R4)

##Looping to reproject the rasters all into the same projection##
allras<-list()
for (i in 1:length(tmp)){
  if(i==1){
    allras[[i]]<-tmp[[i]]
  }else{
    allras[[i]]<-projectRaster(tmp[[i]], crs=crs(tmp[[1]]))
  }
  
}

##Creating a function to make a polygon of each raster's extent##
fxn<-function(ras){
  bb<-bbox(ras)
  bbpoly<-bb_poly(bb)
  st_crs(bbpoly)<-crs(ras)
  return(as_Spatial(bbpoly))
}
ext<-lapply(allras, fxn)

##Aggregating and dissolving all extents to get the full extent of all rasters##
full.ext<-aggregate(do.call(bind, ext), dissolve=TRUE)

##Creating a blank raster with the full extent, the desired final projection, and the desired resolution##
blank<-raster(ext=extent(full.ext), nrow=allras[[1]]@nrows, ncol=allras[[1]]@ncols, crs=allras[[1]]@crs)

##Resampling all rasters in the list to the desired extent and resolution##
rastostack<-lapply(allras, resample, y=blank)

##Stacking the rasters##
Ras<-stack(rastostack)