Skip to content

Instantly share code, notes, and snippets.

@viciana
Created February 20, 2017 15:48
Show Gist options
  • Save viciana/c0d092bc3b4d97466cfb31ed99d1deeb to your computer and use it in GitHub Desktop.
Save viciana/c0d092bc3b4d97466cfb31ed99d1deeb to your computer and use it in GitHub Desktop.
Descarga, lectura y representación de las secciones censales del censo de 2001. Andalucía
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')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment