在 R 中使用多边形裁剪栅格:错误范围不重叠

dta*_*non 3 r crop raster image-processing

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

首先我创建光栅堆栈:

test1 < stack("C:/mydir/test1.tif")
Run Code Online (Sandbox Code Playgroud)

定义投影

myCRS <- test1@crs
Run Code Online (Sandbox Code Playgroud)

然后读取shapefile

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

庄稼

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

我已经寻找了一个解决方案,但我只发现投影可能是问题所在,但是它们肯定都在同一个 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 
Run Code Online (Sandbox Code Playgroud)

Bas*_*ien 5

该消息是不言自明的。你的范围不重叠......这里如何检查它:

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")
Run Code Online (Sandbox Code Playgroud)

我已经根据您问题中的数据创建了范围对象。您会在左上角看到一个范围,而在右下角看到另一个范围。你有没有试过在 QGIS 中阅读这两个文件,你可能很容易看到它。

在此处输入图片说明

如果它们真的应该重叠,那么我会怀疑您阅读 shapefile 的方式。代替

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

用 :

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

  • 这可能是因为,在您调用 `myExtent &lt;- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)` 时,您给了它一个不是真实的投影。陈述投影和变换投影之间有很大的区别。当你玩投影时,让函数从文件中读取投影,而不是总是使用适当的函数(`spTransform`、`projectRaster` 或 `gdalwarp`)进行转换 (2认同)