首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >在R中裁剪带有面的栅格:错误范围不重叠

在R中裁剪带有面的栅格:错误范围不重叠
EN

Stack Overflow用户
提问于 2017-12-19 18:37:55
回答 1查看 7.3K关注 0票数 3

我想使用我在ArcGIS中创建的多边形形状文件来裁剪光栅堆栈,但是我得到了范围不重叠的错误。

首先,我创建栅格堆栈:

代码语言:javascript
运行
复制
test1 < stack("C:/mydir/test1.tif")

定义投影

代码语言:javascript
运行
复制
myCRS <- test1@crs

然后读取shapefile

代码语言:javascript
运行
复制
myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)

裁剪

代码语言:javascript
运行
复制
myCrop <- crop(test1, myExtent)
Error in .local(x, y, ...) : extents do not overlap

我一直在寻找解决方案,但我只发现投影可能是问题所在,然而它们肯定都在同一个CRS中:

代码语言:javascript
运行
复制
> 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 
EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2017-12-19 20:31:49

这个信息是不言而喻的。您的范围不重叠...检查方法如下:

代码语言:javascript
运行
复制
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中读取这两个文件,您可能很容易就能看到它们。

如果它们真的是重叠的,那么我会怀疑你读取shapefile的方式。而不是

代码语言:javascript
运行
复制
myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)

使用:

代码语言:javascript
运行
复制
library(rgdal)
myExtent <- readOGR("C:/mydir","loc1.shp")
myExtent <- spTRansform(myExtent, CRS(proj4string(test1)))
票数 5
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/47885065

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档