CODEWITHSUNDEEP
rclip
Raj
Raj

Reputation: 75

Selection and extraction of points (shapefile) within a polygon in R

I have two shape-files; one is a point file (some information for world) and the other is a 21-countries shape-file.

I need to extract the points which falls with in a country. I have to repeat this step in 21 times in QGIS or ArcGIS.

Is there any way to do this in R and if is in batch processing, that would be really great as I have to repeat this for other 4 data sets as well. Thanks in advance

Upvotes: 2

Views: 7410

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47601

Here is an example polygons shapefile and how you can read it into R.

library(raster)
# example polygons filename
f <- system.file("external/lux.shp", package="raster")
pol <- shapefile(f)

I do not have a matching points file, so I generate some points for this example

pts <- spsample(pol, 5, type="regular")
crs(pts) <- crs(pol)

To see what area the points fall in you can use over or extract. They both return the values in the order of the points, but extract also adds the (sequential) "point id" --- this is primarily for when you have overlapping polygons, such that a point might fall in two polygons.

over(pts, pol)
#  ID_1       NAME_1 ID_2           NAME_2 AREA
#1    3   Luxembourg    9 Esch-sur-Alzette  251
#2    2 Grevenmacher    7           Remich  129
#3    3   Luxembourg   11           Mersch  233
#4    2 Grevenmacher    6       Echternach  188
#5    1     Diekirch    1         Clervaux  312

extract(pol, pts)
#  point.ID poly.ID ID_1       NAME_1 ID_2           NAME_2 AREA
#1        1      10    3   Luxembourg    9 Esch-sur-Alzette  251
#2        2       7    2 Grevenmacher    7           Remich  129
#3        3      12    3   Luxembourg   11           Mersch  233
#4        4       6    2 Grevenmacher    6       Echternach  188
#5        5       1    1     Diekirch    1         Clervaux  312

If you just want the intersection, say of the points that fall in "Grenmacher", like you may want to do with your countries, you can 

g <- pol[pol$NAME_1 == "Grevenmacher", ]
i <- intersect(pts, g)
i
#class       : SpatialPointsDataFrame 
#features    : 2 
#extent      : 6.27111, 6.27111, 49.53585, 49.7887  (xmin, xmax, ymin, ymax)
#crs         : +proj=longlat +datum=WGS84 +no_defs 
#variables   : 5
#names       : ID_1,       NAME_1, ID_2,     NAME_2, AREA 
#min values  :    2, Grevenmacher,    6, Echternach,  129 
#max values  :    2, Grevenmacher,    7,     Remich,  188

And perhaps write that back to disk like this 

shapefile(i, "test.shp")

Upvotes: 3

horseoftheyear
horseoftheyear

Reputation: 915

So as an example how to use over() consider the following simulated data.

library(raster)
library(sp)

#generate point distribution | with (x,y) coordinates
set.seed(42);pts <- SpatialPoints(data.frame(x = runif(30, 1, 5), y = runif(30, 1, 5)))

#create polygons
df<-data.frame(X = c(rep(1,5),rep(2,5),rep(3,5),rep(4,5),rep(5,5)),
               Y = c(rep(seq(1,5,1),5)))
df$cell<-1:nrow(df) #polygon identifier (for later)

#make into spatial object (probably better way to do this)
coordinates(df) <- ~X+Y 
rast <- raster(extent(df), ncol = 5, nrow = 5)
grid<-rasterize(df, rast, df$cell, FUN = max)
grid<-rasterToPolygons(grid) #create polygons

#plot to check
plot(grid); points(pts)

#and extract
pointsinpolygon = over(SpatialPolygons(grid@polygons), SpatialPoints(pts))

The last bit you can also do as

df$pointsinpolygon <- over(SpatialPolygons(grid@polygons), SpatialPoints(pts))

to write the results directly to the dataframe.

Upvotes: 5

Related Questions