当1个多边形与1个点相交时,R栅格包交叉函数中很可能出现错误

问题描述 投票:1回答:2

通过尝试从单个点和单个多边形获得交叉结果,我发现了我认为只能是R栅格包交叉函数中的错误。

我有1个多边形和1个点,并使用交叉如下:

intersect(a_point, a_polygon)

其中a_point包含id属性。这失败并出现错误:

j [,2]出错:尺寸数不正确

但是,如果我反转参数并执行:

intersect(a_polygon, a_point)

它工作正常,但不会返回点形状文件中的id作为我需要的结果的一部分。这是预期的行为,很好,但我需要它以相反的方式工作。

为了排除我的多边形或点数据的某些特性,我创建了一个多边形和单点空间对象并测试了相同的假设,并且使用这些“原始”对象产生了与上述相同的结果。

以下是生成这两个“假”对象的完整性的代码,以便可以重现:

test_list_x = list(530124, 530125)  #For when I use 2 points 
test_list_y =  list(176949, 176950) #For when I use 2 points 

data_frame_object = data.frame(530124, 176950)
names(data_frame_object) = c("Longitude", "Latitude")
coordinates(data_frame_object)=~Longitude+Latitude
proj4string(data_frame_object)=CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894") 
fake_point_shape_object=SpatialPointsDataFrame(data_frame_object, data.frame(id=1:length(data_frame_object)))


coords = matrix(  nrow=5,  ncol=2) 
coords[1,1] =  530106.8
coords[1,2] =  176953.3
coords[2,1] = 530127.5
coords[2,2] =  176953.3
coords[3,1] = 530127.5
coords[3,2] =  176933.3
coords[4,1] = 530106.8
coords[4,2] =  176933.3
coords[5,1] = 530106.8
coords[5,2] = 176953.3
my_fake_polygon = Polygon(coords)

polygon_list = list(my_fake_polygon)

polygon_set <- lapply(seq_along(polygon_list), function(i) Polygons(list(polygon_list[[i]]), i  ))

new_polygons <- SpatialPolygons(polygon_set)
new_polygons@proj4string = CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894") 

df <- data.frame("1")
names(df) = "id"

my_fake_polygon <- SpatialPolygonsDataFrame(new_polygons,df)

现在就是这样,如果我创建彼此相邻的2个点(因此它们都在多边形内)而不是仅仅一个,它工作正常,没有错误。建议在1点和1个多边形之间存在与交叉点相关的错误,此时该点携带要在交叉点过程中返回的属性。

您可能会问为什么实际上如果只有一个点,您确实需要返回属性,这是因为它是一个迭代过程,它可能不是一个点,它可能是无或多个。

我很感激有人解释这个错误或确认我的发现。

r raster r-raster
2个回答
1
投票

以下是更简洁的示例数据。

library(raster)

pnt <- SpatialPoints(cbind(530124, 176950))
pol <- spPolygons(matrix(c(530106.8, 530127.5, 530127.5, 530106.8, 530106.8, 176953.3, 176953.3, 176933.3, 176933.3, 176953.3), ncol=2))

现在说明问题。

intersect(pol, pnt)
#class       : SpatialPolygons 
#features    : 1 
#extent      : 530106.8, 530127.5, 176933.3, 176953.3  (xmin, xmax, ymin, ymax)
#coord. ref. : NA 

# this fails
intersect(pnt, pol)
#Loading required namespace: rgeos
#Error in j[, 2] : incorrect number of dimensions

# but it works with two points!
intersect(bind(pnt, pnt), pol)
#class       : SpatialPoints 
#features    : 2 
#extent      : 530124, 530124, 176950, 176950  (xmin, xmax, ymin, ymax)
#coord. ref. : NA 

这是另一个drop=TRUE错误,当选择单行时,由于将矩阵“丢弃”到向量的R默认值而导致。这在光栅版本2.6-11中修复(尚未在CRAN上)。


0
投票

抱歉,我无法回答您的交叉错误问题,但现在使用sp :: over将多边形属性返回到点可能更简单

# dummy polygon
xym <- as.matrix(data.frame(x=c(16.48438,17.49512,24.74609,22.59277,16.48438),
y=c(59.73633,55.12207,55.03418,61.14258,59.73633)))

# make into SpatialPolygon
p = Polygon(xym)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))

# Promote to SPDF and give an attribute
SPDF = SpatialPolygonsDataFrame(sps, data.frame(N = "hello", row.names = 1))

# make 2 points, one inside the polygon and one outside
p <- data.frame(x=c(16,18),y=c(58,58))
coordinates(p) <- ~x + y

# plot to check
plot(sps)
plot(p,add=T)

# perform the over, returns a named vector for every point in the SpatialPoints
res <- unname(over(p,SPDF))

# promote points to SpatialPointsDataFrame and put in new polygon attribute
data <- data.frame(ID=row.names(p),pol=res)
sp <- SpatialPointsDataFrame(p, data)
© www.soinside.com 2019 - 2024. All rights reserved.