Selects fishing area squares per trip and calculates their centroids

fishing_area.R

# ====================================================================
# Selects fishing area squares per trip and calculates their centroids
# author: Antônio Olinto Ávila da Silva
# creation: 2018-08-24
# last edition: 2018-08-24
# ====================================================================

# prepare workspace
library(rgdal)
library(rgeos)
rm(list = ls())

# import shapes
isobaths <- readOGR("/home/antonio/Documentos/SIG/Ambiental/Batimetria/bathimetry","isobath")
isobaths <- spTransform(isobaths, CRS("+proj=longlat +datum=WGS84"))
proj4string(isobaths)
summary(isobaths)

squares <- readOGR("/home/antonio/Documentos/SIG/BlocosPesca/Quadrados_SSE","QDLon05Lat05")
squares <- spTransform(squares, CRS("+proj=longlat +datum=WGS84"))
proj4string(squares)
summary(squares)

# 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 for the target area
spp.tarea <- SpatialPolygons(list(Polygons(list(Polygon(rbind(limmin,limmax))),ID='1')))
proj4string(spp.tarea) <- CRS("+proj=longlat +datum=WGS84")

plot(spp.tarea,col="chocolate3",lwd=2,add=T)

# select the squares under the polygon and calculate the centroids
tareasq <- squares[spp.tarea,]
tareace <- gCentroid(tareasq,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)