无法更改栅格的范围
Can't change raster's extent
我想裁剪高程栅格以将其添加到栅格堆栈。这很简单，我之前很顺利地做到了这一点，将生态区栅格添加到同一个堆栈中。但是对于海拔高度，这是行不通的。现在，溢出中有几个问题解决了这个问题，我尝试了很多东西...
首先，我们需要这个：
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 形状文件：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
但是高程的范围与预测变量 2 的范围略有不同：
> extent(elevation)
class : Extent
xmin : -120.0001
xmax : -35.00014
ymin : -60.00014
ymax : 34.99986
>
我试图通过各种方式使 then 相等 我在此处阅读的问题...
我试图扩展所以海拔的 ymax 会满足 predictors2 的 ymax
elevation<-extend(elevation,predictors2) #无效，范围保持不变
我尝试了相反的方法...使 predictors2 范围符合高程范围...也没有。
但后来我读到
You might not want to play with setExtent() or extent() <- extent(), as you could end with wrong geographic coordinates of your rasters - @ztl, Jun 29 '15
并且我试图通过这样做
获得我的栅格的最小公共范围，按照@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)
我认为您遇到问题是因为您使用的是不同空间分辨率的栅格。因此，当您将两个栅格裁剪到相同范围时，它们的 actual 范围会因此略有不同。
所以如果你想堆叠光栅，你需要让它们进入相同的分辨率。要么用较粗的分辨率分解栅格（即通过重采样或其他方法增加分辨率），要么用更高分辨率聚合栅格（即降低分辨率，例如取 n 像素的平均值）。
请注意，如果您使用 setExtent(x)、extent(x) <-、res(x) <- 或类似的方式更改范围或分辨率，则 不会 起作用，因为您只是在更改栅格对象中的插槽，而不是实际的基础数据。
因此，要使栅格具有共同的分辨率，您需要更改数据。为此，您可以使用函数（以及其他函数）aggregate、disaggregate 和 resample。但是由于您正在更改数据，因此您需要清楚自己在做什么以及您使用的功能在做什么。
对您来说最方便的方法应该是 resample，您可以在其中将一个栅格重新采样为另一个栅格，以便它们在范围和分辨率上匹配。这将使用定义的方法完成。默认情况下，它使用 nearest neighbor 来计算新值。如果您正在处理诸如高程之类的连续数据，您可能希望选择 bilinear，这是双线性插值。在这种情况下，您实际上是在创建 "new measurements"，需要注意的事项。
如果你的两个分辨率是彼此的倍数，你可以查看 aggregate 和 disaggregate。在 disaggregate 的情况下，您可以将光栅单元按一个因子拆分以获得更高的分辨率（例如，如果您的第一个分辨率是 10 度，而您想要的分辨率是 0.05 度，您可以 disaggregate 使用一个因子200 为每 10 度单元提供 200 个 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
我想裁剪高程栅格以将其添加到栅格堆栈。这很简单，我之前很顺利地做到了这一点，将生态区栅格添加到同一个堆栈中。但是对于海拔高度，这是行不通的。现在，溢出中有几个问题解决了这个问题，我尝试了很多东西...
首先，我们需要这个：
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 形状文件：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
但是高程的范围与预测变量 2 的范围略有不同：
> extent(elevation)
class : Extent
xmin : -120.0001
xmax : -35.00014
ymin : -60.00014
ymax : 34.99986
>
我试图通过各种方式使 then 相等 我在此处阅读的问题... 我试图扩展所以海拔的 ymax 会满足 predictors2 的 ymax elevation<-extend(elevation,predictors2) #无效，范围保持不变
我尝试了相反的方法...使 predictors2 范围符合高程范围...也没有。
但后来我读到
You might not want to play with setExtent() or extent() <- extent(), as you could end with wrong geographic coordinates of your rasters - @ztl, Jun 29 '15
并且我试图通过这样做
获得我的栅格的最小公共范围，按照@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)
我认为您遇到问题是因为您使用的是不同空间分辨率的栅格。因此，当您将两个栅格裁剪到相同范围时，它们的 actual 范围会因此略有不同。
所以如果你想堆叠光栅，你需要让它们进入相同的分辨率。要么用较粗的分辨率分解栅格（即通过重采样或其他方法增加分辨率），要么用更高分辨率聚合栅格（即降低分辨率，例如取 n 像素的平均值）。
请注意，如果您使用 setExtent(x)、extent(x) <-、res(x) <- 或类似的方式更改范围或分辨率，则 不会 起作用，因为您只是在更改栅格对象中的插槽，而不是实际的基础数据。
因此，要使栅格具有共同的分辨率，您需要更改数据。为此，您可以使用函数（以及其他函数）aggregate、disaggregate 和 resample。但是由于您正在更改数据，因此您需要清楚自己在做什么以及您使用的功能在做什么。
对您来说最方便的方法应该是 resample，您可以在其中将一个栅格重新采样为另一个栅格，以便它们在范围和分辨率上匹配。这将使用定义的方法完成。默认情况下，它使用 nearest neighbor 来计算新值。如果您正在处理诸如高程之类的连续数据，您可能希望选择 bilinear，这是双线性插值。在这种情况下，您实际上是在创建 "new measurements"，需要注意的事项。
如果你的两个分辨率是彼此的倍数，你可以查看 aggregate 和 disaggregate。在 disaggregate 的情况下，您可以将光栅单元按一个因子拆分以获得更高的分辨率（例如，如果您的第一个分辨率是 10 度，而您想要的分辨率是 0.05 度，您可以 disaggregate 使用一个因子200 为每 10 度单元提供 200 个 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