Cuadriculas regulares de resolución fija dentro de un polígono delimitado. Regular quadrats with a fixed resolution in a delimited polygon

quadrats.R

#GENERACION DE CUADRICULAS REGULARES (QUADRATS). ES NECESARIO APORTAR RESOLUCION Y POLIGONO DELIMITADOR
#EJEMPLO APLICADO AL CAMPO DE LA UNIVERSIDAD AUTÓNOMA DE SANTO DOMINGO (UASD), CON CUADRICULAS DE 25X25 M Y 50X50 M

#CARGA DE PAQUETES
library(raster) #OBJETOS RASTER
library(ggmap) #MAPAS ESTETICOS. CARGA DE MANERA AUTOMATICA ggplot2. IMPORTANTE; INSTALAR VERSION MAS RECIENTE MEDIANTE install.packages('ggmap')
library(rgeos) #CALCULO DE AREAS
library(OpenStreetMap) #IMAGEN DE SATELITE
library(sp) #OBJESTOS ESPACIALES
library(maptools) #EXPORTACION
library(rgdal) #BIBLIOTECA GDAL
library(plotKML) #GENERAR KML CON NOMBRES

#DIRECTORIO DE TRABAJO
setwd(tempdir()) #FIJAR EL DIRECTORIO DE TRABAJO
getwd() #MUESTRA DONDE ESTA EL DIRECTORIO DE TRABAJO

#DESCARGA POLIGONO. FUNCION shapefile() DEL PAQUETE RASTER, fortify(), get_map(), geom_polygon() DEL PAQUETE ggplot2. bbox() DEL PAQUETE SP
temp <- tempfile() #ASIGNAR EL NOMBRE DE UN ARCHIVO TEMPORAL AL OBJETO temp
download.file("http://www.geografiafisica.org/sem_2016_02/geo131/gis/UASD/perimetro/perimetro.zip",temp) #DESCARGA EL ARCHIVO ZIP Y LO NOMBRA COMO TEMPORAL
unzip(temp) #DESCOMPRIME EN EL DIRECTORIO DE TRABAJO
unlink(temp) #DESVINCULA DEL ARCHIVO TEMPORAL Y LO BORRA
uasd <- shapefile('poligono_uasd.shp')
uasdutm <- shapefile('poligono_uasd_epsg_32619.shp')

#FUNCION DE GENERACION DE CUADRICULA
cuadriculas <- function(resolucion,encuadre){
  require(sp)
  require(raster)
  xmn <- floor(bbox(encuadre)[1,1]/100)*100
  xmx <- ceiling(bbox(encuadre)[1,2]/100)*100
  ymn <- floor(bbox(encuadre)[2,1]/100)*100
  ymx <- ceiling(bbox(encuadre)[2,2]/100)*100
  rutm <- raster(xmn=xmn,xmx=xmx,ymn=ymn,ymx=ymx, resolution = resolucion, crs = CRS("+init=epsg:32619"))
  rutm[] <- 1:ncell(rutm)
  rutmm <- mask(crop(rutm,extent(encuadre)),encuadre)
  rutmm[!is.na(rutmm)] <- 1:length(rutmm[!is.na(rutmm)])
  rgll <- projectRaster(rutmm, crs="+proj=longlat +datum=WGS84", method='ngb')
  return(list(utm=rutmm, latlon=rgll))
}

#BUCLE for() PARA EJECUTAR FUNCION CON CUADRICULAS DE DISTINTOS TAMANOS
#APLICACION A RESOLUCIONES 25 Y 50 M
i <- NULL; res <- c(25,50)
for(i in res){nom <- paste0('c',i,'m'); assign(nom,cuadriculas(i,uasdutm))}
#EXPORTAR CUADRICULAS A FORMATOS SHAPEFILE Y KML. ESTOS ULTIMOS PUEDEN VISUALIZARSE CON GOOGLEEARTH O SUBIRSE A GOOGLEMAPS
lapply(res,function(x){
  y <- cuadriculas(x,uasdutm);
  writeOGR(as(y[[1]],'SpatialPolygonsDataFrame'), dsn=getwd(), layer = paste0('c',x,'m'), driver = 'ESRI Shapefile')
  writeOGR(as(y[[2]],'SpatialPolygonsDataFrame'), dsn=paste0(getwd(),paste0('/c',x,'m'),'.kml'), layer = paste0('c',x,'m'), driver = 'KML')
  plotKML(as(y[[1]],'SpatialPointsDataFrame'), paste0('c',x,'m-rot'), points_names=as(y[[1]],'SpatialPointsDataFrame')@data$layer, shape='')
  }
)


#OPCIONALES. CALCULO DE AREA Y PERIMETRO. VISUALIZACION SOBRE IMAGEN DE SATELITE DENTRO DE R

#AREA, PERIMETRO. FUNCIONES gArea() Y gLength() DEL PAQUETE rgeos
area <- gArea(uasdutm); area
perimetro <- gLength(uasdutm); perimetro

#IMAGEN DE BING, POLIGONO UASD, CUADRICULAS
m <- openmap(c(bbox(c50m[[2]])[2,2],bbox(c50m[[2]])[1,1]),c(bbox(c50m[[2]])[2,1],bbox(c50m[[2]])[1,2]), zoom=18, 'bing')
mutm <- openproj(m,proj4string(uasdutm))
dev.new()
plot(mutm)
plot(uasdutm,add=T, lwd = 3, border = 'green')
plot(as(c50m$utm,'SpatialPixels'), lwd = 3, col = rgb(1,1,1,0.5), add=T)
text(c50m$utm, col = rgb(0,1,0,0.5))

#VISUALIZAR MAPA CON ggmap
uasdf <- fortify(uasd)
lonc <- mean(bbox(uasd)[1,])
latc <- mean(bbox(uasd)[2,])
ptc <- c(lonc,latc)
muasd <- get_map(ptc, zoom = 16, maptype = 'satellite')
dev.new();ggmap(muasd) +
  geom_polygon(data = uasdf, aes(long, lat, group = group), fill='green', size=1, color='black',alpha=0.2) +
  theme(legend.position = 'none')