2

I'm working on converting an old script from using sp() to sf() and some of the computations have to be rewritten, including some operations which (used to) use cbind with over() like so:

cbind(df1data, over(df1 df2))

I've tried a number of combinations of st_intersects, st_join, and merge_cols, and can't seem to get the columns added on from the second sf().

Here's a reproducible example:

# read in relevant polygons, Scottish Index of Multiple deprivation
if (file.exists("data/sc_dz_11.shp") == FALSE) {
download.file("http://simd.scot/2016/data/simd2016_withgeog.zip", 
              destfile = "data/simd2016_withgeog.zip")
unzip("data/simd2016_withgeog.zip", exdir = "data", junkpaths = TRUE)
}
simd_shapes <- st_read("data/sc_dz_11.shp")
simd_indicators <- read.csv("./data/simd2016_withinds.csv")
simd_indicators_min <- simd_indicators[c(1,6:17)]
# Original data is in wgs, so need to convert from wgs to bng to work with analysis across data sets
simd_wgs <- merge(x=simd_shapes, y=simd_indicators_min, by.x = "DataZone", by.y = "Data_Zone")
simd <- simd_wgs %>% st_transform(27700)

c1 <- c("site1","site2","site3")
c2 <- c("357109","350477.6", "257916")
c3 <- c("734754","614926.268", "662462")
df <- data.frame(c1, c2, c3, c3)

df_joined <- cbind(df@data, over(df, simd)) # of course this is the part where I'm struggling
df_joined <- bind_cols(df,st_intersects(df,simd)) # tried this too and it also didn't work

Very grateful if anyone can shine some light on what I'm getting wrong here!

3
  • On this context, sf can work as a regular data.frame. Have you tried cbind(simd, df )/ bind_cols(simd, df )? If both datasets have the same number of rows it should work
    – dieghernan
    Mar 4 at 19:42
  • That's the trouble, simd has far more rows and needs to be filtered by intersection. Sorry, just realised I forgot to add crs to my example df - will fill that in. Mar 4 at 19:44
  • Maybe this may help you (st_intersection) if both simd and df are sf objects (not the case on your example, df is a regular data frame, isn’t it?) stackoverflow.com/questions/71344590/…
    – dieghernan
    Mar 4 at 19:49

1 Answer 1

2

Both sp and sf objects require an indication of which columns are X and Y, in this case that part was not done in the code, once you have a sf object, joins are natural. I did some polishing in the code, to make it reproducible:

library(sf) 
library(dplyr)

dir.create("data") # so all dirs in IF statement work
if (file.exists("data/sc_dz_11.shp") == FALSE) {
  download.file("https://simd.scot/data/simd2016_withgeog.zip",   # URL had changed
                destfile = "data/simd2016_withgeog.zip")
  unzip("data/simd2016_withgeog.zip", exdir = "data", junkpaths = TRUE) 
}


simd_shapes <- st_read("data/sc_dz_11.shp")
simd_indicators <- read.csv("data/simd2016_withinds.csv")
simd_indicators_min <- simd_indicators[c(1,6:17)]

# let's also make this tidyverse 
simd_wgs = left_join(simd_shapes, simd_indicators_min, by = c("DataZone" = "Data_Zone"))

simd <- simd_wgs %>% st_transform(27700)

c1 <- c("site1","site2","site3")
c2 <- c("357109","350477.6", "257916")
c3 <- c("734754","614926.268", "662462")
df <- data.frame(c1, c2, c3, c3)

# convert to sf object
# |> is the R's pipe, for versions > 4
df = df |> st_as_sf(coords = c("c2", "c3")) |> st_set_crs(27700)

# now we join: 
df_joined = st_join(df, simd )

#Simple feature collection with 3 features and 23 fields
#Geometry type: POINT
#Dimension:     XY
#Bounding box:  xmin: 257916 ymin: 614926.3 xmax: 357109 ymax: 734754
#Projected CRS: OSGB 1936 / British National Grid
#     c1       c3.1  DataZone                         Name TotPop2011 #ResPop2011 HHCnt2011 StdAreaHa StdAreaKm2 Shape_Leng Shape_Area
#1 site1     734754 S01007166         Carnoustie East - 04        709        #707       328  20.25832   0.202583   2620.603   202583.2
#2 site2 614926.268 S01012366 Hawick Central - Town Centre        864        #862       513  16.88552   0.168856   2706.169   168855.2
#3 site3     662462 S01009909             Strathbungo - 05        849        #821       471  21.71528   0.217151   4074.232   217152.8


1
  • Yep - that's exactly the ticket. Thanks so much! Also appreciate the tidyverse conversions in there. Mar 5 at 7:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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