# drawing a polygon from several lines

4 messages
Open this post in threaded view
|

## drawing a polygon from several lines

 Hello all, I want to select the cells of a grid that in a range of latitudes and depths. Latitudes limits are indicated by two lines and the depths by two isobaths from a shapefile. The shapefile can be downloaded at https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0Below are the steps I've already taken. I could not find a way to create a polygon from the lines of depth and latitude. With this polygon I could easily select the overlapped cells. I would really appreciate any help. Best regards ---- library(rgeos) library(sp) rm(list = ls()) # import shape isobaths <- readOGR(".","isobath") isobaths <- spTransform(bath, 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)) polys <- as.SpatialPolygons.GridTopology(grd) proj4string(polys) <- CRS("+proj=longlat +datum=WGS84") summary(polys) IDs <- sapply(slot(polys, "polygons"), function(x) slot(x, "ID")) # create spatial lines # southern limit x <- c(-47.5,-45.5) y <- c(-24.8,-24.8) ls <- SpatialLines(list(Lines(Line(cbind(x,y)), ID="b")),proj4string=CRS("+proj=longlat +datum=WGS84")) # northern limit x <- c(-46.2,-44.5) y <- c(-24.02,-24.02) ln <- SpatialLines(list(Lines(Line(cbind(x,y)), ID="c")),proj4string=CRS("+proj=longlat +datum=WGS84")) # plot isolines polygons and lines plot(polys) axis(1);axis(4) text(coordinates(polys), labels=IDs, cex=0.7) plot(isobaths,add=T) plot(ls,add=T,col="blue",lwd=2) plot(ln,add=T,col="orange",lwd=2) # select isobaths summary(isobaths) dmin <- subset(isobaths, ID %in% "-25") proj4string(dmin) <- CRS("+proj=longlat +datum=WGS84") dmax <- subset(isobaths, ID %in% "-60") proj4string(dmax) <- CRS("+proj=longlat +datum=WGS84") plot(dmin,add=T,col="red",lwd=2) plot(dmax,add=T,col="red",lwd=2) # from here I'd like to have a polygon to select the cells. na.omit(over(polys,dmin)) -- Antonio 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: drawing a polygon from several lines

 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=0isobaths <- 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 [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo