processamento_de_mapas.R

rm(list=ls())
graphics.off()
if (.Platform$OS.type == 'windows') windows(record = TRUE)

library(dplyr, quietly=TRUE)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(maptools, quietly=TRUE)

## Warning: package 'maptools' was built under R version 3.2.5

## Checking rgeos availability: TRUE

library(spdep, quietly=TRUE)

## Warning: package 'spdep' was built under R version 3.2.5

## get .shp, .shx, and .dbf files for Brazilian 'microregion' geography

need.map = !('BR.shp' %in% dir() )
if (need.map) {
  for (suffix in c('.shp', '.shx', '.dbf')) {
    url = paste0('http://schmert.net/data/brazil_microregions_1960_2000', suffix)
    local.filename = paste0('BR',suffix)
    download.file(url, destfile = local.filename, mode='wb')
  } # for suffix
} # if need map  

map = readShapePoly('BR')  # this is one of many ways to read shapefiles

# The map@data object is a dataframe that contains the basic info about the polygons.
#  Add a state variable to map@data.

selected.states = 23:29

map@data = map@data %>%
             mutate(state = KREISE6020 %/% 10000)

head(map@data)

##   KREISE6020 NUMBERMCA NUMBERMUN   X_COORD   Y_COORD state
## 1     110006         1        52 -63.28978 -10.83083    11
## 2     120001         1         5 -73.00868  -8.30643    12
## 3     120002         2         3 -70.95150  -8.86857    12
## 4     120004         1        14 -68.86855  -9.86200    12
## 5     130001         4         8 -63.68725  -0.42970    13
## 6     130002         2         2 -66.90512  -1.87192    13

## make a second map that contains only the NE microregion
NEmap = map[ map$state %in% selected.states, ]

## plot the map, then add a layer with the NE regions filled in
plot(map, border='grey')
plot(NEmap, border='red',col='skyblue', add=TRUE)
title('Brazilian microregions (NE highlighted)')

## now plot NE only
plot(NEmap, border='red',col='skyblue')
title('Northeast Brazilian microregions')

## construct a list of adjacent polygons
neighbor.list = poly2nb(NEmap)
head(neighbor.list)  # 1 is adj to {2,3,5,8}; 2 is adj to {1,3}, etc.

## [[1]]
## [1] 2 3 5 8
## 
## [[2]]
## [1] 1 3 5 6
## 
## [[3]]
## [1] 1 2 4 5
## 
## [[4]]
## [1] 3 5
## 
## [[5]]
## [1]  1  2  3  4  6  7  8 10 12
## 
## [[6]]
## [1]  2  5  7 17

## illustrate the adjacencies that we just calculated on the map
plot(NEmap, border='white',col='skyblue')

cent           = coordinates(NEmap)
colnames(cent) = c('long','lat')
head(cent)

##         long       lat
## 81 -40.59787 -3.097523
## 82 -40.99263 -3.884859
## 83 -40.70692 -3.506001
## 84 -40.49399 -3.590043
## 85 -40.21686 -3.741580
## 86 -40.72395 -4.429659

points(cent, pch=16, col='darkblue', cex=0.8)

for (this.region in seq(NEmap)) {
  for (n in neighbor.list[[this.region]]) {
    segments( cent[this.region,'long'], cent[this.region,'lat'],
              cent[n,'long'], cent[n,'lat'], col='red')
  }  
}