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

Question

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)
load('inputs.RData')
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

Answer 1

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

Answer 2

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.