CODEWITHSUNDEEP
r
user33125
user33125

Reputation: 197

Intersect two shapefiles with R

I think I want to do something very basic, but I seem to not find posts on how to do it. In case I am wrong and there has already been a post on this topic, I apologize!

I have two shapefiles which you can find through this link: 

setwd("~/where you saved your data")
nuts <- readOGR(".", layer = "NUTS_RG_60M_2010")
aqueduct <- readOGR(".", layer = "aqueduct_global_dl_20150409")

And to be honest, now I am already stuck. I would like to add the values of all the variables of the aqueduct shapefile, to the intersecting nuts3 regions of the nuts file.

I tried gIntersection, intersect, extract... but without success. Could somebody please help me getting the intersection correct? The end results is a nuts-shapefile with all the variables of the aqueduct-shapefile.

Thank you very much!

Upvotes: 2

Views: 4761

Answers (1)

Michael Dorman
Michael Dorman

Reputation: 1020

This is a spatial join of two polygonal layers. Unless each feature in nuts intersects exactly with one feature of aqueduct, there is no straightforward/single way to do the spatial join.

Instead, you can either obtain a list of rows from the attribute table of aqueduct corresponding to each feature of nuts - 

nuts_over1 = over(nuts, aqueduct, returnList = TRUE)

Or summarize specific attributes with a function, in which case the result can be joined back to the attribute table of nuts. For example, the countries of aqueduct intersecting the six first features are as follows - 

nuts_over2 = over(
  nuts, 
  aqueduct[, "COUNTRY"], 
  fn = function(x) paste(x, collapse = ", ")
)
head(nuts_over2)
                                                          COUNTRY
1                              Austria, Hungary, Austria, Hungary
2 Austria, Hungary, Austria, Hungary, Slovakia, Austria, Slovakia
3                     Austria, Austria, Hungary, Austria, Hungary
4                                                         Austria
5                     Austria, Austria, Austria, Austria, Austria
6                                                Austria, Austria

This information can be joined back to the attribute table of nuts as follows - 

nuts@data = cbind(nuts@data, nuts_over2)
head(nuts@data)
  NUTS_ID STAT_LEVL_ SHAPE_Leng SHAPE_Area
0   AT111          3   1.089017 0.08091455
1   AT112          3   2.257319 0.20926007
2   AT113          3   2.002492 0.17728455
3   AT121          3   3.158370 0.40147321
4   AT122          3   2.956927 0.42675504
5   AT123          3   2.010415 0.14145865
                                                          COUNTRY
0                              Austria, Hungary, Austria, Hungary
1 Austria, Hungary, Austria, Hungary, Slovakia, Austria, Slovakia
2                     Austria, Austria, Hungary, Austria, Hungary
3                                                         Austria
4                     Austria, Austria, Austria, Austria, Austria
5                                                Austria, Austria

Upvotes: 2

Related Questions