用R中的多边形裁剪栅格：错误范围不重合

Question

我希望使用我在ArcGIS中制作的多边形矢量文件来裁剪一个栅格堆栈，然而我收到了“范围不重叠”的错误提示。

首先，我创建了栅格堆栈：

test1 < stack("C:/mydir/test1.tif")

定义投影

myCRS <- test1@crs

然后阅读shapefile。

myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)

裁剪

myCrop <- crop(test1, myExtent)
Error in .local(x, y, ...) : extents do not overlap

我已经寻找解决方案，但我只发现投影可能是问题，然而它们肯定都在相同的CRS中：

> test1$test1.1
class       : RasterLayer 
band        : 1  (of  4  bands)
dimensions  : 10980, 10980, 120560400  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 6e+05, 709800, 5690220, 5800020  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +ellps=WGS84 
+towgs84=0,0,0 
data source : C:\mydir\test1.tif 
names       : test1.1 
values      : 0, 65535  (min, max)

> myExtent
class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 499386.6, 517068.2, 6840730, 6857271  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +ellps=WGS84 
+towgs84=0,0,0 
variables   : 2
names       :    Shape_Leng,    Shape_Area 
min values  : 67444.6461177, 283926851.657 
max values  : 67444.6461177, 283926851.657 

- dtanon

Answer 1

信息已经很清楚了。您的范围不重叠...以下是如何检查它的方法：

library(raster)
ext.ras <- extent(6e+05, 709800, 5690220, 5800020)
ext.pol <- extent(499386.6, 517068.2, 6840730, 6857271)


plot(ext.ras, xlim = c( 499386.6,709800), ylim= c(5690220,6857271), col="red")
plot(ext.pol, add=T, col="blue")

我已经从您的问题中创建了extent对象。您可以在左上角看到一个extent，另一个在右下角。您是否尝试在QGIS中读取两个文件？这样您可能很容易就能看到它们。

如果它们确实应该重叠，那么我会怀疑您读取形状文件的方式。不要使用

，而是使用

myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)

使用：

library(rgdal)
myExtent <- readOGR("C:/mydir","loc1.shp")
myExtent <- spTRansform(myExtent, CRS(proj4string(test1)))

- Bastien

---

谢谢，我按照你的建议使用了readOGR()和spTransform()，效果非常完美。 - dtanon

这可能是因为在你的调用 myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS) 中，你给它一个不是真实投影的投影。声明投影和转换投影之间有很大的区别。当你使用投影时，让函数从文件中读取投影，然后总是使用适当的函数(spTransform, projectRaster 或 gdalwarp)进行转换。 - Bastien