0

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

首先，我们需要这个：

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 海拔<-extend(elevation,predictors2) #没有用，范围保持不变

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

但后来我读到了

您可能不想使用 setExtent() 或 extent() <- extent()，因为您可能以错误的栅格地理坐标结束 - @ztl，2015 年 6 月 29 日

我试图通过这样做来获得我的栅格的最小共同范围，遵循@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)
4

1 回答 1

1

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

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

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

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

如果您的两个分辨率是彼此的倍数，您可以查看aggregateand disaggregate。如果disaggregate您将光栅单元拆分为一个因子以获得更高的分辨率（例如，如果您的第一个分辨率是 10 度，而您想要的分辨率是 0.05 度，那么您可以disaggregate使用 200 的因子为您提供 200 个 0.05 度的单元10 度电池）。这种方法将避免插值。

这是一个工作示例：

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
于 2017-08-12T10:20:55.993 回答