1

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.

1
  • Did the answer below correctly answer your question? Feb 8 at 18:22
2

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.

5
  • Thanks @coconn41. That didn't quite answer my question as I need to create a new variable the the name of the (in your example) country for each point. But, I can definitely see how this approach would work for other analyses.
    – tnt
    Feb 11 at 17:29
  • @tnt I'm not sure you're being clear enough in your question. What exactly are you trying to do? Your question says that you're "...trying to look determine determine which polygons my sampling locations occur in.". Isn't that what I showed you? Or are you trying to take your points dataframe and create a variable within that dataframe that tells you which polygon the point intersects with? Feb 11 at 19:17
  • I'm trying to create a variable in my dataframe with the name of the polygons the points are located within.
    – tnt
    Feb 11 at 20:10
  • @tnt Please see my edited answer, and let me know if you still can't get yours to work. Feb 11 at 20:30
  • Thanks! This worked perfectly.
    – tnt
    Feb 12 at 22:46

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.