Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have a Shapefile of points with certain information. I have another Shapefile with polygons corresponding to the states of the country I am working with (in this case Mexico). Since I only want to work specifically with 3 states, I wanted to know if there is a way to get rid of the points that dont fall into these 3 states.

share|improve this question
add comment

2 Answers

up vote 2 down vote accepted

You can do this type of spatial overlay using the "over" function in the sp library. This library is automatically loaded with libraries like maptools, rgdal, raster and spdep. This is because it defines the primary spatial object classes in R. If you are performing spatial analysis in R, this is a very good library to become familiar with. Almost all the spatial libraries, excepting spatstat, in R use the object classes defined by sp. For more complex geometry operations there is the rgeos library that is a port to the GEOS JTS Topology Suite. This library allows operations like buffering, union and intersect. Here is the example from the over help for what you after.

require(sp)

r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409, 
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676, 
332618, 332413, 332349))

r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437, 
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683, 
331133, 331623, 332152, 332357, 332373))

r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875, 
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004, 
329783, 329665, 329720, 329933, 330478, 331062, 331086))

r4 = cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791))

sr1=Polygons(list(Polygon(r1)),"r1")
sr2=Polygons(list(Polygon(r2)),"r2")
sr3=Polygons(list(Polygon(r3)),"r3")
sr4=Polygons(list(Polygon(r4)),"r4")
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4))
srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2), row.names=c("r1","r2","r3","r4")))

data(meuse)
coordinates(meuse) = ~x+y

plot(meuse)
polygon(r1)
polygon(r2)
polygon(r3)
polygon(r4)
# retrieve mean heavy metal concentrations per polygon:
over(sr, meuse[,1:4], fn = mean)    
share|improve this answer
    
Got it, thank you! –  JEquihua Mar 11 '13 at 19:10
add comment

You don't need to crop anything.

  • First make a subset containing your 3 state polygons (if you haven't done already or if there are more than 3 states)
  • Then select all points that fall within those polygons (Query by location plugin in QGis. Make a new selection containing all points that are within your polygons).
  • Rightclick on your point layer (points still selected) and execute the function "Save selection as", which will just save your selected points somewhere.
share|improve this answer
    
This is not bad at all. Thank you. But I wanted a way to do it directly in R. Is this too cumbersome? –  JEquihua Mar 11 '13 at 18:07
    
This is of course possible but why did you then tag your question with QGis ? I don't have time to write your code, but have a look at the package sp and the functions over and subset. Manual pages should provide you with everything you need –  Curlew Mar 11 '13 at 18:32
    
Wow, that actually is kind of strange. I don't know why I did that. Thank you, I'll look into the Manual pages! –  JEquihua Mar 11 '13 at 19:02
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.