无法更改栅格范围

Question

我想裁剪高程栅格以将其添加到栅格堆栈。这很容易，我之前很顺利地做到了这一点，将生态区域栅格添加到同一个堆栈中。但是对于海拔一号，它是行不通的。现在，这里有几个问题在溢出解决这个问题，我尝试了很多东西......

首先，我们需要这个：

library(rgdal)
library(raster)

我的堆栈是 predictors2：

#Downloading the stack
predictors2_full<-getData('worldclim', var='bio', res=10)

#Cropping it, I don' need the whole world   
xmin=-120; xmax=-35; ymin=-60; ymax=35
limits <- c(xmin, xmax, ymin, ymax)
predictors2 <- crop(predictors2_full,limits)

然后我在这里下载了 terr_ecorregions shapefile：http://maps.tnc.org/files/shp/terr-ecoregions-TNC.zip

setwd("~/ORCHIDACEAE/Ecologicos/w2/layers/terr-ecoregions-TNC")
ecoreg = readOGR("tnc_terr_ecoregions.shp") # I've loaded...
ecoreg2 <- crop(ecoreg,extent(predictors2)) # cropped...
ecoreg2 <- rasterize(ecoreg2, predictors2)  # made the shapefile a raster
predictors4<-addLayer(predictors2,elevation,ecoreg2) # and added the raster
                                                     # to my stack

有了海拔，我就做不到了。数字高程模型基于 GMTED2010，可在此处下载：http://edcintl.cr.usgs.gov/downloads/sciweb1/shared/topo/downloads/GMTED/Grid_ZipFiles/mn30_grd.zip

elevation<-raster("w001001.adf") #I've loaded
elevation<-crop(elevation,predictors2) # and cropped

但是海拔的范围与 predictors2 的范围略有不同：

> extent(elevation)
class       : Extent 
xmin        : -120.0001 
xmax        : -35.00014 
ymin        : -60.00014 
ymax        : 34.99986 
> 

我试图使 then 平等，我在这里读到的问题... 我试图扩展，所以海拔的 ymax 会满足 predictors2 的 ymax 海拔

我尝试了相反的方法...使 predictors2 范围满足海拔的范围...也没有。

但后来我读到了

您可能不想使用 setExtent() 或 extent()

并且我试图通过这样做来获得我的栅格的最小公共范围，遵循@zlt 在另一个范围问题中的回答

# Summing your rasters will only work where they are not NA
r123 = r1+r2+r3 # r123 has the minimal common extent
r1 = crop(r1, r123) # crop to that minimal extent
r2 = crop(r2, r123)
r3 = crop(r3, r123)

为此，首先我必须设置分辨率：

res(elevation)<-res(predictors2) #fixing the resolutions... This one worked.

但是，r123 = r1+r2+r 不起作用：

> r123=elevation+ecoreg2+predictors2
Error in elevation + ecoreg2 : first Raster object has no values

谁能给我一个提示？我真的很想将我的高程添加到栅格中。有趣的是，我有另一个名为 predictors1 的堆栈，其高度范围完全相同......而且我能够裁剪 ecoreg 并将 ecoreg 添加到 predictors1 和 predictors2... 为什么我不能对海拔做同样的事情？ 我对这个世界还很陌生，没有什么想法……我很感激任何提示。

编辑：解决方案，感谢@Val

我明白了：

#Getting the factor to aggregate (rasters are multiples of each other)
res(ecoreg2)/res(elevation)
[1] 20 20 #The factor is 20


elevation2<-aggregate(elevation, fact=20)
elevation2 <- crop(elevation2,extent(predictors2))

#Finally adding the layer:
predictors2_eco<-addLayer(predictors2,elevation2,ecoreg)

新问题，想...

我无法将堆栈写入 geotiff

 writeRaster(predictors2_eco, filename="cropped_predictors2_eco.tif", options="INTERLEAVE=BAND", overwrite=TRUE)

 Error in .checkLevels(levs[[j]], value[[j]]) : 
 new raster attributes (factor values) should be in a data.frame (inside a list)

Answer 1

