无法更改栅格范围 [英] Can't change raster's extent

问题描述

我想裁剪高程栅格以将其添加到栅格堆栈中.这很容易，我之前做得很顺利，在同一堆栈中添加了一个生态区域栅格.但是随着海拔的升高，那是行不通的.现在，这里有几个问题可以解决这个问题，我尝试了很多事情...

I want to crop an elevation raster to add it to a raster stack. It's easy, I did this before smoothly, adding a ecoregions raster to the same stack. But with the elevation one, just doesn't work. Now, there are several questions here in overflow adressing this issue and I tryed a lot of things...

首先，我们需要这个:

library(rgdal)
library(raster)

我的堆栈是预测变量2:

My stack is 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

Then I've downloaded the terr_ecorregions shapefile here: 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

With elevation, I just can't. The Digital elevation model is based in GMTED2010, which can be downloaded here: 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的程度:

But elevation gets a slightly different extent instead of predictors2's extent:

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

我想尽一切办法使我在这里的问题中读到相等... 我试图扩展以使高程的ymax满足预测变量2的ymax 海拔< -extend(elevation，predictors2)#无效，范围保持不变

I tried to make then equal by all means I read about in questions here... I tried to extend so elevation's ymax would meet predictors2's ymax elevation<-extend(elevation,predictors2) #didn't work, extent remains the same

我尝试了相反的做法...使预报器2的范围满足海拔的范围...也一无所获.

I tried the opposite... making predictors2 extent meet elevation's extent... nothing either.

但是我读到了

您可能不希望使用setExtent()或scope()<-range()，因为可能会以错误的栅格地理坐标结尾-@ ztl，2015年6月29日

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回答之后，获取栅格的最小公共范围

And I tried to get the minimal common extent of my rasters, following @zlt answer in another extent question, by doing this

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

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

For that, first I had to set the resolutions:

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" ...为什么我不能对海拔高度做同样的事情? 我对这个世界还很陌生，没有什么主意...感谢任何提示.

Can anyone give me a hint on this? I really would like to add my elevation to the raster. Funny thing is, I have another stack named predictors1 with the exact same elevation's extent... And I was able to crop ecoreg and add ecoreg to both predictors1 and predictors2... Why can't I just do the same to elevation? I'm quite new to this world and runned out of ideas... I appreciate any tips.

我明白了:

#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

New problem, thought...

I can't write stack to a 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像素的均值来降低分辨率).

I think you're having issues because you're working with rasters of different spatial resolutions. So when you crop both rasters to the same extent, they'll have a slightly different actual extent because of that. So if you want to stack rasters, you need to get them into the same resolution. Either you disaggregate the raster with the coarser resolution (i.e. increase the resolution by resampling or other methods) or you aggregate the raster with the higher resolution (i.e. decrease the resolution with for instance taking the mean over n pixel).

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

Please note that if you change the extent or resolution with setExtent(x), extent(x) <-, res(x) <- or similar will NOT work, since you're just changing slots in the raster object, not the actual underlying data.

因此，为了使栅格具有共同的分辨率，您需要更改数据.为此，您可以使用功能aggregate，disaggregate和resample(以及其他功能).但是，由于您要更改数据，因此需要清楚自己所处的位置以及所使用的功能.

So to bring the rasters into a common resolution, you need to change the data. You can use the functions (amongst others) aggregate, disaggregate and resample for that purpose. But since you're changing data, you need to be clear on what you're and the function you use is doing.

最方便的方法应该是resample，在这里您可以将一个栅格重新采样到另一个栅格，以便它们在范围和分辨率上匹配.这将使用定义的方法来完成.默认情况下，它使用nearest neighbor来计算新值.如果要处理诸如海拔等连续数据，则可能要选择bilinear，它是双线性插值.在这种情况下，您实际上是在创建新度量"，这是需要注意的.

The most handy way for you should be resample, where you can resample a raster to another raster so they match in extent and resolution. This will be done using a defined method. Per default it's using nearest neighbor for computing the new values. If you're working with continuous data such as elevation, you might want to opt for bilinear which is bilinear interpolation. In this case you're actually creating "new measurements", something to be aware of.

如果两个分辨率互为倍数，则可以查看aggregate和disaggregate.对于disaggregate，您可以将栅格像元除以一个因子以获得更高的分辨率(例如，如果您的第一个分辨率为10度，而所需的分辨率为0.05度，则您可以disaggregate使用200的因数得到200每10度像元0.05度的像元).这种方法可以避免插值.

If your two resolutions are multiples of each other, you could look into aggregate and disaggregate. In the case of disaggregate you would split a rastercell by a factor to get a higher resolution (e.g. if your first resolution is 10 degrees and your desired resolution is 0.05 degrees, you could disaggregate with a factor of 200 giving you 200 cells of 0.05 degree for every 10 degree cell). This method would avoid interpolation.

这是一个可行的例子:

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 

这篇关于无法更改栅格范围的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！