Ir para Fase 4.

Configuramos o diretório de trabalho e carregamos as bibliotecas necessárias.

setwd("~/OneDrive/r-files/TCC")
library(contoureR)
## Loading required package: geometry
library(geosphere)
library(cluster)
library(sp)
library(rgdal)
## rgdal: version: 1.3-9, (SVN revision 794)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1
library(mapview)
library(RColorBrewer)
library(leaflet.extras)
## Loading required package: leaflet
library(cluster)
library(sp)
library(rgeos)
## rgeos version: 0.4-2, (SVN revision 581)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE

Carregamos o dataset.

load("alvaras.rda")
head(alvaras)

Verificamso os clusters contidos no dataset.

clusters_encontrados = sort(unique(alvaras$cluster))
length(clusters_encontrados)
## [1] 136
clusters_encontrados
##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [18]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34
##  [35]  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51
##  [52]  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68
##  [69]  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85
##  [86]  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 101 102
## [103] 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
## [120] 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136

Criamos um dataframe “poly” e o povoamos com os dados que desejamos que estejam contidos no data frame de polígono espacial.

load("exemplars.rda")
poly = data.frame()
for (i in clusters_encontrados){
  temp = alvaras[alvaras$cluster == i,  ]
  ch1 = convexHullAM_Indexes(temp[,1],temp[,2], includeColinear=FALSE,zeroBased = FALSE)
  poligono = temp[ch1, 1:2 ]
  area <- geosphere::areaPolygon(x = poligono)
  medoide_lat=pam(temp, 1)$medoids[,1]
  medoide_lon=pam(temp, 1)$medoids[,2]
  pol = temp
  coordinates(pol) = ~V1+V2
  centr_lat=gCentroid(pol, byid=FALSE)$x
  centr_lon=gCentroid(pol, byid=FALSE)$y
  centr_id=pam(temp, 1)$medoids[,6]
  exemplars_lat = exemplars[exemplars$Cluster == i, 1]
  exemplars_lon = exemplars[exemplars$Cluster == i, 2]
  for (ii in ch1) {
    polying = temp[ii,]
    polying$area = area/1000000 #km2
    polying$ocorrencias = length(temp$cluster)
    polying$centroide_lat = centr_lat
    polying$centroide_lon = centr_lon
    polying$medoide_lat = medoide_lat
    polying$medoide_lon = medoide_lon
    polying$exemplars_lat = exemplars_lat
    polying$exemplars_lon = exemplars_lon
    polying$centr_id = centr_id
    poly = rbind(poly, polying)
  }
}
head(poly)
tail(poly)
dim(poly)
## [1] 1837   15
dados = poly
names(dados)
##  [1] "V1"            "V2"            "V3"            "V4"           
##  [5] "V5"            "cluster"       "area"          "ocorrencias"  
##  [9] "centroide_lat" "centroide_lon" "medoide_lat"   "medoide_lon"  
## [13] "exemplars_lat" "exemplars_lon" "centr_id"

Renomeamos os atributos do data frame e criamos novas variáveis necessárias.

names(dados) = c("lat", "lon", "logr", "group", "segmento", "box_id", "area", "ocorrencias", "centroide_lat", "centroide_lon", "medoide_lat", "medoide_lon", "exemplars_lat", "exemplars_lon","centr_id")
dados$group = dados$box_id
dados$cluster = dados$box_id
dados$id = (dados$ocorrencias / dados$area)
dados$box_id = dados$id
head(dados)
rm(poly)

Criamos o data frame temporário e enxuto.

dadostemp = dados
dadostemp$lat = NULL; dadostemp$lon = NULL; dadostemp$logr = NULL; dadostemp$segmento = NULL; dadostemp$box_id = NULL
dadostemp = unique(dadostemp)
dim(dadostemp)
## [1] 136  12
head(dadostemp)
save(dadostemp, file = "clusters_centroids.rda")

Transformamos o data set original em um objeto da classe ‘SpatialPointsDataFrame’.

