Descarga, lectura y representación de las secciones censales del censo de 2001. Andalucía

LeeSHPFileSeccCen2001.R

library(rgdal)
require(ggmap)
rm(list = ls())

##  Descarga las secciones censales del censo de 2001, previa preparacion
##  De la pagina:
## 'http://www.juntadeandalucia.es/institutodeestadisticaycartografia/clientedescarga/
##  Hay que seleccionar manualmente la capa : "Secciones Censales noviembre de de 2001 (censo)"  y 
##  y marcar  "extensión completa" esperar unos minutos, que se genere la petición y descargarla. 
##  O mejor copiar el nombre del fichero ZIP (con shapefiles) y  modificar las dos linead dabajo de esta
if (! dir.exists('SC2001A') ) dir.create('SC2001A')
download.file('http://www.juntadeandalucia.es/institutodeestadisticaycartografia/clientedescarga/proxyDownload?id=1487591252579CLIPED_FEATURESresult-15051480-6154-4ce4-b838-cc71948adddb&status=1&format=ZIP', 
                  destfile = 'SC2001A/Secc2001And.zip')

unzip('SC2001A/Secc2001And.zip', exdir = './SC2001A/')  # descomprime
(unzip('SC2001A/Secc2001And.zip', list = T) -> nn)  
nn <- paste0('./SC2001A/',nn$Name) ; 
gsub('shp[0-9]+','Secc2001And',nn)->nn2 
file.rename(nn,nn2) # Cambia Nombre aleatorio
file.remove('SC2001A/Secc2001And.zip')
dir(path = 'SC2001A')
rm(nn,nn2)
##
scA  <- readOGR('SC2001A/','Secc2001And')
summary(scA)
## No ha leido el DATUM ??? hay que ponerlo a mano
##------------------------------------------------------------
# PROJCS["ED50 / UTM zone 30N", GEOGCS["ED50", DATUM["European Datum 1950", SPHEROID["International 1924", 6378388.0, 297.0, AUTHORITY["EPSG","7022"]], TOWGS84[-157.89, -17.16, -78.41, 2.118, 2.697, -1.434, -1.1097046576093785], AUTHORITY["EPSG","6230"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4230"]], PROJECTION["Transverse Mercator", AUTHORITY["EPSG","9807"]], PARAMETER["central_meridian", -3.0], PARAMETER["latitude_of_origin", 0.0], PARAMETER["scale_factor", 0.9996], PARAMETER["false_easting", 500000.0], PARAMETER["false_northing", 0.0], UNIT["m", 1.0], AXIS["Easting", EAST], AXIS["Northing", NORTH], AUTHORITY["EPSG","23030"]]
##------------------------------------------------------------
#.. CRS: Coordinate Reference System
kk <- make_EPSG()
kk <- kk[grep("ED50", kk$note), ]
kk <- kk[grep("^# ED50 / UTM zone 30N$", kk$note), ]
myCRS<- CRS("+init=epsg:23030") ; rm(kk)
proj4string(scA) <- myCRS
## Transforma DATUM
scA <- spTransform(scA, CRSobj = CRS("+proj=longlat +datum=WGS84"))

## Busca centroides de las secciones censales asocia datos a dichos puntos.
centro <- t( sapply(scA@polygons, function(e) {e@labpt} ) )
rownames(centro) <-sapply(scA@polygons, function(e) {e@ID} )
colnames(centro) <- c('long','latit')
# head(rownames (centro)) ; head(rownames(scA@data)) # los rownames son los @ID !!
centro <- SpatialPoints(centro, proj4string = CRS("+proj=longlat +datum=WGS84") )
centro<- SpatialPointsDataFrame(centro,scA@data, match.ID=TRUE)

### Representa para comprobación
area <- "Andalucia"
qmap(area, zoom = 7) + 
  geom_polygon(aes(x=long, y=lat, group=group), 
               data=scA,colour='black',fill='white',alpha=.4,size=.2) -> mymap
mymap  + geom_point(data = data.frame(coordinates(centro)),
                    aes(x=long, y=latit), 
                    colour='red',fill='blue',alpha=.5,size=.1)

save(scA,centro,file='mapSC01A.RData')