在 R 中使用多边形裁剪栅格:错误范围不重叠
Crop raster with polygon in R: Error extent does not overlap
我想使用在 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
该消息非常不言自明。你的范围不重叠......这里如何检查它:
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 的方式。而不是
myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)
使用:
library(rgdal)
myExtent <- readOGR("C:/mydir","loc1.shp")
myExtent <- spTRansform(myExtent, CRS(proj4string(test1)))
我想使用在 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
该消息非常不言自明。你的范围不重叠......这里如何检查它:
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 的方式。而不是
myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)
使用:
library(rgdal)
myExtent <- readOGR("C:/mydir","loc1.shp")
myExtent <- spTRansform(myExtent, CRS(proj4string(test1)))