coordinates(dados)=c("lat","lon")
df = dados
class(df)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
str(df)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  1837 obs. of  15 variables:
##   .. ..$ logr         : Factor w/ 60123 levels "90020060","90020100",..: 19679 30342 34537 50317 30221 19495 37247 46342 34789 33042 ...
##   .. ..$ group        : num [1:1837] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ segmento     : Factor w/ 2 levels "Comercio","Servico": 2 1 2 2 2 2 2 2 1 2 ...
##   .. ..$ box_id       : num [1:1837] 3535 3535 3535 3535 3535 ...
##   .. ..$ area         : num [1:1837] 0.65 0.65 0.65 0.65 0.65 ...
##   .. ..$ ocorrencias  : int [1:1837] 2299 2299 2299 2299 2299 2299 2299 2299 2299 2299 ...
##   .. ..$ centroide_lat: num [1:1837] -51.2 -51.2 -51.2 -51.2 -51.2 ...
##   .. ..$ centroide_lon: num [1:1837] -30 -30 -30 -30 -30 ...
##   .. ..$ medoide_lat  : num [1:1837] -51.2 -51.2 -51.2 -51.2 -51.2 ...
##   .. ..$ medoide_lon  : num [1:1837] -30 -30 -30 -30 -30 ...
##   .. ..$ exemplars_lat: num [1:1837] -51.2 -51.2 -51.2 -51.2 -51.2 ...
##   .. ..$ exemplars_lon: num [1:1837] -30 -30 -30 -30 -30 ...
##   .. ..$ centr_id     : num [1:1837] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ cluster      : num [1:1837] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ id           : num [1:1837] 3535 3535 3535 3535 3535 ...
##   ..@ coords.nrs : int [1:2] 1 2
##   ..@ coords     : num [1:1837, 1:2] -51.2 -51.2 -51.2 -51.2 -51.2 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:1837] "86026" "124694" "78386" "12638" ...
##   .. .. ..$ : chr [1:2] "lat" "lon"
##   ..@ bbox       : num [1:2, 1:2] -51.3 -30.2 -51.1 -30
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lat" "lon"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

ReCriamos uma função que fará o bind dos objetos ‘SpatialPointsDataFrame’ e ‘data.frame’, criando um objeto do tipo ‘SpatialPolygonsDataFrame’.

data <- data.frame(box_id=unique(df$box_id),row.names=unique(df$id))
head(data)
data = cbind(data, dadostemp)
head(data)
dim(data)
## [1] 136  13
points2polygons <- function(df,data) {
  get.grpPoly <- function(group,ID,df) {
    Polygon(coordinates(df[df$id==ID & df$group==group,]))
  }
  get.spPoly  <- function(ID,df) {
    Polygons(lapply(unique(df[df$id==ID,]$group),get.grpPoly,ID,df),ID)
  }
  spPolygons  <- SpatialPolygons(lapply(unique(df$id),get.spPoly,df))
  SpatialPolygonsDataFrame(spPolygons,match.ID=T,data=data)
}

Criamos o objeto ‘SpatialPolygonsDataFrame’.

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -51.18699, -51.1777, -30.02746, -30.01799  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 13
## names       :           box_id, group,              area, ocorrencias, centroide_lat, centroide_lon, medoide_lat, medoide_lon, exemplars_lat, exemplars_lon, centr_id, cluster,               id 
## value       : 3534.73117030931,     1, 0.650403068643782,        2299,  -51.18344652,  -30.02297009, -51.1789439, -30.0203054,   -51.1834357,   -30.0231583,        1,       1, 3534.73117030931
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
A caption

A caption

Salvamos o objeto ‘SpatialPolygonsDataFrame’ em um banco de dados NoSql e não tabular e sim do tipo GeoJson.

library(rgdal)
rgdal::writeOGR(obj = spDF,
                dsn = "alvaraspoly.json",
                layer = "alvaraspoly",
                driver = "GeoJSON",
                overwrite_layer = TRUE)

Carregamos o arquivo GeoJson salvo previamente.

Alvaras <- geojsonio::geojson_read("alvaraspoly.json", what = "sp")
head(Alvaras)
## class       : SpatialPolygonsDataFrame 
## features    : 6 
## extent      : -51.26295, -51.05223, -30.24431, -29.99994  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       :           box_id, group,             area, ocorrencias, centroide_lat, centroide_lon, medoide_lat, medoide_lon, exemplars_lat, exemplars_lon, centr_id, cluster,               id 
## min values  : 1.86752510668252,    26, 2.20888327306396,          30,  -51.25912282,  -30.21957912, -51.2594808,  -30.225552,   -51.2591895,   -30.2229995,       26,      26, 1.86752510668252 
## max values  : 25.8048945795741,   134, 16.6517856537877,         202,  -51.08494452,  -30.01609689, -51.0742668, -30.0155273,   -51.0805274,   -30.0168615,      134,     134, 25.8048945795741

Plotamos os dados com a biblioteca mapview.
Mapa de densidade de alvarás.

Mapa de área total.

Mapa de números de ocorrências de alvarás.

## [1] 136  12
##                                                                            
## 1 function (x, ..., degree = 1, coefs = NULL, raw = FALSE, simple = FALSE) 
## 2 {                                                                        
## 3     dots <- list(...)                                                    
## 4     if (nd <- length(dots)) {                                            
## 5         dots_deg <- nd == 1L && length(dots[[1L]]) == 1L                 
## 6         if (dots_deg)

Ir para Fase 6.

.