Cria grade de polígonos de tamanho personalizados e exporta como shapefile

grade.R

# ================================================================================
# Cria grade de polígonos de tamanho personalizados e exporta como shapefile
# autor: Antônio Olinto Ávila da Silva
# criação: 2018-08-16
# última edição: 2018-08-16
# ================================================================================

# ------------------------
# prepara área de trabalho
# ------------------------
#~ no Linux: sudo apt update && sudo apt install libgdal-dev libproj-dev
library(rgdal)
library(sp)
setwd("...") # caminho da pasta de trabalho
rm(list = ls())

# ---------------------------------
# parâmetros customizáveis da grade
# ---------------------------------

# faixa de longitude
xmin <- -50
xmax <- -41

# faixa de latitude
ymin <- -29
ymax <- -22

# intervalo em milhas náuticas
xint <- 5
yint <- 10

# --------------------
# variáveis calculadas
# --------------------

# número de intervalos
xnin <- 60/xint
ynin <- 60/yint

# intervalo em graus
xstep <- xint/60
ystep <- yint/60

# letras do código
xlet <- LETTERS[seq(1:(60/xint))]
ylet <- LETTERS[seq(1:(60/yint))]

# sequência de long/lat para quadrados
x <- seq(xmin+xint/60,xmax,xstep)
y <- seq(ymin+yint/60,ymax,ystep)

# número de quadrados por eixo
xn <- length(x)
yn <- length(y)

# ---------------------------------
# criação dos polígonos e da camada 
# ---------------------------------

# tabela com a lista de pontos
grid.xy <- expand.grid(x[1:xn],y[1:yn])
grid.xy$id <- paste(
abs(as.integer(grid.xy[,1])),
xlet[round(abs(grid.xy[,1]-as.integer(grid.xy[,1]))*xnin+1,0)],
abs(as.integer(grid.xy[,2])),
ylet[round(abs(grid.xy[,2]-as.integer(grid.xy[,2]))*ynin+1,0)],
rep(sprintf("-WS%02d%02d",xint,yint,sep=""),xn*yn),sep="")
names(grid.xy) <- c("lon","lat","id")
#~ head(grid.xy)
#~ plot(grid.xy[,1],grid.xy[,2])

# tabela com a grade
vertices <- data.frame(rep(grid.xy[,1],each=4),rep(grid.xy[,2],each=4),rep(1:4,nrow(grid.xy)),rep(grid.xy[,3],each=4))
names(vertices) <- c("lon","lat","ver","id")
#~ head(vertices)
for (i in seq(1,nrow(vertices),4)) {
	vertices[i+1,1] <- vertices[i,1]-xstep
	vertices[i+2,1] <- vertices[i,1]-xstep
	vertices[i+2,2] <- vertices[i,2]-ystep
	vertices[i+3,2] <- vertices[i,2]-ystep
}
#~ head(vertices)

# polígonos
polys <- SpatialPolygons(lapply(unique(vertices$id),
                         function(x) {
                           Polygons(list(Polygon(vertices[vertices$id==x,1:2])), ID=x)
                         }),proj4string=CRS("+proj=longlat +datum=WGS84"))
#~ summary(polys)
#~ plot(polys)

# cria dataframe com ids
polys.df <- data.frame(ID=grid.xy[,3])
#~ head(polys.df)

# extrai o id dos polígonos
pid <- sapply(slot(polys,"polygons"),function(x) slot(x, "ID"))

# cria dataframe com ids como nome de linha e coluna 
p.df <- data.frame( ID=pid,row.names = pid)

# cria spatial polygon dataframe
p <- SpatialPolygonsDataFrame(polys, p.df)
#~ class(p) 

# exporta como shapefile
writeOGR(p,dsn='.',layer=sprintf("QDLon%02dLat%02d",xint,yint),overwrite_layer=TRUE,driver="ESRI Shapefile")