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