# Most efficient way to find points in a polygon (polygon always rectangular)

This question has been asked a lot by a lot of different people I am sure. I did some searching and found some blog posts and potential ways forward, but before I dove off into the deep end I thought I would come here to get the community consensus.

• Problem Statement

I have a series of points (approximately 1.5 million) which I need to test for inclusion in rectangular polygon. When I profile my code, this is the longest pole in the tent when it comes to execution time. I would like to see if there is something that I can do which would reduce the execution time.

• Minimal code example and Offending Code Line
``````library(sf)
points_to_include = unlist(sf::st_contains(profile_segments\$buffer[[1]], bathy_data))
``````

(The above code and the input.RData file can be found at this GitHub Gist)

This code is currently taking ~ 6 seconds to run based on `microbenchmark`. While that does not seem like a long time, doing this over and over really begins to add up.

``````Unit: seconds
expr      min       lq     mean   median       uq      max neval
points_to_include 5.494459 5.532229 5.985662 5.591614 6.060509 7.249499     5
``````
• Potential Solutions

I looked at the blog post https://www.r-bloggers.com/speeding-up-spatial-analyses-by-integrating-sf-and-data-table-a-test-case-2/ for some guidance. There hinted that possibly I could use some mashup between data.table, parallel execution, and chunking the points. I would like to avoid the use of data.table if possible (however, the speed improvements identified on that site are fairly amazing).

One possibility I thought of was to go to a UTM style coordinate system and apply a transformation (rotation) which would place my polygon (always a rectangle, however arbitrarily sized and oriented) at the coordinate 0,0 and then just do a simple filter on xmin, xmax and ymin, ymax. However, while the search is very fast in that sense, the transformation would take some time.

Like I mentioned before, I wanted to get peoples opinion on which way I should go so as to not spin my wheels. Thanks for any info on what I can do with my specific use case and the current state of the tools.

• Platform
``````R: 3.5.2
RStudio: 1.3.66
SF: 0.9-2
``````
• Your polygons are rectangles without angle right? In this case do not test `point_in_polygon` but bounding box with simple coordinate comparisons. This is much faster. And if you polygon are not perfect rectangle use a bounding box test first.
– JRR
May 5, 2020 at 21:27
• Rectangle can have arbitrary orientation and position. But will always be rectangle. May 5, 2020 at 21:35
• Extracting community consensus isn't really what GIS SE is about. Achieving that goal might be better directed to Geographic Information Systems Chat (if Chat weren't sparsely visited). May 6, 2020 at 3:03

You can first use a simple bounding box comparison to quickly remove the points out of the bounding box without performing the costly `point_in_polygon`. Then you perform `point_in_polygon` with the subset only

``````pts = st_coordinates(bathy_data)
bbox = st_bbox(profile_segments[1,]\$buffer)

inbbox = pts[,1] >= bbox[1] & pts[,1] <= bbox[3] & pts[,2] >= bbox[2] & pts[,2] <= bbox[4]
sub = bathy_data[inbbox,]
points_to_include = unlist(sf::st_contains(profile_segments\$buffer[[1]], sub))
``````

This performs in 0.5 second on my computer against 5 seconds for your options

• Thanks for the suggestion. I tried it and it did in-fact work better for the specific case that I provided. When I tried it on the more general case (a slanted rectangle), there was no performance gain. I ended up with a different approach which I think works well in my case. I will provide that as a separate solution. May 6, 2020 at 1:56

Because my polygons are always rectangles, this problem becomes much easier than I expected it to be. My steps for solving the problem with high speed are the following:

1. Find the centroid of the bathymetric grid that I extracted from the database
2. Generate a local coordinate system based on that centroid
3. Re-project the bathymetric grid coordinates to the local grid
4. Re-project the profile_segment to the new local grid
5. Transform and then rotate the transformed bathymetric grid data
6. Use a simple box selection to determin which point indicies would have been inside the untransformed box
7. Rinse and repeat...

Here is the code that I used:

``````# The following is done once
bathy_bbox = sf::st_bbox(bathy_data)
bathy_center = c(mean(bathy_bbox[c(1,3)]), mean(bathy_bbox[c(2,4)]))
bathy_local_format_string = '+proj=aeqd +lat_0=%f +lon_0=%f +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
bathy_local_crs = sprintf(bathy_local_format_string,
bathy_center[2],
bathy_center[1])
pts = sf::sf_project(from=sf::st_crs(bathy_data), to=bathy_local_crs, sf::st_coordinates(bathy_data))

# The following is done for each segment

# This would of course iterate through a for statement.
# Just hard setting it here for the stack overflow post.
i = 1
segment_local = sf::sf_project(from=sf::st_crs(bathy_data), to=bathy_local_crs, profile_segments\$segments[[i]])
angle = (profile_segments\$bearing[i] * (pi/180))
M = matrix(c(cos(angle), -sin(angle), sin(angle), cos(angle)), ncol=2, byrow=TRUE)
pts_transformed = t(M %*% (t(pts) - segment_local[1,]))
points_to_include = which(
(pts_transformed[,1] >= -profile_segments\$buffer_width[i]) &
(pts_transformed[,1] <= profile_segments\$buffer_width[i]) &
(pts_transformed[,2] >= -profile_segments\$buffer_width[i]) &
(pts_transformed[,2] <= profile_segments\$length[i]+profile_segments\$buffer_width[i])
)

tmp_grid_points = bathy_data[points_to_include,]

``````

One of the great things about this solution is the expensive operation is taken out of the loop for multiple profiles. The individual profiles should be pretty quick.

• Instead of transforming all the points with a costly rotation, you can test if the points are inside the rectangle using the equation of its edges. It is like testing the bounding box but instead of `Y >= bbox[2] && Y <= bbox[4] ...` you have something like `Y >= a1*X+b1 && Y <= a2*X+b2 ...`
– JRR
May 6, 2020 at 10:12
• Interesting. I will look into that. Seems like that would be an improvement. However, trying to figure out the box edge slopes, y-intercepts, and test orientation might be fragile. But one just has to be careful. May 6, 2020 at 11:04
• Imo it is not fragile. If you are sure your polygons are real rectangles there is one limit case which is infinite slope (regular rectangle).
– JRR
May 6, 2020 at 12:10
• I am irking out the solution on this. As soon as I get it together I will update my answer. And yea, I will just have an if statement to handle the special case. May 6, 2020 at 12:13
• By the way, what does “lmo” mean? May 6, 2020 at 12:14