仅仅因为你只有一把锤子，并不意味着每个问题都是钉子。
将您的数据加载到PostGIS中，为多边形构建一个空间索引，然后进行单个SQL空间叠加。将结果导出回R。
顺便说一下，说shapefile的大小是8GB并不是非常有用的信息。Shapefile至少由三个文件组成，.shp是几何图形，.dbf是数据库，.shx连接两者。如果您的.dbf是8GB，那么您可以通过替换它来轻松读取形状本身。即使.shp是8GB，它也可能只包含三个多边形，这种情况下简化它们可能很容易。您有多少个多边形，.shp部分的大小是多少？
我很想了解这个问题并想知道你在这方面是否有所进展。既然你问了这个问题，我想你的电脑硬件和软件应该已经得到了提升，可以用来进行相对简单的操作，尽管处理百万个点可能需要较长时间。你可以尝试使用以下方法：
# Load relevant libraries
library(sp)
library(maptools)
library(spatstat)
# Read shapefile. Hopefully you have a .prj file with your .shp file
# otherwise you need to set the proj4string argument. Don't inlcude
# the .shp extension in the filename. I also assume that this will
# create a SpatialPolygonsDataFrame with the "Severity" attribute
# attached (from your .dbf file).
myshapefile <- readShapePoly("myshapefile_without_extension", proj4string=CRS("+proj=latlong +datum=WGS84"))
# Read your location data in. Here I assume your data has two columns X and Y giving locations
# Remeber that your points need to be in the same projection as your shapefile. If they aren't
# you should look into using spTransform() on your shapefile first.
mylocs.df <- read.table(mypoints.csv, sep=",", h=TRUE)
# Coerce X and Y coordinates to a spatial object and set projection to be the same as
# your shapefile (WARNING: only do this if you know your points and shapefile are in
# the same format).
mylocs.sp <- SpatialPoints(cbind(mylocs.df$X,mylocs.df$Y), proj4string=CRS(proj4string(myshapefile))
# Use over() to return a dataframe equal to nrows of your mylocs.df
# with each row corresponding to a point with the attributes from the
# poylgon in which it fell.
severity.df <- over(mylocs.sp, myshapefile)
希望这个框架能够满足您的需求。但是，您现在可用的计算机/内存是否可以完成此操作是另一回事！