[R-sig-Geo] drawing a polygon from several lines

[Vijay Lulla]

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 using 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 using 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 using iupui.edu
ORCID: https://orcid.org/0000-0002-0823-2522
Webpage: http://vijaylulla.com

	[[alternative HTML version deleted]]