When joining attributes by location with the over function in R from the sp package, only the values for one of the polygons that intersects is retained. Is there a way in R that allows to take attributes of more than located feature? That is, to get ALL qualitative values the polygons that intersect/overlap. With an attribute column for the value of the first located feature, another for the second located feature, etc.

I have also asked in a separate question about Joining attributes by location for more than one located feature using QGIS

Sample code to perform simple spatial joins

#set up sample data
    install.packages('tmap')#use this library to upload sample spatial polygons
    data(Europe) # countries
    data(metro) # cities
    #transform sample data of points into polygons with same crs/projection for this example (trivial example because point.in.poly suited to this problem but just for the sake of an example
    metro.pts=spTransform(metro, crs(Europe))
    metro.poly=gBuffer(metro.pts, width=50000, byid=TRUE)
#load libraries for spatial join
 #spatial join of the two polygon shapefiles
        #"find the cities in each country"
        joined_one=cbind(over(Europe, metro.poly), as.data.frame(Europe))
    #but this only joins one of the cities in each country even when more than one city is present
    #if we check for spain, we only have Barcelona and not Madrid even though it is part of the city file

    #using the `fn` option (function) we can calculate arithmetic for quantitative variables as to get values for more than one polygon that intersects
    joined_qn=cbind(over(Europe, metro.poly[,5:8], fn=mean), as.data.frame(Europe))

2 Answers 2


You can use raster::union

Example data:

# simplifying a bit
Eu <- Europe[Europe$iso_a3 %in% c('ESP', 'PRT', 'FRA'),1:2]
metro.pts <- spTransform(metro[, 'name'], crs(Europe))
metro.poly <- gBuffer(metro.pts, width=50000, byid=TRUE)

Original solution

u <- union(Eu, metro.poly)


i <- intersect(Eu, metro.poly) 

**Update, now I understand the question better **

You can returnList = TRUE in over:

v <- over(Eu, metro.poly, returnList=TRUE)

You could now do:

a <- lapply(1:length(v), function(i) cbind(iso_a3=Eu$iso_a3[i], v[[i]]))
b <- do.call(rbind, a)
m <- merge(Eu, b, by="iso_a3", all.x=TRUE, duplicateGeoms=TRUE)

In merge, the argument duplicateGeoms=TRUE is necessary because otherwise merge cannot handle having more than one (i.e. duplicate) value to assign to each feature.
To get the following:

#  iso_a3   name.x    name.y
#1    ESP    Spain Barcelona
#2    ESP    Spain    Madrid
#3    FRA   France     Lille
#4    FRA   France      Lyon
#5    FRA   France Marseille
#6    FRA   France     Paris
#7    PRT Portugal    Lisbon
#8    PRT Portugal     Porto

I think that this is the same as what Tiernan showed in perhaps more elegant sf code, but I do not think that this is what you asked for.

For what you want, we need to create a data.frame from a list with varying number of elements. There might be a better way, but here is one:

m <- max(sapply(v, nrow))
iso <- as.character(Eu$iso_a3)

x <- sapply(1:length(v), 
       function(i) c(iso[i], c(v[[i]][,1], rep(NA, m))[1:m])

x <- data.frame(t(x), stringsAsFactors=FALSE)
colnames(x) <- c("iso_a3", paste0("metro_", 1:m))

Merge with the original polygons:

r <- merge(Eu, x, by="iso_a3", all.x=TRUE)

#  iso_a3     name   metro_1 metro_2   metro_3 metro_4
#1    ESP    Spain Barcelona  Madrid      <NA>    <NA>
#2    FRA   France     Lille    Lyon Marseille   Paris
#3    PRT Portugal    Lisbon   Porto      <NA>    <NA>
  • I have updated my answer. I believe that the answer by Tiernan, that you marked as correct does not seem to do what you ask for, as it creates new records for each country/metro combination? Also, I assume that with "shapefile" you mean polygons Feb 8, 2018 at 23:26
  • You're new solution does indeed fulfill my request perfectly I did ask for columns and not additional rows, but the last piece of converting to multiple rows to additional columns (missing in @Tiernan 's answer) was something I knew how to do. I realize now that I need to learn how to manipulate lists and sapply. I tend to get overwhelmed with a list output but your answer illustrates how to manipulate them. Feb 9, 2018 at 14:46

