subsetting a spatial polygons

4 messages
Open this post in threaded view
|

subsetting a spatial polygons

 Dear list users I have a SpatialPolygons with several squares. How to subset it to have only the squares between given latitudes and longitudes? In the example library(sp) library(raster) grd <- GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12)) polys <- as.SpatialPolygons.GridTopology(grd) proj4string(polys) <- CRS("+proj=longlat +datum=WGS84") plot(polys,axes=T) How to select only the squares, let's say, between 24-25°S and 45-46°W? The farthest I went was: e <- extent(-45.9,-45.1,-24.9,-24.1) # which is not very elegant mask <- crop(polys,e) polys2 <- polys[mask,] plot(polys2,add=T,col="green") Thanks a lot. Best regards -- Antônio Olinto Ávila da Silva Instituto de Pesca (Fisheries Institute) São Paulo, Brasil         [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Open this post in threaded view
|

Re: subsetting a spatial polygons

 Antonio, I am not sure why you think that your solution is not very elegant. In case you want to have more visual control over the subsetting, you could try mapedit: library(mapedit) myselection = selectFeatures(polys, mode = "draw") which will let you draw a e.g. rectangle and only return those features that intersect it. Best Tim On 09/12/2018 04:18 AM, Antonio Silva wrote: > Dear list users > > I have a SpatialPolygons with several squares. How to subset it to have > only the squares between given latitudes and longitudes? > > In the example > > library(sp) > library(raster) > grd <- > GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12)) > polys <- as.SpatialPolygons.GridTopology(grd) > proj4string(polys) <- CRS("+proj=longlat +datum=WGS84") > plot(polys,axes=T) > > How to select only the squares, let's say, between 24-25°S and 45-46°W? > > The farthest I went was: > > e <- extent(-45.9,-45.1,-24.9,-24.1) # which is not very elegant > mask <- crop(polys,e) > polys2 <- polys[mask,] > plot(polys2,add=T,col="green") > > Thanks a lot. Best regards > _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo