drawing a polygon from several lines
Very nice, Antonio!! Thanks for sharing your script. The only minor
suggestion I have to offer is that you can make your "cut isobath" portion
easier to read (at least IMO) by using `with` and `between` (either
dplyr::between or data.table::between). Below is how I would've done it.
limmin <- with(list(coords = coordinates(isomin)[[1]][[1]]),
subset(coords, between(coords[, 1], lonmin, lonmax) &
between(coords[, 2], latmin, latmax)))
Again, thanks for sharing your problem and solution. I found it
instructive!
Cordially,
Vijay.
On Fri, Aug 24, 2018 at 9:16 AM Antonio Silva <aolinto.lst at gmail.com> wrote:
Dear colleagues I finally could select the squares (cells) under the polygon. I extract the points from the selected isobaths that fell within lat / lon limits and with them I created a spatial polygon. Bellow I share the script I wrote. Considering that I don't know much about mapping in R, I would appreciate to hear suggestions to improve the code. Best regards. -- Ant?nio Olinto ?vila da Silva Instituto de Pesca (Fisheries Institute) S?o Paulo, Brasil =================== library(rgdal) library(rgeos) rm(list = ls()) # import shapes # isobath https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0 isobaths <- readOGR(".","isobath") isobaths <- spTransform(isobaths, CRS("+proj=longlat +datum=WGS84")) proj4string(isobaths) summary(isobaths) # create a spatial grid and polygons grd <- GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12)) squares <- as.SpatialPolygons.GridTopology(grd) proj4string(squares) <- CRS("+proj=longlat +datum=WGS84") summary(squares) IDs <- sapply(slot(squares, "polygons"), function(x) slot(x, "ID")) # set limits latmin <- -24.8 latmax <- -24.02 lonmin <- -47.2 lonmax <- -44.8 depmin <- -25 depmax <- -60 # bounding box mat.area <- matrix(c(lonmin,latmax,lonmax,latmax,lonmax,latmin,lonmin,latmin),ncol=2,byrow=T) spp.area <- SpatialPolygons(list(Polygons(list(Polygon(mat.area)) , ID='1'))) proj4string(spp.area) <- CRS("+proj=longlat +datum=WGS84") # plot polygons isolines and bounding box plot(squares, axes=T) #~ text(coordinates(squares),labels=IDs, cex=0.7) plot(isobaths,add=T) plot(spp.area,border="red",add=T) # select isobaths summary(isobaths) isomin <- subset(isobaths, ID %in% depmin) proj4string(isomin) <- CRS("+proj=longlat +datum=WGS84") isomax <- subset(isobaths, ID %in% depmax) proj4string(isomax) <- CRS("+proj=longlat +datum=WGS84") plot(isomin,add=T,col="deepskyblue3",lwd=2) plot(isomax,add=T,col="dodgerblue4",lwd=2) # cut isobaths limmin <- subset(coordinates(isomin)[[1]][[1]], coordinates(isomin)[[1]][[1]][,2]<=latmax & coordinates(isomin)[[1]][[1]][,2]>=latmin & coordinates(isomin)[[1]][[1]][,1]<=lonmax & coordinates(isomin)[[1]][[1]][,1]>=lonmin) limmax <- subset(coordinates(isomax)[[1]][[1]], coordinates(isomax)[[1]][[1]][,2]<=latmax & coordinates(isomax)[[1]][[1]][,2]>=latmin & coordinates(isomax)[[1]][[1]][,1]<=lonmax & coordinates(isomax)[[1]][[1]][,1]>=lonmin) points(limmin,col="chartreuse3") points(limmax,col="chartreuse4") # create the polygon spp.farea <- SpatialPolygons(list(Polygons(list(Polygon(rbind(limmin,limmax))),ID='1'))) proj4string(spp.farea) <- CRS("+proj=longlat +datum=WGS84") plot(spp.farea,col="chocolate3",lwd=2,add=T) # select the squares under the polygon fareasq <- squares[spp.farea,] fareace <- gCentroid(fareasq,byid=TRUE) # final plot plot(fareasq,axes=T) points(fareace,pch=10,col="darkgreen",cex=4) text(coordinates(squares), labels=IDs, cex=0.7) plot(isomin,add=T,col="deepskyblue3",lwd=2) plot(isomax,add=T,col="dodgerblue4",lwd=2) plot(spp.area,border="red",add=T,lty=2) [[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Vijay Lulla Assistant Professor, Dept. of Geography, IUPUI 425 University Blvd, CA-207C. Indianapolis, IN-46202 vlulla at iupui.edu ORCID: https://orcid.org/0000-0002-0823-2522 Webpage: http://vijaylulla.com [[alternative HTML version deleted]]