我认为您遇到了问题，因为您正在使用不同空间分辨率的栅格。因此，当您将两个栅格裁剪到相同的范围时，它们的 实际 范围会因此而略有不同。 因此，如果要堆叠栅格，则需要使它们具有相同的分辨率。要么用较粗的分辨率分解栅格（即通过重新采样或其他方法增加分辨率），要么用较高的分辨率聚合栅格（即降低分辨率，例如采用n 像素的平均值）。

请注意，如果您使用setExtent(x)、extent(x) <-、res(x) <- 或类似名称更改范围或分辨率，则不会起作用，因为您只是在更改栅格对象中的插槽，不是实际的基础数据。

因此，要使栅格达到共同的分辨率，您需要更改数据。为此，您可以使用函数（以及其他函数）aggregate、disaggregate 和 resample。但是由于您正在更改数据，因此您需要清楚自己是什么以及您使用的功能正在做什么。

对您来说最方便的方法应该是resample，您可以在其中将一个栅格重新采样到另一个栅格，以便它们在范围和分辨率上匹配。这将使用定义的方法完成。默认情况下，它使用nearest neighbor 来计算新值。如果您正在处理诸如高程之类的连续数据，您可能需要选择bilinear，它是双线性插值。在这种情况下，您实际上是在创建“新测量”，这是需要注意的。

如果您的两个分辨率是彼此的倍数，您可以查看aggregate 和disaggregate。在disaggregate 的情况下，您可以将光栅单元拆分一个因子以获得更高的分辨率（例如，如果您的第一个分辨率是 10 度，而您想要的分辨率是 0.05 度，您可以将 disaggregate 用 200 的因子得到 200每 10 度的单元格 0.05 度的单元格）。这种方法可以避免插值。

这里有一个小例子：

library(raster)
library(rgeos)

shp <- getData(country='AUT',level=0)

# get centroid for downloading eco and dem data
centroid <- coordinates(gCentroid(shp))

# download 10 degree tmin
ecovar <- getData('worldclim', var='tmin', res=10, lon=centroid[,1], lat=centroid[,2])
ecovar_crop <- crop(ecovar,shp)

# output
> ecovar_crop
class       : RasterBrick 
dimensions  : 16, 46, 736, 12  (nrow, ncol, ncell, nlayers)
resolution  : 0.1666667, 0.1666667  (x, y)
extent      : 9.5, 17.16667, 46.33333, 49  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : tmin1, tmin2, tmin3, tmin4, tmin5, tmin6, tmin7, tmin8, tmin9, tmin10, tmin11, tmin12 
min values  :  -126,  -125,  -102,   -77,   -33,    -2,    19,    20,     5,    -30,    -74,   -107 
max values  :   -31,   -21,     9,    51,    94,   131,   144,   137,   106,     60,     18,    -17 


# download SRTM elevation - 90m resolution at eqt
elev <- getData('SRTM',lon=centroid[,1], lat=centroid[,2])
elev_crop <- crop(elev, shp)

# output

> elev_crop
class       : RasterLayer 
dimensions  : 3171, 6001, 19029171  (nrow, ncol, ncell)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : 9.999584, 15.00042, 46.37458, 49.01708  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : srtm_39_03 
values      : 198, 3865  (min, max)

# won't work because of different resolutions (stack is equal to addLayer)
ecoelev <- stack(ecovar_crop,elev_crop)

# resample
elev_crop_RS <- resample(elev_crop,ecovar_crop,method = 'bilinear')

# works now
ecoelev <- stack(ecovar_crop,elev_crop_RS)

# output
> ecoelev
class       : RasterStack 
dimensions  : 16, 46, 736, 13  (nrow, ncol, ncell, nlayers)
resolution  : 0.1666667, 0.1666667  (x, y)
extent      : 9.5, 17.16667, 46.33333, 49  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       :     tmin1,     tmin2,     tmin3,     tmin4,     tmin5,     tmin6,     tmin7,     tmin8,     tmin9,    tmin10,    tmin11,    tmin12, srtm_39_03 
min values  : -126.0000, -125.0000, -102.0000,  -77.0000,  -33.0000,   -2.0000,   19.0000,   20.0000,    5.0000,  -30.0000,  -74.0000, -107.0000,   311.7438 
max values  :   -31.000,   -21.000,     9.000,    51.000,    94.000,   131.000,   144.000,   137.000,   106.000,    60.000,    18.000,   -17.000,   3006.011