##################################################################
# /////////////////////////////////CARGAR ESTOS PROGRAMAS
Packages <- c("sp","sf","rgdal","maptools","raster","spatialEco","dplyr", "ggplot2")
lapply(Packages, library, character.only = TRUE)
library(extrafont)
library(xkcd)
library(ggthemes)
library(scales)
# //////////////////////////////////////////CARGAR DATOS
# ARCHIVOS EN FORMATO ESRI-SHP DESCARGADOS DEL SITIO DE CONABIO
# SISTEMA DE COORDENADAS GEOGRÁFICAS LAT-LONG WGS84
proy_GCS<- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
# /////////////////////////////////////////////////////////////Climas
setwd('/Datos GIS varios/clima1mgw') # USA LA RUTA A TU FICHERO
climas<- readShapePoly("clima1mgw.shp", proj4string=CRS("+proj=longlat"))
projection(climas) <- proy_GCS
# ////////////////////////////////////////////////Poblaciones (2010)
setwd('/Datos GIS varios/unikloc10gw') # USA LA RUTA A TU FICHERO
CitiesPts<- readShapePoints("unikloc10gw.shp", proj4string=CRS("+proj=longlat"))
projection(CitiesPts) <- proy_GCS
# ////////////////////////////////////////////////México país polígono
setwd('/Datos GIS varios/contdv250kgw') # USA LA RUTA A TU FICHERO
mexico<- readShapePoly("contdv250kgw.shp", proj4string=CRS("+proj=longlat"))
projection(mexico) <- proy_GCS
# ////////Modificamos la ubicación de estas localidades
shapefile <- add.xy(shapefile)
# Mexicali, -115.46742,32.65982 ; -114.76273,32.47288 San Luis Ro; Saltillo -101.00058,25.43285
which(CitiesPts@data$NOM_LOC == "Mexicali")
CitiesPts@coords[3699,1] = -115.46742 ; CitiesPts@coords[3699,2] = 32.65982
# San Luis Río Colorado
CitiesPts@data[151824,]; CitiesPts@coords[151824,1] = -114.76273 ; CitiesPts@coords[151824,2] = 32.47288
# Saltillo
CitiesPts@data[14855,];CitiesPts@coords[14855,1] = -101.00058 ; CitiesPts@coords[14855,2] = 25.43285
# ////////////////////////////////////// ANALIZAMOS CLIMAS /////////
tipos_clima= unique(climas@data$CLIMA_TIPO);sort(tipos_clima) #Factor
tipos_clima=as.character(tipos_clima); sort(tipos_clima) # As string
# Climas BW y BS
pattern <- "BW|BS"
arido_semiaridos=tipos_clima[grep(pattern, tipos_clima)]
length(arido_semiaridos) # 25 categorías
### Solo BWh o Árido cálido
pattern <- "^BWh|^BW\\(h\\'"
climas_arido_calido=tipos_clima[grep(pattern, tipos_clima)]; climas_arido_calido
length(climas_arido_calido) # 6 categorías
### Solo BWk o Árido templado
pattern <- "^BWk"
climas_arido_frio=tipos_clima[grep(pattern, tipos_clima)]; climas_arido_frio
length(climas_arido_frio) # 3 categorías
### Solo BSk o Semi-árido templado
pattern <-"^(?=.*BS)(?!.*h)"
climas_semiarido_k=tipos_clima[grep(pattern, tipos_clima, perl=TRUE)]; climas_semiarido_k
length(climas_semiarido_k) # 7
### Solo BSh o Semi-árido cálido
pattern <-"^(?=.*BS)(?!.*k)"
climas_semiarido_c=tipos_clima[grep(pattern, tipos_clima, perl=TRUE)]; climas_semiarido_c
length(climas_semiarido_c) # 9
# /////////////////////////////MAPA: SUBCONJUNTO POR TIPO DE CLIMAS/////////
arido_semiaridos # Tipos de clima árido y semiárido
# BWh, BWk, BSh y BSk, respectivamente.
climas_arido_calido; climas_arido_frio; climas_semiarido_c; climas_semiarido_k
# /////////////////////////////// MAPA 1 CLIMAS ARIDOS Y SEMIARIDOS /////////
climas_aridos=climas[which(climas@data$CLIMA_TIPO %in% arido_semiaridos),]
names(climas_aridos@data) ; climas_aridos@data$CLIMA_TIPO
# //////CREAMOS UNA VARIABLE CON LA CLASIFICACIÓN SIMPLE DE Köppen–Geiger
climas_aridos@data$tipoKG = ifelse(climas_aridos@data$CLIMA_TIPO %in% climas_arido_frio, "BWk",
ifelse(climas_aridos@data$CLIMA_TIPO %in% climas_semiarido_c, "BSh",
ifelse(climas_aridos@data$CLIMA_TIPO %in% climas_semiarido_k, "BSk",
"BWh")))
# VERIFICAMOS ALGUNOS
climas_aridos@data$tipoKG
climas_aridos@data$tipoKG[307]; climas_aridos@data$CLIMA_TIPO[307]
climas_aridos@data$tipoKG[260]; climas_aridos@data$CLIMA_TIPO[260]
climas_aridos@data$tipoKG[140]; climas_aridos@data$CLIMA_TIPO[140]
climas_aridos@data$tipoKG[10]; climas_aridos@data$CLIMA_TIPO[10]
# BWh #ed462f BWk #f29494 BSh #f5a331 BSk #fddb62
climas_aridos@data$colores = ifelse(climas_aridos@data$CLIMA_TIPO %in% climas_arido_frio, "#f29494",
ifelse(climas_aridos@data$CLIMA_TIPO %in% climas_semiarido_c, "#f5a331",
ifelse(climas_aridos@data$CLIMA_TIPO %in% climas_semiarido_k, "#fddb62",
"#ed462f")))
# /////////////////////////////////////////MAPA 1 /////////
plot(mexico, axes=T, col=NA, border="gray")
plot(climas_aridos, axes=T, col=climas_aridos@data$colores, border=NA, add = TRUE)
legend("topright", as.character(unique(climas_aridos@data$tipoKG)), fill = unique(climas_aridos@data$colores),
box.col="gray", inset=.02, title="Clima (Köppen–Geiger):")
# /////////////////////////////////////////TERMINA MAPA 1 /////////
# //////////////////////////////////ANALIZAMOS EL ARCHIVO DE CIUDADES /////////
names(CitiesPts@data)
head(CitiesPts@data)
unique(CitiesPts@data$Urb_Rur)
# [1] Urbana >= 15000; Rural (< 5000 habitantes)
# 5000 <= Mixta rural<10000; 10000 <= Mixta urbana < 15000
CitiesPts@data[1,]
urbanas <- subset(CitiesPts, Urb_Rur == "Urbana" | Urb_Rur =="Mixta urbana")
names(urbanas@data); unique(urbanas@data$Urb_Rur)
table(urbanas@data$Urb_Rur) # Mixta Urbana 300; Urbana 630
aggregate(urbanas@data$POBTOT, by=list(Category=urbanas@data$Urb_Rur), FUN=sum)
# Category x
# 1 Mixta urbana 3,664,946
# 2 Urbana 70,179,136
70179136+3664946
# 73844082 o casi 74 millones de personas, vivían en zonas urbanas el año 2010
urbanas@data$punto_t = ifelse(urbanas@data$Urb_Rur == "Urbana", 19, 18)
urbanas@data$colores = ifelse(urbanas@data$Urb_Rur == "Urbana", "#663300", "#cc00cc")
# /////////////////////////////////////////MAPA 2 /////////
plot(mexico, axes=T, col=NA, border="gray")
plot(urbanas, pch=urbanas@data$punto_t, cex= 0.5, col = urbanas@data$colores, add = TRUE)
legend("topright", as.character(unique(urbanas@data$Urb_Rur)),
col =unique(urbanas@data$colores), pch=c(19,18),
box.col="gray", inset=.02, title="Localidad")
# /////////////////////////////////////////TERMINA MAPA 2 /////////
# /////////////////////////////UN MAPA PARA CADA TIPO DE CLIMA /////////
# BWh
climas_arido_calido
climas_BWh=climas_aridos[which(climas_aridos@data$CLIMA_TIPO %in% climas_arido_calido),]
urbanas_BWh = urbanas[!is.na(over(urbanas,as(climas_BWh,"SpatialPolygons"))),]
names(urbanas_BWh@data)
urbanas_BWh@data$POBTOT[47]
sum(urbanas_BWh@data$POBTOT)
table(urbanas_BWh@data$Urb_Rur) # Mixta Urbana 12; Urbana 37
aggregate(urbanas_BWh@data$POBTOT, by=list(Category=urbanas_BWh@data$Urb_Rur), FUN=sum)
# Category x
# 1 Mixta urbana 149611
# 2 Urbana 4620228
# /////////////////////////////////////////MAPA 3 /////////
plot(mexico, axes=T, col=NA, border="gray")
plot(climas_BWh, axes=T, col=unique(climas_BWh@data$colores), border=NA, add = TRUE)
plot(urbanas_BWh, col=c("#663300", "#000000"), pch=unique(urbanas_BWh@data$punto_t), cex=0.5, add = TRUE)
title(main = "Zonas urbanas con clima \u{E1}rido c\u{E1}lido")
legend("topright",
legend = c(as.character(unique(urbanas@data$Urb_Rur)), "BWh"),
col =c(unique(urbanas@data$colores), unique(climas_BWh@data$colores)),
pch=c(19,18, NA),
fill =c(NA, NA,unique(climas_BWh@data$colores) ),
border = c(NA,NA,NA),
box.col="gray", inset=.02, title="Simbología")
# /////////////////////////////////////////TERMINA MAPA 3 /////////
# BWk
climas_arido_frio
climas_BWk=climas_aridos[which(climas_aridos@data$CLIMA_TIPO %in% climas_arido_frio),]
names(climas_BWk@data)
urbanas_BWk = urbanas[!is.na(over(urbanas,as(climas_BWk,"SpatialPolygons"))),]
dim(urbanas_BSk) # 95
which(urbanas@data$NOM_LOC == "Ensenada")
urbanas@data[8,]
names(urbanas_BWk@data)
urbanas_BWk@data$NOM_LOC
table(urbanas_BWk@data$Urb_Rur) # Mixta Urbana 0; Urbana 2
aggregate(urbanas_BWk@data$POBTOT, by=list(Category=urbanas_BWk@data$Urb_Rur), FUN=sum)
# Category x
# 1 Mixta urbana 0
# 2 Urbana 1337298
# /////////////////////////////////////////MAPA 4 /////////
plot(mexico, axes=T, col=NA, border="gray")
plot(climas_BWk, axes=T, col=unique(climas_BWk@data$colores), border=NA, add = TRUE)
plot(urbanas_BWk, col=c("#663300", "#000000"), pch=unique(urbanas_BWh@data$punto_t), cex=0.5, add = TRUE)
title(main = "Zonas urbanas con clima \u{E1}rido templado")
legend("topright",
legend = c(as.character(unique(urbanas@data$Urb_Rur)), "BWk"),
col =c(unique(urbanas@data$colores), unique(climas_BWk@data$colores)),
pch=c(19,18, NA),
fill =c(NA, NA,unique(climas_BWk@data$colores) ),
border = c(NA,NA,NA),
box.col="gray", inset=.02, title="Simbología")
# /////////////////////////////////////////TERMINA MAPA 4 /////////
# BSh
climas_semiarido_c
climas_BSh=climas_aridos[which(climas_aridos@data$CLIMA_TIPO %in% climas_semiarido_c),]
names(climas_BSh@data)
urbanas_BSh = urbanas[!is.na(over(urbanas,as(climas_BSh,"SpatialPolygons"))),]
table(urbanas_BSh@data$Urb_Rur) # Mixta Urbana 36; Urbana 76
aggregate(urbanas_BSh@data$POBTOT, by=list(Category=urbanas_BSh@data$Urb_Rur), FUN=sum)
# Category x
# 1 Mixta urbana 432833
# 2 Urbana 8803147
# /////////////////////////////////////////MAPA 5 /////////
plot(mexico, axes=T, col=NA, border="gray")
plot(climas_BSh, axes=T, col=unique(climas_BSh@data$colores), border=NA, add = TRUE)
plot(urbanas_BSh, col=c("#663300", "#000000"), pch=unique(urbanas_BSh@data$punto_t), cex=0.7, add = TRUE)
title(main = "Zonas urbanas con clima Semi\u{E1}rido c\u{E1}lido ")
legend("topright",
legend = c(as.character(unique(urbanas@data$Urb_Rur)), "BSh"),
col =c(unique(urbanas@data$colores), unique(climas_BSh@data$colores)),
pch=c(19,18, NA),
fill =c(NA, NA,unique(climas_BSh@data$colores) ),
border = c(NA,NA,NA),
box.col="gray", inset=.02, title="Simbología")
# /////////////////////////////////////////TERMINA MAPA 5 /////////
# BSk
climas_semiarido_k
climas_BSk=climas_aridos[which(climas_aridos@data$CLIMA_TIPO %in% climas_semiarido_k),]
names(climas_BSk@data)
urbanas_BSk = urbanas[!is.na(over(urbanas,as(climas_BSk,"SpatialPolygons"))),]
table(urbanas_BSk@data$Urb_Rur) # Mixta Urbana 27; Urbana 68
aggregate(urbanas_BSk@data$POBTOT, by=list(Category=urbanas_BSk@data$Urb_Rur), FUN=sum)
# Category x
# 1 Mixta urbana 330038
# 2 Urbana 11150847
# /////////////////////////////////////////MAPA 6 /////////
plot(mexico, axes=T, col=NA, border="gray")
plot(climas_BSk, axes=T, col=unique(climas_BSk@data$colores), border=NA, add = TRUE)
plot(urbanas_BSk, col=c("#663300", "#000000"), pch=unique(urbanas_BSk@data$punto_t), cex=0.7, add = TRUE)
title(main = "Zonas urbanas con clima Semi\u{E1}rido templado")
legend("topright",
legend = c(as.character(unique(urbanas@data$Urb_Rur)), "BSk"),
col =c(unique(urbanas@data$colores), unique(climas_BSk@data$colores)),
pch=c(19,18, NA),
fill =c(NA, NA,unique(climas_BSk@data$colores) ),
border = c(NA,NA,NA),
box.col="gray", inset=.02, title="Simbología")
# /////////////////////////////////////////TERMINA MAPA 6 /////////
### TABLA 1
# BWh Las diez ciudades con mayor población
BWh_df <- as(urbanas_BWh, "data.frame")
class(BWh_df); names(BWh_df); dim(BWh_df) # 49 Ciudades
unique(BWh_df$Urb_Rur)
sub1_cities_BWh= BWh_df %>%
filter (Urb_Rur == "Urbana") %>%
select(NOM_ENT, NOM_MUN, NOM_LOC,ALTITUD, POBTOT) %>%
arrange_(~ desc(POBTOT)) %>%
slice(1:10)
sub1_cities_BWh
setwd('/Users/COLSON/Desktop')
write.csv(sub1_cities_BWh, file = "sub1_cities_BWh.csv")
######################################
### TABLA 1
# BWk Las diez ciudades con mayor población
BWk_df <- as(urbanas_BWk, "data.frame")
class(BWk_df); names(BWk_df); dim(BWk_df) # 2
unique(BWk_df$Urb_Rur)
sub1_cities_BWk= BWk_df %>%
filter (Urb_Rur == "Urbana") %>%
select(NOM_ENT, NOM_MUN, NOM_LOC,ALTITUD, POBTOT) %>%
arrange_(~ desc(POBTOT)) %>%
slice(1:10)
sub1_cities_BWk
setwd('/Users/COLSON/Desktop')
write.csv(sub1_cities_BWk, file = "sub1_cities_BWk.csv")
######################################
### TABLA 2
# BSh Las diez ciudades con mayor población
BSh_df <- as(urbanas_BSh, "data.frame")
class(BSh_df); names(BSh_df); dim(BSh_df) # 112
unique(BSh_df$Urb_Rur)
sub1_cities_BSh= BSh_df %>%
filter (Urb_Rur == "Urbana" | Urb_Rur == "Mixta urbana") %>%
select(NOM_ENT, NOM_MUN, NOM_LOC,ALTITUD, POBTOT) %>%
arrange_(~ desc(POBTOT)) %>%
slice(1:10)
sub1_cities_BSh
setwd('/Users/COLSON/Desktop')
write.csv(sub1_cities_BSh, file = "sub1_cities_BSh.csv")
######################################
### TABLA 3
# BSk Las diez ciudades con mayor población
BSk_df <- as(urbanas_BSk, "data.frame")
class(BSk_df); names(BSk_df); dim(BSk_df) # 95
unique(BSk_df$Urb_Rur)
sub1_cities_BSk= BSk_df %>%
filter (Urb_Rur == "Urbana" | Urb_Rur == "Mixta urbana") %>%
select(NOM_ENT, NOM_MUN, NOM_LOC,ALTITUD, POBTOT) %>%
arrange_(~ desc(POBTOT)) %>%
slice(1:10)
sub1_cities_BSk
setwd('/Users/COLSON/Desktop')
write.csv(sub1_cities_BSk, file = "sub1_cities_BSk.csv")
######################################
BWh_df$clima="BWh"
BWk_df$clima="BWk"
BSh_df$clima="BSh"
BSk_df$clima="BSk"
pob_arids= rbind(BWh_df,BWk_df, BSh_df,BSk_df)
sum(pob_arids$POBTOT)
max(pob_arids$POBTOT); min(pob_arids$POBTOT)
dim(pob_arids)
df <- aggregate(formula=POBTOT~clima, data=pob_arids, FUN=sum)
df$colores=unique(climas_aridos@data$colores)
df$Pob_000=df$POBTOT/1000
df$Pob_000
### Eliminamos notación científica
options(scipen=999)
### XKCD theme
### Creamos el estilo de la gráfica
######## DESCARGAR FONTS PARA MAC
download.file("http://simonsoftware.se/other/xkcd.ttf", dest="xkcd.ttf", mode="wb")
system("mkdir ~/.fonts")
system("cp xkcd.ttf ~/.fonts")
font_import(paths = "~/.fonts", pattern = "[X/x]kcd", prompt=FALSE)
fonts()
fonttable()
if(.Platform$OS.type != "unix") {
## Register fonts for Windows bitmap output
loadfonts(device="win")
} else {
loadfonts()
}
### Humor Sans font descargado de: http://antiyawn.com/uploads/humorsans.html
### Se deben de instalar los Fonts en tu Mac, dejar archivo en el directorio de trabajo
### Explora en Internet la forma y procedimientos para instalar y usar los "fonts" de
### acuerdo a tu sistema operativo. Hoy en día es fácil conseguir gratis muchos tipos
### de "fonts"
font_import(".")
##################################
theme_xkcd <- theme(
panel.background = element_rect(fill="white"),
axis.ticks = element_line(colour=NA),
panel.grid = element_line(colour="white"),
axis.text.y = element_text(colour= "black"),
axis.text.x = element_text(colour="black"),
text = element_text(size=16, family="Humor Sans"),
plot.margin = margin(1, 1, 1, 1, "cm")
)
################################## Gráfica 1
ggplot(data=df, aes(reorder(clima, Pob_000),y=Pob_000)) +
geom_bar(stat="identity", color=NA, fill=df$colores[order(-df$Pob_000)]) + coord_flip() +
ggtitle("Habitantes en zonas urbanas por tipo de clima (2010)") +
xlab("Climas") + ylab("Personas x 1000") + scale_y_continuous(labels = comma) +
geom_text(aes(label = comma(Pob_000), y = Pob_000/2), size = 5, family="Humor Sans") +
theme_xkcd
################################## Gráfica 2
df1=aggregate(formula=POBTOT~clima+Urb_Rur, data=pob_arids, FUN=sum, na.rm=TRUE)
df1$Pob_000=df1$POBTOT/1000
df1
ggplot(data=df1, aes(reorder(clima, Pob_000),y=Pob_000, fill = factor(Urb_Rur))) +
geom_bar(stat="identity", color=NA, position=position_dodge()) + coord_flip() +
scale_fill_manual(values=c('#999999','#E69F00'), labels=c("Mixta-urbana", "Urbana"))+
ggtitle("Habitantes en zonas urbanas/mixtas por tipo de clima (2010)") +
guides(fill=guide_legend(title="tipo"))+
xlab("Climas") + ylab("Personas x 1000") + scale_y_continuous(labels = comma) +
geom_text(aes(label = comma(Pob_000), y = Pob_000 + 500), position = position_dodge(0.9),
vjust = -0.25, color = "black", family="Humor Sans", size = 5)+ theme_xkcd
#####################################################################
Ciudades áridas y semiáridas de méxico
Dr. Luis Alan Navarro Navarro
Centro de Estudios en Gobierno y Asuntos Públicos
El Colegio de Sonora
