Analisis descriptivo de datos espaciales con R: Segunda Parte
María Guzmán
7, 8 y 9 de noviembre, 2022
Librerías necesarias para el análisis
rm(list = ls())
library(readxl) # read_xlsx
library(ggplot2) # Gráficos
library(tidyverse) # %>%
library(knitr) # kable
library(dplyr) # mutate, rename, summarise, group_by
library(classInt) # classIntervals
library(spdep) # poly2nb
library(rgdal) # readOGR
library(raster) # crs
library(mxmaps) # df_mxmunicipio
library(sp) # coordinates
library(ggrepel) # geom_text_repel
library(spatialEco) # point.in.poly
library(raster) # shapefile
library(sf) # read_sf1 Introducción
Vamos a trabajar con el shapefile de los municipios de Guerrero y con variables climáticas del estado de Guerrero.
2 Lectura de un shapefile con read_sf()
ruta_shape<-"C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/guerrero_12/conjunto_de_datos"
guerrero.shp_sf<-read_sf(ruta_shape, layer="12mun")
class(guerrero.shp_sf) # "sf" "tbl_df" "tbl" "data.frame"-> [1] "sf" "tbl_df" "tbl" "data.frame"head(guerrero.shp_sf)-> Simple feature collection with 6 features and 4 fields
-> Geometry type: MULTIPOLYGON
-> Dimension: XY
-> Bounding box: xmin: 2620010 ymin: 524239.3 xmax: 2904337 ymax: 699944.3
-> Projected CRS: MEXICO_ITRF_2008_LCC
-> # A tibble: 6 × 5
-> CVEGEO CVE_ENT CVE_MUN NOMGEO geometry
-> <chr> <chr> <chr> <chr> <MULTIPOLYGON [m]>
-> 1 12001 12 001 "Acapulco de Ju\xe1rez" (((2723458 539108.7, 2723…
-> 2 12002 12 002 "Ahuacuotzingo" (((2818281 663676.7, 2818…
-> 3 12003 12 003 "Ajuchitl\xe1n del Progreso" (((2646281 699141.5, 2646…
-> 4 12004 12 004 "Alcozauca de Guerrero" (((2882116 623227.3, 2881…
-> 5 12005 12 005 "Alpoyeca" (((2872644 637916.8, 2872…
-> 6 12006 12 006 "Apaxtla" (((2717134 695405.2, 2717…3 Lectura de un shapefile con shapefile()
carpeta<-"C:\\Documetos MGM_OS_Dell\\Documentos Maria_Dell\\Cursos_Taller_Impartidos\\Analisis descriptivo de datos espaciales con R_43 Aniversario FM\\guerrero_12\\conjunto_de_datos\\"
guerrero.shp_shapefile<-shapefile(paste0(carpeta,"12mun.shp"))
class(guerrero.shp_shapefile)-> [1] "SpatialPolygonsDataFrame"
-> attr(,"package")
-> [1] "sp"head(guerrero.shp_shapefile)-> CVEGEO CVE_ENT CVE_MUN NOMGEO
-> 0 12001 12 001 Acapulco de Ju\xe1rez
-> 1 12002 12 002 Ahuacuotzingo
-> 2 12003 12 003 Ajuchitl\xe1n del Progreso
-> 3 12004 12 004 Alcozauca de Guerrero
-> 4 12005 12 005 Alpoyeca
-> 5 12006 12 006 Apaxtla4 Lectura de un shapefile con readOGR()
Los archivos shapefile, se van a descargar del Marco Geoestadístico. Censo de Población y Vivienda 2020, en el siguiente enlace:
4.1 Estado: Guerrero
Con la función readOGR() se puede leer el archivo
shapefile del estado de Guerrero.
Ruta del archivo
ruta_shape<-"C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/guerrero_12/conjunto_de_datos"Lectura de un archivo shapefile 12mun.shp:
guerrero.shp <-readOGR(ruta_shape, layer="12mun")-> OGR data source with driver: ESRI Shapefile
-> Source: "C:\Documetos MGM_OS_Dell\Documentos Maria_Dell\Cursos_Taller_Impartidos\Analisis descriptivo de datos espaciales con R_43 Aniversario FM\guerrero_12\conjunto_de_datos", layer: "12mun"
-> with 81 features
-> It has 4 fieldshead(guerrero.shp)-> CVEGEO CVE_ENT CVE_MUN NOMGEO
-> 0 12001 12 001 Acapulco de Ju\xe1rez
-> 1 12002 12 002 Ahuacuotzingo
-> 2 12003 12 003 Ajuchitl\xe1n del Progreso
-> 3 12004 12 004 Alcozauca de Guerrero
-> 4 12005 12 005 Alpoyeca
-> 5 12006 12 006 Apaxtlaguerrero.shp <-readOGR("C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/12_guerrero/conjunto_de_datos/12mun.shp")Gráfico del estado de Guerrero
plot(guerrero.shp,main="Estado de Guerrero")
El shapefile en R es un objeto del tipo SpatialPolygonsDataFrame:
class(guerrero.shp)-> [1] "SpatialPolygonsDataFrame"
-> attr(,"package")
-> [1] "sp"Características del objeto espacial:
crs(guerrero.shp)-> Coordinate Reference System:
-> Deprecated Proj.4 representation:
-> +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
-> +y_0=0 +ellps=GRS80 +units=m +no_defs
-> WKT2 2019 representation:
-> PROJCRS["MEXICO_ITRF_2008_LCC",
-> BASEGEOGCRS["MEXICO_ITRF_2008",
-> DATUM["International Terrestrial Reference Frame 2008",
-> ELLIPSOID["GRS 1980",6378137,298.257222101,
-> LENGTHUNIT["metre",1]],
-> ID["EPSG",1061]],
-> PRIMEM["Greenwich",0,
-> ANGLEUNIT["Degree",0.0174532925199433]]],
-> CONVERSION["unnamed",
-> METHOD["Lambert Conic Conformal (2SP)",
-> ID["EPSG",9802]],
-> PARAMETER["Latitude of false origin",12,
-> ANGLEUNIT["Degree",0.0174532925199433],
-> ID["EPSG",8821]],
-> PARAMETER["Longitude of false origin",-102,
-> ANGLEUNIT["Degree",0.0174532925199433],
-> ID["EPSG",8822]],
-> PARAMETER["Latitude of 1st standard parallel",17.5,
-> ANGLEUNIT["Degree",0.0174532925199433],
-> ID["EPSG",8823]],
-> PARAMETER["Latitude of 2nd standard parallel",29.5,
-> ANGLEUNIT["Degree",0.0174532925199433],
-> ID["EPSG",8824]],
-> PARAMETER["Easting at false origin",2500000,
-> LENGTHUNIT["metre",1],
-> ID["EPSG",8826]],
-> PARAMETER["Northing at false origin",0,
-> LENGTHUNIT["metre",1],
-> ID["EPSG",8827]]],
-> CS[Cartesian,2],
-> AXIS["(E)",east,
-> ORDER[1],
-> LENGTHUNIT["metre",1,
-> ID["EPSG",9001]]],
-> AXIS["(N)",north,
-> ORDER[2],
-> LENGTHUNIT["metre",1,
-> ID["EPSG",9001]]]]El objeto guerrero.shp tiene una proyección:
+proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs
4.2 Cambio de coordenadas
4.2.1 WGS84
Cambio del Sistema de Coordenadas de Referencia (CRS) del objeto guerrero.shp, se le va asignar el CRS WGS84.
guerrero.shp_wgs84<- sp::spTransform(guerrero.shp,
sp::CRS("+proj=longlat +datum=WGS84 +no_defs"))
crs(guerrero.shp_wgs84) # sistema de referencia de coordenadas (proyección)-> Coordinate Reference System:
-> Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
-> WKT2 2019 representation:
-> GEOGCRS["unknown",
-> DATUM["World Geodetic System 1984",
-> ELLIPSOID["WGS 84",6378137,298.257223563,
-> LENGTHUNIT["metre",1]],
-> ID["EPSG",6326]],
-> PRIMEM["Greenwich",0,
-> ANGLEUNIT["degree",0.0174532925199433],
-> ID["EPSG",8901]],
-> CS[ellipsoidal,2],
-> AXIS["longitude",east,
-> ORDER[1],
-> ANGLEUNIT["degree",0.0174532925199433,
-> ID["EPSG",9122]]],
-> AXIS["latitude",north,
-> ORDER[2],
-> ANGLEUNIT["degree",0.0174532925199433,
-> ID["EPSG",9122]]]]guerrero.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white", color="grey30")+ # Shapefile
labs(x="Longitud",y="Latitud", size=0.5)+ # Nombres ejes
theme_minimal()
4.2.2 UTM
SRC México: ITRF2008/UTM zone
México ITRF2008 / UTM zone 11N (EPSG:6366)
México ITRF2008 / UTM zone 12N (EPSG:6367)
México ITRF2008 / UTM zone 13N (EPSG:6368)
México ITRF2008 / UTM zone 14N (EPSG:6369)
México ITRF2008 / UTM zone 15N (EPSG:6370)
México ITRF2008 / UTM zone 16N (EPSG:6371)
Guerrero está en la zona 14
guerreo.shp_utm<- sp::spTransform(guerrero.shp,
sp::CRS("+proj=utm +zone=14+south=T ellps=WGS84"))
crs(guerreo.shp_utm) # sistema de referencia de coordenadas (proyección)-> Coordinate Reference System:
-> Deprecated Proj.4 representation:
-> +proj=utm +zone=14 +ellps=WGS84 +units=m +no_defs
-> WKT2 2019 representation:
-> PROJCRS["unknown",
-> BASEGEOGCRS["unknown",
-> DATUM["Unknown based on WGS84 ellipsoid",
-> ELLIPSOID["WGS 84",6378137,298.257223563,
-> LENGTHUNIT["metre",1],
-> ID["EPSG",7030]]],
-> PRIMEM["Greenwich",0,
-> ANGLEUNIT["degree",0.0174532925199433],
-> ID["EPSG",8901]]],
-> CONVERSION["UTM zone 14N",
-> METHOD["Transverse Mercator",
-> ID["EPSG",9807]],
-> PARAMETER["Latitude of natural origin",0,
-> ANGLEUNIT["degree",0.0174532925199433],
-> ID["EPSG",8801]],
-> PARAMETER["Longitude of natural origin",-99,
-> ANGLEUNIT["degree",0.0174532925199433],
-> ID["EPSG",8802]],
-> PARAMETER["Scale factor at natural origin",0.9996,
-> SCALEUNIT["unity",1],
-> ID["EPSG",8805]],
-> PARAMETER["False easting",500000,
-> LENGTHUNIT["metre",1],
-> ID["EPSG",8806]],
-> PARAMETER["False northing",0,
-> LENGTHUNIT["metre",1],
-> ID["EPSG",8807]],
-> ID["EPSG",16014]],
-> CS[Cartesian,2],
-> AXIS["(E)",east,
-> ORDER[1],
-> LENGTHUNIT["metre",1,
-> ID["EPSG",9001]]],
-> AXIS["(N)",north,
-> ORDER[2],
-> LENGTHUNIT["metre",1,
-> ID["EPSG",9001]]]]guerreo.shp_utm%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white", color="grey30")+ # Shapefile
labs(x="X",y="Y", size=0.5)+ # Nombres ejes
theme_minimal()
4.3 Selección de un conjunto de polígonos
Guerrero se divide en 7 regiones.
Creación de los vectores de las 7 regiones de Guerrero, usando la clave de municipio.
Acapulco <-c("001") # CVE_MUN
CostaChica<-c("012","013","018","025","023","030","036","080","077","046",
"052","053","056","062","071")
CostaGrande<-c("011","014","016","021","068","048","057","038")
Centro<-c("002","028","029","075","032","079","039","040","042","044",
"051","061","074")
LaMontana<-c("076","004","005","009","010","078","020","024","033","081","041","043","045","063","065","066","069","070","072")
RegionNorte<-c("006","008","015","017","019","026","034","035","031","037",
"047","049","055","058","059","060")
TierraCaliente<-c("003","007","022","027","050","054","064","067","073")Selección de los municipios de Costa Chica:
guerrero.shp_wgs84_CostaChica<-guerrero.shp_wgs84[guerrero.shp_wgs84$CVE_MUN%in%CostaChica,]
class(guerrero.shp_wgs84_CostaChica)-> [1] "SpatialPolygonsDataFrame"
-> attr(,"package")
-> [1] "sp"guerrero.shp_wgs84_CostaChica%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white", color="grey30")+ # Shapefile
labs(x="Longitud",y="Latitud", size=0.5)+ # Nombres ejes
theme_minimal()
4.4 Etiquetado de poligonos
Información del Censo del 2015:
head(df_mxmunicipio)-> state_code municipio_code region state_name state_name_official
-> 1 01 001 01001 Aguascalientes Aguascalientes
-> 2 01 002 01002 Aguascalientes Aguascalientes
-> 3 01 003 01003 Aguascalientes Aguascalientes
-> 4 01 004 01004 Aguascalientes Aguascalientes
-> 5 01 005 01005 Aguascalientes Aguascalientes
-> 6 01 006 01006 Aguascalientes Aguascalientes
-> state_abbr state_abbr_official municipio_name year pop pop_male
-> 1 AGS Ags. Aguascalientes 2015 877190 425731
-> 2 AGS Ags. Asientos 2015 46464 22745
-> 3 AGS Ags. Calvillo 2015 56048 27298
-> 4 AGS Ags. Cosío 2015 15577 7552
-> 5 AGS Ags. Jesús María 2015 120405 60135
-> 6 AGS Ags. Pabellón de Arteaga 2015 46473 22490
-> pop_female afromexican part_afromexican indigenous part_indigenous
-> 1 451459 532 2791 104125 14209
-> 2 23719 3 130 1691 92
-> 3 28750 10 167 7358 2223
-> 4 8025 0 67 2213 191
-> 5 60270 32 219 8679 649
-> 6 23983 3 74 6232 251
-> metro_area long lat
-> 1 Aguascalientes -102.2960 21.87982
-> 2 <NA> -102.0893 22.23832
-> 3 <NA> -102.7188 21.84691
-> 4 <NA> -102.3000 22.36641
-> 5 Aguascalientes -102.3434 21.96127
-> 6 <NA> -102.2765 22.14920Creación de un data.frame con los centroides de los municipios de Costa Chica
# Obtención de los centroides de los municipios de Costa Chica:
centroides_constachica<-data.frame(coordinates(guerrero.shp_wgs84_CostaChica), # Centroides de Costa Chica
CVE_MUN=guerrero.shp_wgs84_CostaChica$CVE_MUN)%>%
rename(long=X1, lat=X2) # Re-nombrando las coordenadas de los centroides
head(centroides_constachica)-> long lat CVE_MUN
-> 11 -99.06754 16.95727 012
-> 12 -98.57240 16.68360 013
-> 17 -98.93060 16.60803 018
-> 22 -98.51767 16.46197 023
-> 23 -98.75067 16.87683 052
-> 24 -99.40904 16.80869 053#--- data.frame con la información de los municipios de Costa Chica
mun_CostaChica<-df_mxmunicipio%>%
filter(state_name=="Guerrero", municipio_code%in%CostaChica) %>%
dplyr::select(municipio_code, municipio_name) %>%
rename(CVE_MUN=municipio_code,# Se renombre la variable municipio_code con CVE_MUN
lab_mun=municipio_name)
head(mun_CostaChica)-> CVE_MUN lab_mun
-> 1 012 Ayutla de los Libres
-> 2 013 Azoyú
-> 3 018 Copala
-> 4 023 Cuajinicuilapa
-> 5 025 Cuautepec
-> 6 030 Florencio Villarreal#--- Creación del data.frame con los nombre de los municipios de Costa Chica
lab_CostaChica<-merge(centroides_constachica, mun_CostaChica, by="CVE_MUN")
head(lab_CostaChica)-> CVE_MUN long lat lab_mun
-> 1 012 -99.06754 16.95727 Ayutla de los Libres
-> 2 013 -98.57240 16.68360 Azoyú
-> 3 018 -98.93060 16.60803 Copala
-> 4 023 -98.51767 16.46197 Cuajinicuilapa
-> 5 025 -98.96371 16.71854 Cuautepec
-> 6 030 -99.13422 16.69293 Florencio VillarrealPara el siguiente gráfico se ocupa:
El shapefile de la región Costa Chica: guerrero.shp_wgs84_CostaChica
El data.frame de las etiquetas: lab_CostaChica
guerrero.shp_wgs84_CostaChica%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white", color="grey30")+ # Shapefile
geom_text_repel(data=lab_CostaChica, aes(long, lat, label=lab_mun),size=2)+ # data.frame lab
labs(x="Longitud",y="Latitud", size=1)+ # Nombres ejes
theme_minimal()
La función geom_text_repel() ayuda que las etiquetas de
los municipios no se encimen.
5 Estaciones de monitoreo
5.1 Creación de un data.frame
Estaciones meteorológicas de Guerrero
estaciones<-read_excel("C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/Estaciones_guerrero.xlsx")
head(estaciones)-> # A tibble: 6 × 14
-> ESTACION NOMBRE MUNIC…¹ LAT_G LAT_M LAT_S LON_G LON_M LON_S ALTITUD
-> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
-> 1 12001 ACAPETLAHUAYA TELOLO… 18 20 7 100 4 21 1292
-> 2 12003 AGUA SALADA (CFE) ACAPUL… 17 10 43 99 37 52 233
-> 3 12004 AHUEHUEPAN IGUALA… 18 20 17 99 38 48 760
-> 4 12005 ALCOZAUCA (SMN) ALCOZA… 17 27 0 98 23 0 1462
-> 5 12006 APANGO (CFE) MARTIR… 17 44 21 99 19 46 1065
-> 6 12007 ARATICHANGUIO ZIRAND… 18 28 33 101 21 33 226
-> # … with 4 more variables: Fecha_INICIO <dttm>, Fecha_FIN <chr>,
-> # PORCENTAJE <dbl>, DATOS <dbl>, and abbreviated variable name ¹MUNICIPIO5.2 Creación de la variable fecha
En excel, la variable de fecha debe de estar capturada con dia/mes/año, una forma de hacerlo es la siguiente:
estaciones$F_Fin<-paste(estaciones$Fecha_FIN,"/01", sep="") # Para agregar el díaDefiniendo las variables de tiempo con la función
as.Date()
est_2000_2017<-estaciones %>%
mutate(F_Inicio=as.Date(Fecha_INICIO), F_Fin=as.Date(F_Fin))%>%
dplyr:: filter(F_Fin>= "2000-01-01")
head(est_2000_2017)-> # A tibble: 6 × 16
-> ESTACION NOMBRE MUNIC…¹ LAT_G LAT_M LAT_S LON_G LON_M LON_S ALTITUD
-> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
-> 1 12001 ACAPETLAHUAYA TELOLO… 18 20 7 100 4 21 1292
-> 2 12003 AGUA SALADA (CFE) ACAPUL… 17 10 43 99 37 52 233
-> 3 12004 AHUEHUEPAN IGUALA… 18 20 17 99 38 48 760
-> 4 12006 APANGO (CFE) MARTIR… 17 44 21 99 19 46 1065
-> 5 12007 ARATICHANGUIO ZIRAND… 18 28 33 101 21 33 226
-> 6 12008 ARCELIA ARCELIA 18 19 1 100 16 48 414
-> # … with 6 more variables: Fecha_INICIO <dttm>, Fecha_FIN <chr>,
-> # PORCENTAJE <dbl>, DATOS <dbl>, F_Fin <date>, F_Inicio <date>, and
-> # abbreviated variable name ¹MUNICIPIOSon 183 estaciones que siguen funcionando después del año 2000:
dim(est_2000_2017) -> [1] 183 165.3 Conversión de coordenadas
Las coordenadas GPS están en formato de grados (d), minutos (m) y segundos (s); para convertirlo a DD ( decimal degree) se debe puede usar la siguiente fórmula:
\[ DD=D+\frac{M}{60}+\frac{S}{3600} \]
attach(est_2000_2017) # Adjunta la base de datos
est_2000_2017$lat <- (LAT_G+LAT_M/60+ LAT_S/3600) # Creación de la variable lat
est_2000_2017$long<- -1*(LON_G+LON_M/60+ LON_S/3600) # Creación de la variable long
head(est_2000_2017)-> # A tibble: 6 × 18
-> ESTACION NOMBRE MUNIC…¹ LAT_G LAT_M LAT_S LON_G LON_M LON_S ALTITUD
-> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
-> 1 12001 ACAPETLAHUAYA TELOLO… 18 20 7 100 4 21 1292
-> 2 12003 AGUA SALADA (CFE) ACAPUL… 17 10 43 99 37 52 233
-> 3 12004 AHUEHUEPAN IGUALA… 18 20 17 99 38 48 760
-> 4 12006 APANGO (CFE) MARTIR… 17 44 21 99 19 46 1065
-> 5 12007 ARATICHANGUIO ZIRAND… 18 28 33 101 21 33 226
-> 6 12008 ARCELIA ARCELIA 18 19 1 100 16 48 414
-> # … with 8 more variables: Fecha_INICIO <dttm>, Fecha_FIN <chr>,
-> # PORCENTAJE <dbl>, DATOS <dbl>, F_Fin <date>, F_Inicio <date>, lat <dbl>,
-> # long <dbl>, and abbreviated variable name ¹MUNICIPIOPara el siguiente gráfico se ocupa:
shapefile: guerrero.shp_wgs84
data.frame de las estaciones: est_2000_2017
guerrero.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white",color="grey30")+
geom_point(data=est_2000_2017, aes(x=long,y=lat, fill=ESTACION), # Estaciones
shape = 21,colour="black",fill = "grey", size = 2)+
labs(x="Longitud",y="Latitud",title = "183 estaciones seguian funcionando después del año 2000.")+
theme_minimal()
Etiquetado de las estaciones:
CVE_Est<-substr(est_2000_2017$ESTACION,start=3, stop=5) # Extracción de las claves de los municipios
est_2000_2017$CVE_Est<-as.numeric(CVE_Est)
head(est_2000_2017)-> # A tibble: 6 × 19
-> ESTACION NOMBRE MUNIC…¹ LAT_G LAT_M LAT_S LON_G LON_M LON_S ALTITUD
-> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
-> 1 12001 ACAPETLAHUAYA TELOLO… 18 20 7 100 4 21 1292
-> 2 12003 AGUA SALADA (CFE) ACAPUL… 17 10 43 99 37 52 233
-> 3 12004 AHUEHUEPAN IGUALA… 18 20 17 99 38 48 760
-> 4 12006 APANGO (CFE) MARTIR… 17 44 21 99 19 46 1065
-> 5 12007 ARATICHANGUIO ZIRAND… 18 28 33 101 21 33 226
-> 6 12008 ARCELIA ARCELIA 18 19 1 100 16 48 414
-> # … with 9 more variables: Fecha_INICIO <dttm>, Fecha_FIN <chr>,
-> # PORCENTAJE <dbl>, DATOS <dbl>, F_Fin <date>, F_Inicio <date>, lat <dbl>,
-> # long <dbl>, CVE_Est <dbl>, and abbreviated variable name ¹MUNICIPIOPara el siguiente gráfico se ocupa:
shapefile: guerrero.shp_wgs84
data.frame: dat_2000_2017
guerrero.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group),
fill="white",color="grey30")+
geom_point(data=est_2000_2017, # datos del data.frame
aes(x=long,y=lat, fill=CVE_Est),
shape = 21,colour="black",fill = "pink", size = 2)+
geom_text(data=est_2000_2017, aes(long,lat, label =CVE_Est), size=2)+ # data.frame etiq
labs(x="Longitud",y="Latitud",title = "Estaciones meteorológicas en Guerrero")+
theme_minimal()
5.3.1 Intersección de puntos y polígonos
Selección de las variables:
est_2000_2017_SPointsDF<-est_2000_2017 %>%
dplyr::select(long,lat,ESTACION)
head(est_2000_2017_SPointsDF)-> # A tibble: 6 × 3
-> long lat ESTACION
-> <dbl> <dbl> <dbl>
-> 1 -100. 18.3 12001
-> 2 -99.6 17.2 12003
-> 3 -99.6 18.3 12004
-> 4 -99.3 17.7 12006
-> 5 -101. 18.5 12007
-> 6 -100. 18.3 12008Creación del objeto SpatialPointsDataFrame:
coordinates(est_2000_2017_SPointsDF) = ~long+lat # Creación de SpatialPointsDataFrame
class(est_2000_2017_SPointsDF)-> [1] "SpatialPointsDataFrame"
-> attr(,"package")
-> [1] "sp"Asignación de el sistema de coordenadas WGS84 al objeto SpatialPointsDataFrame, est_2000_2017_SPointsDF:
proj4string(est_2000_2017_SPointsDF)=CRS("+proj=longlat +datum=WGS84 +no_defs") # WGS84La función point.in.poly() permite intersectar un
shapefile con un SpatialPointsDataFrame
pts.poly <- point.in.poly(est_2000_2017_SPointsDF, guerrero.shp_wgs84)
class(pts.poly) #SpatialPointsDataFrame-> [1] "SpatialPointsDataFrame"
-> attr(,"package")
-> [1] "sp"View(pts.poly@data)Con el código anterior, ya se puede identificar a que municipio pertenece cada una de las estaciones meteorológicas.
pts.poly.df<-as.data.frame(pts.poly)
head(pts.poly.df)-> ESTACION CVEGEO CVE_ENT CVE_MUN NOMGEO coords.x1
-> 1 12001 12058 12 058 Teloloapan -100.07250
-> 2 12003 12001 12 001 Acapulco de Ju\xe1rez -99.63111
-> 3 12004 12035 12 035 Iguala de la Independencia -99.64667
-> 4 12006 12042 12 042 M\xe1rtir de Cuilapan -99.32944
-> 5 12007 12073 12 073 Zir\xe1ndaro -101.35917
-> 6 12008 12007 12 007 Arcelia -100.28000
-> coords.x2
-> 1 18.33528
-> 2 17.17861
-> 3 18.33806
-> 4 17.73917
-> 5 18.47583
-> 6 18.31694dim(pts.poly.df)-> [1] 183 7Identificación de las estaciones que están en cada uno de los municipios:
Est_Mun<-pts.poly.df %>%
dplyr::select(ESTACION,CVE_MUN)
head(Est_Mun)-> ESTACION CVE_MUN
-> 1 12001 058
-> 2 12003 001
-> 3 12004 035
-> 4 12006 042
-> 5 12007 073
-> 6 12008 0075.4 Calculo de las estadísticas desciptivas
datos<-read_excel("C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/GRO_DLY_Variables climaticas.xlsx")Creación de la variable Fecha:
datos$Fecha<-paste(datos$YEAR_MONTH,"/01", sep="") # Se agrega el día a la fechaSelección de la variable Temperatura máxima \(°C\), VALUE_1 , y definición de la variable
tiempo con la función as.Date(),
dat_sel<-datos%>%
mutate(Fecha=as.Date(Fecha),ESTACION=Station_ID)%>%
dplyr::select("ESTACION","ELEMENT_CODE","Fecha","VALUE_1") %>%
filter(VALUE_1>0,ELEMENT_CODE==2,Fecha>= "2000-01-01")
dim(dat_sel) -> [1] 23595 4head(dat_sel)-> # A tibble: 6 × 4
-> ESTACION ELEMENT_CODE Fecha VALUE_1
-> <dbl> <dbl> <date> <dbl>
-> 1 12001 2 2000-01-01 23
-> 2 12001 2 2000-02-01 22.5
-> 3 12001 2 2000-03-01 30
-> 4 12001 2 2000-05-01 31
-> 5 12001 2 2000-06-01 31
-> 6 12001 2 2000-08-01 29Calculo del promedio de las precipitaciones:
dat_mean<-dat_sel %>%
group_by(ESTACION) %>%
summarise(Promedio.Tem.Max=mean(VALUE_1))
head(dat_mean)-> # A tibble: 6 × 2
-> ESTACION Promedio.Tem.Max
-> <dbl> <dbl>
-> 1 12001 32.9
-> 2 12003 35.8
-> 3 12004 33.7
-> 4 12006 30.7
-> 5 12007 38.8
-> 6 12008 36.5Asignación de las coordenadas a las estaciones
Est_coord<-est_2000_2017 %>%
dplyr::select(long,lat,ESTACION)
dat_mean_coord<-merge(dat_mean, Est_coord, by="ESTACION")
dim(dat_mean_coord) # Número de Estaciones-> [1] 181 4head(dat_mean_coord)-> ESTACION Promedio.Tem.Max long lat
-> 1 12001 32.85054 -100.07250 18.33528
-> 2 12003 35.83333 -99.63111 17.17861
-> 3 12004 33.71494 -99.64667 18.33806
-> 4 12006 30.66667 -99.32944 17.73917
-> 5 12007 38.81117 -101.35917 18.47583
-> 6 12008 36.47409 -100.28000 18.31694Se identifican las estaciones con sus municipios:
dat_mean_coord.mun<-merge(dat_mean_coord, Est_Mun, by="ESTACION")
dim(dat_mean_coord.mun) # Número de estaciones-> [1] 181 5head(dat_mean_coord.mun)-> ESTACION Promedio.Tem.Max long lat CVE_MUN
-> 1 12001 32.85054 -100.07250 18.33528 058
-> 2 12003 35.83333 -99.63111 17.17861 001
-> 3 12004 33.71494 -99.64667 18.33806 035
-> 4 12006 30.66667 -99.32944 17.73917 042
-> 5 12007 38.81117 -101.35917 18.47583 073
-> 6 12008 36.47409 -100.28000 18.31694 0075.4.1 Creación: SpatialPointsDataFrame
Creación de un objeto espacial de puntos; SpatialPointsDataFrame
dat_mean_coord.mun_SPointsDF<-dat_mean_coord.mun
coordinates(dat_mean_coord.mun_SPointsDF) = ~long+lat # Creación de SpatialPointsDataFrame
class(dat_mean_coord.mun_SPointsDF)-> [1] "SpatialPointsDataFrame"
-> attr(,"package")
-> [1] "sp"6 Generación de gráficos
6.1 Polígono vs valores
Para la generación del siguiente gráfico se ocupa:
shapefile: guerrero.shp_wgs84
SpatialPointsDataFrame: dat_mean_coord_SPointsDF
plot(guerrero.shp_wgs84, axes=T, xlab="X Coord", ylab="Y Coord", main=" ",lty=1)
plot(dat_mean_coord.mun_SPointsDF, add=T, pch=21,col="blue") # Ubicación de las estaciones
Para el siguiente gráfico se ocupa:
shapefile: guerrero.shp_wgs84
data.frame: dat_mean_coord.mun
guerrero.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group),
fill="white",color="grey30")+
geom_point(data=dat_mean_coord.mun, # datos del data.frame
aes(x=long,y=lat, fill=Promedio.Tem.Max),
shape = 21,colour="black",fill = "blue", size = 2)+
labs(x="Longitud",y="Latitud")+
theme_minimal()
6.2 Puntos de colores
Para el siguiente gráfico se ocupa:
Un shapefile: guerrero.shp_wgs84
Un data.frame: dat_mean_coord.mun
guerrero.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white",color="grey30")+ # Shapefile
geom_point(data=dat_mean_coord.mun, aes(x=long,y=lat, color=Promedio.Tem.Max), size = 2)+ # data.frame
coord_fixed(ratio = 1) + # Colores
scale_color_gradient(low = "blue", high = "orange") +
#geom_text_repel(data=lab_gerrero, aes(long, lat, label=lab),size=2.5)+ # Etiquetas
labs(x="Longitud",y="Latitud",color="Prom. Tem. Max.", title = "")+
theme_minimal()
En el mapa podemos observar que los niveles de precipitación más altos
sueles ocurrir en los municipios que están en las costas.
6.3 Polígonos de colores
Conversión de un objeto SpatialPolygonsDataFrame a un objeto data.frame:
guerrero.shp_wgs84_d.f<-ggplot2::fortify(guerrero.shp_wgs84, region="CVE_MUN") # Conversión de SpatialPolygonsDataFrame a data.frame
guerrero.shp_wgs84_d.f$CVE_MUN<-guerrero.shp_wgs84_d.f$id # Creación de la variable CVE_MUNPara la generación del gráfico se ocupa:
El data.frame: dat_mean_coord.mun
El data.frame de un shapefile: guerrero.shp_wgs84_d.f
ggplot(dat_mean_coord.mun, aes(map_id = CVE_MUN)) + # map_id es la columna con la clave de unión. Nos ahorra el join previo.
geom_map(aes(fill=Promedio.Tem.Max), map = guerrero.shp_wgs84_d.f) + # objeto SPDF convertido a data.frame.
expand_limits(x =guerrero.shp_wgs84_d.f$long, y = guerrero.shp_wgs84_d.f$lat)+# Se indica el tamaño del mapa al máximo de lat y long.
labs(x="Longitud",y="Latitud", fill="Lluvia promedio")+
theme_bw()
Los municipios en blanco, fue porque no hubo estaciones con registros.
Gráfico coloreado por rangos:
Val.intervalo: variable creada con las función
classIntervals()data.frame: dat_mean_coord.mun
Shapefile convertido a data.frame: guerrero.shp_wgs84_d.f
#--- Creación de los intervalos
ngrupos<- 3
intervalo <- classIntervals(dat_mean_coord.mun$Promedio.Tem.Max, n =ngrupos, style = "quantile")
punt.cor<-intervalo$brks
offs <- 0.0000001 # para que no aparezca NA
punt.cor[1] <- punt.cor[1] - offs
punt.cor[length(punt.cor)] <- punt.cor[length(punt.cor)] + offs
# --- Creación de la variable
dat_mean_coord.mun <- dat_mean_coord.mun%>%
mutate(Val.intervalo = cut(Promedio.Tem.Max, punt.cor,dig.lab=8)) # con dig.lab se indican el número se cifras
#--- Generación del gráfico
ggplot(dat_mean_coord.mun, # datos,
aes(map_id = CVE_MUN)) +# map_id es la columna con la clave de unión. Nos ahorra el join previo.
geom_map(aes(fill = Val.intervalo ),color="black",map = guerrero.shp_wgs84_d.f) + # map= objeto SPDF convertido a data.frame.
scale_fill_brewer(palette = "red") +
expand_limits(x =guerrero.shp_wgs84_d.f$long, y = guerrero.shp_wgs84_d.f$lat)+ #Especificamos el tamaño del mapa al máximo de lat y long.
labs(x="Longitud",y="Latitud", fill="Prom. Tem. max", title = "")+
theme_minimal()+
theme(legend.position = c(.3, .4),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6))+
guides(colour = guide_legend(nrow = 1, override.aes = list(size = 4)))
6.4 Grafico con etiquetas
Creación del data.frame de las etiquetas:
# Obtención de los centroides de los municipios de guerreto:
centroides_guerrero<-data.frame(coordinates(guerrero.shp_wgs84), # Centroides de los Estados
CVE_MUN=guerrero.shp_wgs84$CVE_MUN)%>%
rename(long=X1, lat=X2) # Re-nombrando las coordenadas de los centroides
head(centroides_guerrero)-> long lat CVE_MUN
-> 0 -99.71424 16.95725 001
-> 1 -98.95137 17.75227 002
-> 2 -100.60773 17.91902 003
-> 3 -98.36611 17.37590 004
-> 4 -98.51057 17.64224 005
-> 5 -99.98090 18.09123 006#--- Información de los municipios de Guerrero
mun_guerrero<-df_mxmunicipio%>%
filter(state_code=="12")%>%
rename(CVE_MUN=municipio_code,# Se renombre la variable municipio_code con CVE_MUN
NOMGEO=municipio_name)
head(mun_guerrero)-> state_code CVE_MUN region state_name state_name_official state_abbr
-> 1 12 001 12001 Guerrero Guerrero GRO
-> 2 12 002 12002 Guerrero Guerrero GRO
-> 3 12 003 12003 Guerrero Guerrero GRO
-> 4 12 004 12004 Guerrero Guerrero GRO
-> 5 12 005 12005 Guerrero Guerrero GRO
-> 6 12 006 12006 Guerrero Guerrero GRO
-> state_abbr_official NOMGEO year pop pop_male pop_female
-> 1 Gro. Acapulco de Juárez 2015 810669 385812 424857
-> 2 Gro. Ahuacuotzingo 2015 26858 12543 14315
-> 3 Gro. Ajuchitlán del Progreso 2015 38134 18592 19542
-> 4 Gro. Alcozauca de Guerrero 2015 19368 9102 10266
-> 5 Gro. Alpoyeca 2015 6657 3163 3494
-> 6 Gro. Apaxtla 2015 11159 5358 5801
-> afromexican part_afromexican indigenous part_indigenous metro_area long
-> 1 77837 8320 164595 13867 Acapulco -99.88656
-> 2 1128 52 22378 497 <NA> -98.93491
-> 3 20 48 903 162 <NA> -100.49606
-> 4 103 17 18792 15 <NA> -98.38389
-> 5 232 20 4275 22 <NA> -98.51324
-> 6 133 192 1506 275 <NA> -99.93231
-> lat
-> 1 16.86187
-> 2 17.71409
-> 3 18.15956
-> 4 17.46352
-> 5 17.67006
-> 6 18.13008#--- Creación del data.frame con los nombre de los municipios de Guerrero
lab_guerrero<-merge(centroides_guerrero, mun_guerrero%>%dplyr::select(CVE_MUN, NOMGEO), by="CVE_MUN")
head(lab_guerrero)-> CVE_MUN long lat NOMGEO
-> 1 001 -99.71424 16.95725 Acapulco de Juárez
-> 2 002 -98.95137 17.75227 Ahuacuotzingo
-> 3 003 -100.60773 17.91902 Ajuchitlán del Progreso
-> 4 004 -98.36611 17.37590 Alcozauca de Guerrero
-> 5 005 -98.51057 17.64224 Alpoyeca
-> 6 006 -99.98090 18.09123 ApaxtlaPara el siguiente gráfico se ocupa:
shapefile: guerrero.shp_wgs84
data.frame de etiquetas: lab_guerrero
guerrero.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white", color="grey30")+ # Shapefile
geom_text(data=lab_guerrero, aes(long, lat, label =NOMGEO), size=1.5)+ # data.frame etiq
labs(x="Longitud",y="Latitud", size=0.5)+ # Nombres ejes
theme_minimal()
#--- Creación de los intervalos
ngrupos<- 3
intervalo <- classIntervals(dat_mean_coord.mun$Promedio.Tem.Max, n =ngrupos, style = "quantile")
punt.cor<-intervalo$brks
offs <- 0.0000001 # para que no aparezca NA
punt.cor[1] <- punt.cor[1] - offs
punt.cor[length(punt.cor)] <- punt.cor[length(punt.cor)] + offs
# --- Creación de la variable
dat_mean_coord.mun <- dat_mean_coord.mun%>%
mutate(Val.intervalo = cut(Promedio.Tem.Max, punt.cor,dig.lab=8)) # con dig.lab se indican el número se cifras
#--- Generación del gráfico
ggplot(dat_mean_coord.mun, # datos,
aes(map_id = CVE_MUN)) +# map_id es la columna con la clave de unión. Nos ahorra el join previo.
geom_map(aes(fill = Val.intervalo ),color="black",map = guerrero.shp_wgs84_d.f) + # map= objeto SPDF convertido a data.frame.
scale_fill_brewer(palette = "red") +
expand_limits(x =guerrero.shp_wgs84_d.f$long, y = guerrero.shp_wgs84_d.f$lat)+ # Especificamos el tamaño del mapa al máximo de lat y long.
geom_text(data=lab_guerrero, aes(long, lat, label =NOMGEO), size=1.5)+ # data.frame label
labs(x="Longitud",y="Latitud", fill="Pro. Tem. Max", title = "Temperatura promedio máxima del 2000 al 2017")+
theme_minimal()+
theme(legend.position = c(1, 1),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6))+
guides(colour = guide_legend(nrow = 1, override.aes = list(size = 4)))6.5 Grafico de regiones
mun_guerrero<-df_mxmunicipio%>%
filter(state_code=="12")
head(mun_guerrero)-> state_code municipio_code region state_name state_name_official state_abbr
-> 1 12 001 12001 Guerrero Guerrero GRO
-> 2 12 002 12002 Guerrero Guerrero GRO
-> 3 12 003 12003 Guerrero Guerrero GRO
-> 4 12 004 12004 Guerrero Guerrero GRO
-> 5 12 005 12005 Guerrero Guerrero GRO
-> 6 12 006 12006 Guerrero Guerrero GRO
-> state_abbr_official municipio_name year pop pop_male pop_female
-> 1 Gro. Acapulco de Juárez 2015 810669 385812 424857
-> 2 Gro. Ahuacuotzingo 2015 26858 12543 14315
-> 3 Gro. Ajuchitlán del Progreso 2015 38134 18592 19542
-> 4 Gro. Alcozauca de Guerrero 2015 19368 9102 10266
-> 5 Gro. Alpoyeca 2015 6657 3163 3494
-> 6 Gro. Apaxtla 2015 11159 5358 5801
-> afromexican part_afromexican indigenous part_indigenous metro_area long
-> 1 77837 8320 164595 13867 Acapulco -99.88656
-> 2 1128 52 22378 497 <NA> -98.93491
-> 3 20 48 903 162 <NA> -100.49606
-> 4 103 17 18792 15 <NA> -98.38389
-> 5 232 20 4275 22 <NA> -98.51324
-> 6 133 192 1506 275 <NA> -99.93231
-> lat
-> 1 16.86187
-> 2 17.71409
-> 3 18.15956
-> 4 17.46352
-> 5 17.67006
-> 6 18.13008Creación de una variable Region7 que contiene las 7 regiones de Guerrero:
# Creación de la variable Region7 en el data.frame lab_guerrero
lab_guerrero$Region7<-lab_guerrero$CVE_MUN%>% # data.frame con la variable a re-categorizar CVE_MUN
fct_collapse("Acapúlco"=Acapulco, "Costa Chica"=CostaChica,
"Costa Grande"=CostaGrande, Centro=Centro,
"La Montaña"=LaMontana, "Región Norte"=RegionNorte,
"Tierra Caliente"=TierraCaliente, "Tierra Caliente"=TierraCaliente)
head(lab_guerrero)-> CVE_MUN long lat NOMGEO Region7
-> 1 001 -99.71424 16.95725 Acapulco de Juárez Acapúlco
-> 2 002 -98.95137 17.75227 Ahuacuotzingo Centro
-> 3 003 -100.60773 17.91902 Ajuchitlán del Progreso Tierra Caliente
-> 4 004 -98.36611 17.37590 Alcozauca de Guerrero La Montaña
-> 5 005 -98.51057 17.64224 Alpoyeca La Montaña
-> 6 006 -99.98090 18.09123 Apaxtla Región NortePara el gráfico de las regiones, se ocupa:
data.frame: lab_guerrero
Un shapefle convertido a data.frame: gerrero.shp_wgs84_d.f
ggplot(lab_guerrero,
aes(map_id=CVE_MUN))+
geom_map(aes(fill=Region7), # Variable a graficar
map=guerrero.shp_wgs84_d.f, color="grey30")+
expand_limits(x=guerrero.shp_wgs84_d.f$long,
y=guerrero.shp_wgs84_d.f$lat)+
labs(x="Longitud",y="Latitud",fill="Regiones",title = "Regiones de México")+
theme_minimal()+
theme(
legend.position = c(.2, .7),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)
)
7 Apéndice
7.1 Uso de teclas
La tecla shift es una fecha hacia arriba en el teclado.
Combinaciones de teclas en RStudio:
Ejecución de código R: ctrl + enter
Ayuda: F1 (fn+F1)
Abrir el visualizador,
View(): ctrl + click (con mause, es el izquierdo)Insertar el símbolo de asignación
<-, presionar: alt + F1Inserta el operador de asignación
<-, presionar: alt + -Inserta un pipe
%>%, presionar: ctrl+ shift + mCambiar entre las ventanas en RStudio: ctrl+ 1:8
Navegar entre pestañas (.Rmd, .R, tablas, etc): ctrl + tabulador (<- ->)
Renombrar un objeto creado, en barra de herramientas: Code-> Rename in Scope
Muestra el panel de combinaciones de teclas: alt +shift +k
Se arrastran líneas: alt+up (up=re pág)
Pedir ayuda al R: fn+F1
Cursor multilínea: ctrl+alt +arriba; se termina con la tecla esc
Modificar selección, presionar: ctrl+alt+shift+m
7.2 Paquetes de R
El paquete rgdal, tiene funciones para leer y escribir objetos vectoriales;
SpatialPointsDataFrame,
SpatialLinesDataFrame
SpacialPolygonsDataFrame
7.3 Funciones de R
- La función coordinates() transforma un data.frame en un objeto espacial.
7.4 Descarga de R y RStudio
La página oficial del software estadístico R es : https://www.r-project.org/
Antes de realizar la instalación del R, ir a Sistema del equipo de Computo y verificar el número de bits que se tiene 64 0 32, ya que este dato se pregunta durante la instalación de R.
La instalación de R, RTools y RStudio, se realiza en el siguiente orden:
Instalar R versión 4.2.1 para Windows: https://cran.r-project.org/bin/windows/base/
Instalar RTools versión 4.2: https://cran.r-project.org/bin/windows/Rtools/
Instalar RStudio: https://www.rstudio.com/products/rstudio/download/
Nota:
Rtools es un paquete de cadena de herramientas que se utiliza para compilar paquetes R desde el origen (aquellos que necesitan compilación de código C/C++ o Fortran) y para compilar R en sí.
Cuando se utiliza la ubicación de instalación predeterminada, R y Rtools42 se pueden instalar en cualquier orden y Rtools42 se puede instalar cuando R ya se está ejecutando
8 Información
INEGI
Censo de Población y Vivienda 2020, base de datos: https://www.inegi.org.mx/programas/ccpv/2020/#Datos_abiertos
Principales resultados por AGEB y manzana urbana 2020: https://www.inegi.org.mx/app/scitel/Default?ev=10
Mortalidad: https://www.inegi.org.mx/programas/mortalidad/?ps=Microdatos#Tabulados
IMSS:
- Bases dedatos: http://datos.imss.gob.mx/dataset/informacion-en-salud/resource/5c865fa6-f632-48e4-8858-54d30350f107
Enlaces para análisis espaciales:
- Ecología, estadística y ciencia de datos con R: https://cmerow.github.io/RDataScience/04_Spatial.html
+Introducción a GIS y análisis espacial: https://mgimond.github.io/Spatial/coordinate-systems-in-r.html
An Introduction to R for Spatial Analysis and Mapping: https://study.sagepub.com/brunsdon2e
Geospatial Data Science in R: https://zia207.github.io/geospatial-r-github.io/spatial-autocorrelation.html