R中具有多边形的裁剪栅格：错误范围不重叠 [英] Crop raster with polygon in R: Error extent does not overlap

问题描述

我想使用我在ArcGIS中制作的多边形shapefile裁剪栅格堆栈，但是遇到范围不重叠的错误。

I want to crop a raster stack using a polygon shapefile i made in ArcGIS, however I get error that extent does not overlap.

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

First I create the raster stack:

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

定义投影

myCRS <- test1@crs

然后读取shapefile

then read 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中完全相同：

I have searched for a solution, but i only find that projection can be the problem, however they are definetly both in the same 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 

推荐答案

该消息很容易解释。您的范围不重叠...在这里如何检查它：

The message is pretty self explanatory. Your extent do not overlap... here how to check it:

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

我已经根据问题中的数据创建了范围对象。您在左上角看到一个范围，在右下角看到另一个。您尝试读取QGIS中的两个文件，您可能会很容易看到它。

I've created extent object from data in your question. You see one extent in the top left corner and the other in the bottom right. Have you tried reading both files in QGIS, you could probably easily see it.

如果它们确实应该重叠，那么我会怀疑您读取shapefile的方式。而不是

If they really are suppose to overlap, than I would suspect the way you read your shapefile. Instead of

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

使用：

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

这篇关于R中具有多边形的裁剪栅格：错误范围不重叠的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！