title: A function for determining whether a series of points and a series polygons intersect
body: <p>I have a list of points and a list of polygons, and I wish to find the polygon each point locates in. Though it is straightforward to use the command <code>st_intersects</code> in <code>sf</code> package, I find that is not good here since both the points and polygons are so many that the resulting matrix is too big for my computer to work.</p>
  
  <p>I try to get around this by making a function, which can look through all the polygons to find the right one for each point. I made one but it did not work well and I have had no idea what is the problem with it. So I need some help! Thanks!</p>
  
  <p>The following is how I create the list of points</p>
  
  <pre><code>pts = st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5, 
  2.5))) %&gt;% st_sf() %&gt;% mutate(namepoint=c(1,2,3))
  
  Simple feature collection with 3 features and 1 field
  geometry type:  POINT
  dimension:      XY
  bbox:           xmin: 0.5 ymin: 0.5 xmax: 2.5 ymax: 2.5
  epsg (SRID):    NA
  proj4string:    NA
     namepoint        geometry
  1         1 POINT (0.5 0.5)
  2         2 POINT (1.5 1.5)
  3         3 POINT (2.5 2.5)
  </code></pre>
  
  <p>The following is how I create the list of polygons</p>
  
  <pre><code>pol_down = st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0))))
  pol_up=pol_down+c(2,1) 
  pol_add &lt;- list(pol_down,pol_up) %&gt;% st_sfc() %&gt;% st_sf %&gt;% 
  mutate(namepgn=c('A','B'))
  
  Simple feature collection with 2 features and 1 field
  geometry type:  POLYGON
  dimension:      XY
  bbox:           xmin: 0 ymin: 0 xmax: 4 ymax: 3
  epsg (SRID):    NA
  proj4string:    NA
     namepgn                       geometry
  1       A POLYGON ((0 0, 2 0, 2 2, 0 ...
  2       B POLYGON ((2 1, 4 1, 4 3, 2 ...
  </code></pre>
  
  <p>The function I made is </p>
  
  <pre><code> fun_pgn &lt;- function(x,y){
                  location_m&lt;- st_intersects(x,
                                             y,
                                             sparse = F)
                  name &lt;- pull(y,1) %&gt;% .[location_m]
                  return(name)
               }
  </code></pre>
  
  <p>Use the command <code>st_intersects</code> to get the their relations</p>
  
  <pre><code>st_intersects(pts,pol,sparse = F)
  
        [,1]  [,2]
  [1,]  TRUE FALSE
  [2,]  TRUE FALSE
  [3,] FALSE  TRUE
  </code></pre>
  
  <p>Then what I want is to add a variable,like <code>pgnname</code> in the data <code>pts</code> to record the polygon for each point,</p>
  
  <pre><code>        namepoint    pgnname    geometry
  1               1       "A"     POINT (0.5 0.5)
  2               2       "A"     POINT (1.5 1.5)
  3               3       "B"     POINT (2.5 2.5)
  </code></pre>
  
  <p>However my function does not work. After running the following codes</p>
  
  <pre><code>pts %&gt;% mutate(pgnname=fun_pgn(geometry,pol))
  </code></pre>
  
  <p>it follows </p>
  
  <pre><code>Simple feature collection with 3 features and 2 fields
  geometry type:  POINT
  dimension:      XY
  bbox:           xmin: 0.5 ymin: 0.5 xmax: 2.5 ymax: 2.5
  epsg (SRID):    NA
  proj4string:    NA
    namepoint pgnname        geometry
  1         1       A POINT (0.5 0.5)
  2         2       B POINT (1.5 1.5)
  3         3    &lt;NA&gt; POINT (2.5 2.5)   
  </code></pre>
  
  <p>The names of polygons are incorrect! But the key part of my function works well. For example  </p>
  
  <pre><code>&gt; pull(pol,1) %&gt;% .[st_intersects(pts[3,],pol,sparse = F)]
  [1] "B"
  </code></pre>
  
  <p>So I'm confused. What happens when I run my functions?</p>