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

[Antonio Silva]

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]]