I recommend you consider transitioning from sp to the sf package (this handy guide will help you along the way).

Here's an approach that uses sf::st_join() to achieve the spatial join you're looking for:



europe_sf <- st_as_sf(Europe)

metro_pts_sf <- metro %>%
  st_as_sf() %>%

metro_poly_sf <- st_buffer(metro_pts_sf, dist = 5e4)

joined_sf <- st_join(europe_sf, metro_poly_sf)

Let's check the two Spanish cities you mention in the question:

joined_sf %>%
  filter(name.x %in% "Spain") %>%
  select(matches("name"), everything())

## Simple feature collection with 2 features and 28 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 2766833 ymin: 1584304 xmax: 3833320 ymax: 2463100
## epsg (SRID):    NA
## proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
##   name.x    name.y name_long iso_a3.x sovereignt continent            part
## 1  Spain Barcelona Barcelona      ESP      Spain    Europe Southern Europe
## 2  Spain    Madrid    Madrid      ESP      Spain    Europe Southern Europe
##   EU_Schengen   area  pop_est pop_est_dens gdp_md_est gdp_cap_est
## 1 EU Schengen 498800 40525002     81.24499    1403000     34620.6
## 2 EU Schengen 498800 40525002     81.24499    1403000     34620.6
##                      economy           income_grp life_exp well_being
## 1 2. Developed region: nonG7 1. High income: OECD     81.4   6.188263
## 2 2. Developed region: nonG7 1. High income: OECD     81.4   6.188263
##        HPI iso_a3.y pop1950 pop1960 pop1970 pop1980 pop1990 pop2000
## 1 44.06279      ESP 1809390 2467926 3482047 3836761 4100808 4355425
## 2 44.06279      ESP 1699752 2391543 3520861 4253190 4413870 5014411
##   pop2010 pop2020 pop2030                       geometry
## 1 4933548 5478238 5685152 MULTIPOLYGON (((3830463 187...
## 2 5787392 6476334 6707421 MULTIPOLYGON (((3830463 187...
  • 1
    It is fair to recommend using both packages however, sf is not a replacement for sp just yet. Especially, if you consider critical dependencies such as spdep. Whereas, it is quite easy to coerce back-and-forth there is an associated overhead in doing so. Feb 7, 2018 at 23:38
  • 1
    Not sure what your point is... is it poor form to suggest sf when a user poses a question using sp? This seems like a good opportunity to nudge a user toward the more modern implementation of spatial vector data (i.e., simple features); one which is being rapidly implemented by a huge number of other supporting packages, albeit not spdep. ¯_(ツ)_/¯
    – Tiernan
    Feb 7, 2018 at 23:50
  • @Jeffrey Evans I find it aggravating to get used to one set of functions (in this case sp) only to need to switch to another set of functions with their own language and syntax (in this case sf but it also applies to my frustration trying to use ggplots2). BUT, if the functionality I need is only available in one language... then I will make the effort to use another set. If spatialEco allowed me to extract more than one species name for each forest polygon then I would just keep to sp... hint hint Feb 8, 2018 at 16:01
  • I feel that the tidyverse path is the way forward but, recoding entire packages to leverage new object classes is not trivial. In the spatial packages there are actually very few that have moved to an sf class and the ones that have are doing it through object coercion as(sf "spatial"), which is how I am currently dealing with it so sf objects do not fail. I am looking at this function today and will address the multi-polygon problem but, keep in mind that it is a many to one match. In this case I would likely output the polygon IDS as columns and not individual attributes. Feb 8, 2018 at 16:18
  • 2
    @user3386170, hint taken, I just rewrote the spatialEco::point.in.poly function so that it accounts for multiple polygon intersections. In these cases the point data is recycled so, more points are returned than go in to allow for merging of attributes. These duplicated points represent areas where more than one polygon intersects a given point. So, if three polygons intersect a given point then that point will occur in the returned data three times, with each associated polygon attribute. I should have the new package version up on CRAN in a few days. Feb 8, 2018 at 18:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.