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

Posted 2023-03-23

【中文标题】应用于 R 中的列表对象时，projectRaster 无法更改 crs

【英文标题】：projectRaster fails to change crs when applied to a list object in R

【发布时间】：2022-01-18 06:02:36

【问题描述】：

我想在一个名为 allrasters 的列表中堆叠 6 个栅格，但首先必须修复 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 

【参考方案1】：

我假设你所追求的是：

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

也就是你忘了分配projectRaster的输出

【参考方案2】：

我认为您需要几个步骤：

您需要在同一投影中获取所有栅格 您需要找到所有栅格的完整范围，就像它们被镶嵌在一起一样 您需要对栅格重新采样，以使它们具有相同的范围、分辨率和投影 您将堆叠您的栅格。

这是我使用一些虚假数据创建的示例，可以帮助您完成此操作：

##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)