Intersect points and polygons

Question

I have a dataset where I'm trying to look determine determine which polygons my sampling locations occur in. I have been following the example here:

https://gis.stackexchange.com/questions/282750/identify-polygon-containing-point-with-r-sf-package

But I'm running into error messages I don't understand, nor do I follow all the code.

#polygon shapefile information
bcr
Simple feature collection with 378 features and 1 field
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -179.1413 ymin: 14.53287 xmax: -52.63629 ymax: 83.11063
geographic CRS: NAD83
First 10 features:

                       BCRNAME                       geometry
1  ALEUTIAN/BERING SEA ISLANDS MULTIPOLYGON (((-179.1036 5...
2               WESTERN ALASKA MULTIPOLYGON (((-162.4558 5...
3  ARCTIC PLAINS AND MOUNTAINS MULTIPOLYGON (((-162.6464 6...
4  ARCTIC PLAINS AND MOUNTAINS MULTIPOLYGON (((-94.83377 5...
5  ARCTIC PLAINS AND MOUNTAINS MULTIPOLYGON (((-61.63853 5...
6  ARCTIC PLAINS AND MOUNTAINS MULTIPOLYGON (((-101.9977 6...
7  ARCTIC PLAINS AND MOUNTAINS MULTIPOLYGON (((-102.0227 6...
8  ARCTIC PLAINS AND MOUNTAINS MULTIPOLYGON (((-101.9978 6...
9  ARCTIC PLAINS AND MOUNTAINS MULTIPOLYGON (((-79.70593 5...
10 ARCTIC PLAINS AND MOUNTAINS MULTIPOLYGON (((-69.77826 5...

#select distinct set of points
pnts <- df %>% select(y = Latitude, x = Long) %>% distinct()

#create a points collection 
pnts_sf <- do.call("st_sfc", 
                   c(lapply(1:nrow(pnts), 
                            function(i) {st_point(as.numeric(pnts[i, ]))}), list("crs" = 4326))) 
#My question: what is this doing?? what does 4326 come from?

pnts_trans <- st_transform(pnts_sf, 2136) # apply transformation to pnts sf
bcr_trans <- st_transform(bcr, 2163) #apply transformation to polygon sf

#intersect and extract BCR name
pnts$bcrName <- apply(st_intersects(bcr_trans, pnts_trans, sparse = FALSE), 2, 
                     function(col) { 
                       bcr_trans[which(col), ]$BCRNAME
                     })
#My question: So this is where we're looking at the interaction between the points and polygons, but what is 2 for and what is the function doing?

When I try to run everything, I get the following error message:

Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared,  : 
  st_crs(x) == st_crs(y) is not TRUE

Unclear if it's coming from a variable name I didn't change or something else entirely.

Answer 1

I'm afraid you're making this more complicated for yourself. I have created reproducible code that illustrates what you should do step by step. See below:

EDIT: After seeing what you were really looking for, which is essentially the answer noted in the other stack exchange post you linked here, I have added a reproducible example that should be able to guide you. See below. Note, I am leaving my previous answer below this one just in case it is helpful.

library(sf)
library(tmap)
library(dplyr)

data(World) # free dataset that will load
countries = World %>%
  filter(name == "Afghanistan" | 
           name == "Turkmenistan" |
           name == "Iran") #Picking three countries in the world

pnts = data.frame(name=c("Point1","Point2","Point3","Point4"),latitude = c(34,35,36,37),longitude=c(63,65,62.5,66.5)) #Picking random points

sp_points = st_as_sf(pnts,coords = c('longitude',"latitude"))#make points spatial
st_crs(sp_points)= 4326 # Give the points a coordinate reference system (CRS)
sp_points=st_transform(sp_points,crs = st_crs(countries)) # Match the point and polygon CRS

tm_shape(countries)+
  tm_borders() +
  tm_text('name')+
  tm_shape(sp_points)+
  tm_dots(col='red',size=1.2)+
  tm_text('name',ymod = 1) #Creates the map below

sp_points$country_name <- apply(st_intersects(countries, sp_points, sparse = FALSE), 2, 
                     function(col) {countries[which(col), ]$name}) 
# Just make sure you're correctly placing your polygon objects and point objects. 
# In this case, my point object is sp_points, and my polygon object is countries

sp_points
#> Simple feature collection with 4 features and 2 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 5413377 ymin: 4321480 xmax: 5732135 ymax: 4670058
#> CRS:            +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#>     name                geometry country_name
#> 1 Point1 POINT (5507002 4321480)  Afghanistan
#> 2 Point2 POINT (5656241 4438573)  Afghanistan
#> 3 Point3 POINT (5413377 4554775) Turkmenistan
#> 4 Point4 POINT (5732135 4670058)  Afghanistan

Old answer below, which is a way to visualize the same operation, with a completely reproducible example.

library(sf)
library(tmap)
library(dplyr)
data(World)
countries = World %>%
    filter(name == "Afghanistan" | 
             name == "Turkmenistan" |
                name == "Iran") #Picking three countries in the world

pnts = data.frame(latitude = c(34,35,36,37),longitude=c(63,65,62.5,66.5)) #Picking random points

sp_points = st_as_sf(pnts,coords = c('longitude',"latitude"))#make points spatial
st_crs(sp_points)= 4326 # Give the points a coordinate reference system (CRS)
sp_points=st_transform(sp_points,crs = st_crs(countries)) # Match the point and polygon CRS

tm_shape(countries)+
  tm_borders() +
tm_shape(sp_points)+
  tm_dots(col='red',size=1.2)

#Notice that Iran does not contain a dot

newcountries=st_intersects(countries,sp_points) #preforms the intersection
newc_logical=lengths(newcountries)>0 #Finds only countries where intersection is TRUE
newcountries2=countries[newc_logical,] #Applies it to the original countries
#The three above lines were found from: 
#https://cran.r-project.org/web/packages/sf/vignettes/sf3.html


tm_shape(newcountries2)+
  tm_borders() #Leaves you with the df newcountries2 without the country containing a dot

To answer your question about 4326, it is the coordinate reference system. This coordinate reference system is probably the most standard one for mapping, most GIS generally revert to it. The number is specifically called an EPSG code.

What you should focus on is getting both the polygon and points you are working with on the same coordinate reference system, and make sure that they are both spatial features (hence sf). The code I have showed is the most simple version of an intersection. Hopefully your problem will be simple enough to be handled this way.