Bug très probable dans la fonction d'intersection du paquet R raster en croisant 1 polygone avec 1 point

Question

En essayant d'obtenir un résultat d'intersection à partir d'un seul point et d'un seul polygone, j'ai trouvé ce que je crois ne peut être qu'un bogue dans le raster fonction d'intersection de paquet.

J'ai 1 polygone et 1 point, et utiliser Intersection comme suit:

intersect(a_point, a_polygon) 

Où a_point contient un attribut id. Cela échoue avec l'erreur:

Error in j[, 2] : incorrect number of dimensions

Cependant, si je renverse les arguments et faire:

intersect(a_polygon, a_point) 

Il fonctionne très bien, mais ne renvoie pas l'identifiant du fichier de forme de points dans le cadre du résultat dont j'ai besoin. C'est un comportement attendu, tellement bien mais j'ai besoin de travailler dans l'autre sens. Pour exclure certaines particularités de mon polygone ou de mes données ponctuelles, j'ai créé un seul polygone et un seul objet spatial ponctuel et testé la même hypothèse, et le même résultat s'est produit comme ci-dessus avec ces objets «bruts».

Ce qui suit est le code pour générer ces deux objets « faux » pour l'exhaustivité et pour qu'il puisse être reproduit:

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) 

Maintenant, voici la chose, si je crée 2 points à côté de l'autre (si ils sont tous les deux dans le polygone) au lieu d'un seul, ça marche bien, pas d'erreur. Suggérant qu'il y a un bug associé à l'intersection entre 1 point et 1 polygone, QUAND le point porte un attribut à retourner dans le processus d'intersection. Vous pourriez vous demander pourquoi vous avez réellement besoin d'avoir l'attribut retourné s'il y a juste un point, c'est parce que c'est un processus itératif dans lequel il peut ne pas être juste un point, il pourrait être aucun ou plusieurs.

J'apprécierais que quelqu'un explique cette erreur ou confirme mes conclusions.

Answer 1

Désolé je ne peux pas répondre à votre question bug Intersection, mais il pourrait être plus simple pour l'instant d'utiliser sp :: plus revenir polygone attributs aux points

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

Answer 2

Voici vos données d'exemple dans un plus concis façon.

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

Maintenant, illustrez le problème.

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 

Ce fut un autre bug drop=TRUE causé par le défaut R de matrices « dropping » à des vecteurs lorsqu'une seule ligne est sélectionnée. Cela a été corrigé dans la version raster 2.6-11 (pas encore sur CRAN).