Análisis descriptivo de datos espaciales con R: Parte I
María Guzmán
7, 8 y 9 de noviembre, 2022
Librerías necesarias para el análisis
rm(list = ls()) # Eliminar objetos creados
library(rgdal) # readOGR Para importar shapefiles.
library(raster) # crs
library(ggplot2) # Gráficos
library(tidyverse) # %>%
library(knitr) # kable
library(dplyr) # rename, group_by, summarise,
library(sf) # st_crs
library(classInt) # classIntervals
library(sp) # coordinates
library(readxl) # read_xlsx
library(writexl) # write_xlsx
library(forcats) # fct_collapse
library(maptools) # unionSpatialPolygons
library(spatialEco) # point.in.poly
library(geoR) # locations.inside
library(mxmaps) # df_mxmunicipio
library(rgeos) # gUnaryUnion1 Introducción
1.1 Resumen
El objetivo del Curso-Taller es utilizar el R como un sistema de información geográfica (SIG), también citado como GIS Geographical Information System; para calcular las principales estadísticas descriptivas de centralidad y dispersión de un conjunto de datos dado, y mostrar los resultados en un mapa, archivos shapefile. Entre los objetos de R que aprenderemos a manejar están los SpatialPointsDataFrame y SpatialPolygonsDataFrame.
1.2 SIG
Un Sistema de Información Geográfica (SIG) es una serie de procedimientos que integra tecnología informática, personas e información geográfica, cuya principal función es procesar, analizar y representar datos georreferenciados, es decir, que están asociados a un lugar específico (Olaya, 2014).
1.3 Software SIG
Software para los SIG disponibles:
ArcGIS: El entorno de escritorio de ArcGIS abarca un conjunto de aplicaciones que incluyen ArcMap, ArcCatalog, ArcScene y ArcGlobe. ArcGIS viene en tres niveles de licencia diferentes (básico, estándar y avanzado).
QGIS: Es un software GIS de código abierto, abarca la mayor parte de la funcionalidad incluida en ArcGIS.
R
1.4 Sistemas de referencia de coordenadas
Sistemas de referencia de coordenadas (SRC):
WGS84: Sistema geodésico mundial de 1984. (EPSG:4326)
ETRS89: Sistema de referencia terrestre europeo de 1989 muy similar al WGS84.
NAD83: Datum estadounidense de 1983 el cual es muy similar al WGS84. (EPSG:4269)
PSAD56: Datum provisional sudamericano de 1956. (4248)
SIRGAS: Sistema de Referencia Geocéntrico para las Américas.
UTM: Universal Transverse Mercator.
1.5 Shapefile
Un shapefile es un formato de representación vectorial desarrollado por ESRI (Enviromental Systems Research Institute). Consta de un número variable de archivos, en los que se almacena digitalmente la localización de los elementos geográficos (archivo shape .shp) junto con sus atributos o características ( .dbf)
La lista de archivos que definen un shapefile, se muestra a continuación:
.shp: La geometría
.dbf: Información de atributos o base de datos
.shx: almacena el índice de las geometrías
.aih (.ain): Índice de atributo
.prj: Información del sistema de coordenadas,
.sbn (.sbx): Archivo de índice espacial
.cpg: Opcional, usado para especificar la página de código
2 Lectura de un shapefile
Los archivos shapefile, se van a descargar del Marco Geoestadístico. Censo de Población y Vivienda 2020, en el sigiente enlace:
Los shapefile de INEGI, traen 5 archivos:
.dbf
.prj
.shp
.shx
.cpg
2.1 País: México
Cargando un shapefile con la función con la función
readOGR()
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/889463807469_s_shp_nacional/MG_2020_Integrado/conjunto_de_datos"Lectura de un archivo shapefile (00ent):
mexico.shp <-readOGR(ruta_shape, layer="00ent")-> 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\889463807469_s_shp_nacional\MG_2020_Integrado\conjunto_de_datos", layer: "00ent"
-> with 32 features
-> It has 3 fieldshead(mexico.shp@data)-> CVEGEO CVE_ENT NOMGEO
-> 0 01 01 Aguascalientes
-> 1 02 02 Baja California
-> 2 03 03 Baja California Sur
-> 3 04 04 Campeche
-> 4 05 05 Coahuila de Zaragoza
-> 5 06 06 Colimaplot(mexico.shp, main="Mapa de México")
mexico.shp <-readOGR("C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/889463807469_s_shp_nacional/MG_2020_Integrado/conjunto_de_datos/00ent.shp")# Lectura de Shapefile
mexico.shp <-readOGR("C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Investigacion_MGM/Bases de Datos e Informacion/Shapefile_archivos/Municipios_2016/conjunto_de_datos/integracion_territorial.shp")El shapefile en R es un objeto del tipo SpatialPolygonsDataFrame:
class(mexico.shp)-> [1] "SpatialPolygonsDataFrame"
-> attr(,"package")
-> [1] "sp"Características del objeto espacial:
crs(mexico.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 mexico.shp tiene una proyección:
- proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5
Atributos del objeto espacial, shapefile:
head(mexico.shp)-> CVEGEO CVE_ENT NOMGEO
-> 0 01 01 Aguascalientes
-> 1 02 02 Baja California
-> 2 03 03 Baja California Sur
-> 3 04 04 Campeche
-> 4 05 05 Coahuila de Zaragoza
-> 5 06 06 ColimaOtra forma de ver los atributos del objeto shapefile:
head(mexico.shp@data)-> CVEGEO CVE_ENT NOMGEO
-> 0 01 01 Aguascalientes
-> 1 02 02 Baja California
-> 2 03 03 Baja California Sur
-> 3 04 04 Campeche
-> 4 05 05 Coahuila de Zaragoza
-> 5 06 06 Colima2.2 Cambio de coordenadas
Con la función spTransform() se cambia el Sistema de
Coordenadas de Referencia del objeto mexico.shp:
mexico.shp_wgs84<- sp::spTransform(mexico.shp,
sp::CRS("+proj=longlat +datum=WGS84 +no_defs"))Las funciones st_crs() y crs() muestra el
sistema de coordenadas de referencias del objeto shapefile:`
st_crs(mexico.shp_wgs84) # sistema de referencia de coordenadas (proyección)-> Coordinate Reference System:
-> User input: +proj=longlat +datum=WGS84 +no_defs
-> wkt:
-> 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]]]]crs(mexico.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]]]]2.3 Selección de un conjunto de polígonos
México se divide en ocho regiones.
Noroeste:
| Estado | CVE_ENT |
|---|---|
| Baja California | 02 |
| Baja California Sur | 03 |
| Chihuahua | 08 |
| Durango | 10 |
| Sinaloa | 25 |
| Sonora | 26 |
Noreste:
| Estado | CVE_ENT |
|---|---|
| Coahuila | 05 |
| Nuevo León | 19 |
| Tamaulipas | 28 |
Oeste:
| Estado | CVE_ENT |
|---|---|
| Colima | 06 |
| Jalisco | 14 |
| Michoacán | 16 |
| Nayarit | 18 |
Oriente:
| Estado | CVE_ENT |
|---|---|
| Hidalgo | 13 |
| Puebla | 21 |
| Tlaxcala | 29 |
| Veracruz | 30 |
Centro-Norte:
| Estado | CVE_ENT |
|---|---|
| Aguascalientes | 01 |
| Guanajuato | 11 |
| Querétaro | 22 |
| San Luis Potosí | 24 |
| Zacatecas | 32 |
Centro-Sur:
| Estado | CVE_ENT |
|---|---|
| Ciudad de México | 09 |
| México | 15 |
| Morelos | 17 |
Suroeste:
| Estado | CVE_ENT |
|---|---|
| Chiapas | 07 |
| Guerrero | 12 |
| Oaxaca | 20 |
Sureste:
| Estado | CVE_ENT |
|---|---|
| Campeche | 04 |
| Quintana Roo | 23 |
| Tabasco | 27 |
| Yucatán | 31 |
Creación de los vectores de las regiones
Nor.Oeste_cla <-c("02", "03", "08", "10", "25", "26")
Nor.Este_cla <-c("05", "19", "28")
Oeste_cla <-c("06", "14", "16", "18")
Oriente_cla <-c("13", "21", "29", "30")
Centro.Norte_cla <-c("01", "11", "22", "24", "32")
Centro.Sur_cla <-c("09","15","17")
Sur.Oeste_cla <-c("07","12","20")
Sur.Este_cla <-c("04","23","27","31")Selección de la región Suroeste de México:
mexico.shp_wgs84_suroeste<-mexico.shp_wgs84[mexico.shp_wgs84$CVE_ENT%in%Sur.Oeste_cla,]
plot(mexico.shp_wgs84_suroeste, main="Región Suroeste de México")
2.4 Etiquetado de poligonos
Para la región del Suroeste, se seleccionan las etiquetas
centroides_suroeste<-coordinates(mexico.shp_wgs84_suroeste) # Centroides de las 32 entidades de México
# Creación del data.frame con las etiquetes
head(mexico.shp_wgs84) # En la tercer columna están los nombres de los Estados-> CVEGEO CVE_ENT NOMGEO
-> 0 01 01 Aguascalientes
-> 1 02 02 Baja California
-> 2 03 03 Baja California Sur
-> 3 04 04 Campeche
-> 4 05 05 Coahuila de Zaragoza
-> 5 06 06 Colimalab_suroeste<-data.frame(long=centroides_suroeste[,1],
lat=centroides_suroeste[,2],
lab_NOMGEO=mexico.shp_wgs84_suroeste$NOMGEO
)
head(lab_suroeste)-> long lat lab_NOMGEO
-> 6 -92.47291 16.48524 Chiapas
-> 11 -99.92183 17.66801 Guerrero
-> 19 -96.43013 16.96145 OaxacaPara el siguiente mapa se ocupa:
El shapefile de la región Suroeste: mexico.shp_wgs84_suroeste
El data.frame de las etiquetas: lab_suroeste
mexico.shp_wgs84_suroeste%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white", color="grey30")+ # Shapefile
geom_text(data=lab_suroeste, aes(long, lat, label =lab_NOMGEO), size=2.5)+ # data.frame lab
labs(x="Longitud",y="Latitud", size=0.5)+ # Nombres ejes
theme_minimal()
3 Base de datos
3.1 Estadísticas desciptivas
Dado \(n\) datos el promedio está dado por: \[ \bar{x}=\frac{1}{n}\sum_{i=1}^{n} x_{i} \] La varianza muestral está dada por
\[ S^2=\frac{1}{n-1}\sum_{i=1}^{n} (x_{i}-\bar{x})^2 \]
El máximos de un conjunto de datos es el valor más grande; mientras que el mínimo de un conjunto de datos es el valor más pequeño.
3.2 Creación de un data.frame
# Datos -----
# Población del Censo de 1990
dat_1990<-df_mxstate_1990_2010 %>%
filter(year=="1990" & sex=="Total")%>%
dplyr::select("region", "state_name_official", "pop") %>%
rename(pop.1990=pop,state_name=state_name_official,CVE_ENT=region)
# Población del Censo de 1995
dat_1995<-df_mxstate_1990_2010 %>%
filter(year=="1995" & sex=="Total")%>%
dplyr::select("region", "state_name_official", "pop")
# Población del Censo del 2000
dat_2000<-df_mxstate_1990_2010 %>%
filter(year=="2000" & sex=="Total")%>%
dplyr::select("region", "state_name_official", "pop")
# Población del Censo del 2005
dat_2005<-df_mxstate_1990_2010 %>%
filter(year=="2005" & sex=="Total")%>%
dplyr::select("region", "state_name_official", "pop")
# Población del Censo del 2010
dat_2010<-df_mxstate_1990_2010 %>%
filter(year=="2010" & sex=="Total")%>%
dplyr::select("region", "state_name_official", "pop")
# Población del Censo del 2015
dat_2015<-df_mxstate_2015%>%
dplyr::select("region", "state_name", "pop")
# Población del Censo del 2020
dat_2020<-df_mxstate_2020%>%
dplyr::select("region", "state_name", "pop")Guardemos las bases de datos
write_xlsx(df_mxstate_1990_2010, "C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/df_mxstate_1990_2010.xlsx")
write_xlsx(df_mxstate_2015, "C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/df_mxstate_2015.xlsx")
write_xlsx(df_mxstate_2020, "C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Cursos_Taller_Impartidos/Analisis descriptivo de datos espaciales con R_43 Aniversario FM/df_mxstate_2020.xlsx")Población de los censos de INEGI:
# Final----
dat_1990_2020<-data.frame(dat_1990,
pop.1995=dat_1995$pop,
pop.2000=dat_2000$pop,
pop.2005=dat_2005$pop,
pop.2010=dat_2010$pop,
pop.2015=dat_2015$pop,
pop.2020=dat_2020$pop
)
View(dat_1990_2020)
head(dat_1990_2020)-> CVE_ENT state_name pop.1990 pop.1995 pop.2000 pop.2005 pop.2010
-> 1 01 Aguascalientes 719659 862720 944285 1065416 1184996
-> 2 02 Baja California 1660855 2112140 2487367 2844469 3155070
-> 3 03 Baja California Sur 317764 375494 424041 512170 637026
-> 4 04 Campeche 535185 642516 690689 754730 822441
-> 5 05 Coahuila de Zaragoza 1972340 2173775 2298070 2495200 2748391
-> 6 06 Colima 428510 488028 542627 567996 650555
-> pop.2015 pop.2020
-> 1 1312544 1425607
-> 2 3315766 3769020
-> 3 712029 798447
-> 4 899931 928363
-> 5 2954915 3146771
-> 6 711235 731391dim(dat_1990_2020) # Dimensiones del objeto creado-> [1] 32 9#dat_1990_2020$Region8<-Region8 # Asignación de las regiones al data.frameNombres de los estados de México
centroides_mexico<-coordinates(mexico.shp_wgs84)
lab_mexico<-data.frame(long=centroides_mexico[,1],
lat=centroides_mexico[,2],
CVE_ENT=mexico.shp_wgs84@data$CVE_ENT,
lab_NOMGEO=mexico.shp_wgs84@data$NOMGEO)
head(lab_mexico)-> long lat CVE_ENT lab_NOMGEO
-> 0 -102.36194 22.00644 01 Aguascalientes
-> 1 -115.11002 30.58375 02 Baja California
-> 2 -112.07892 25.92532 03 Baja California Sur
-> 3 -90.35982 18.84017 04 Campeche
-> 4 -102.04404 27.29544 05 Coahuila de Zaragoza
-> 5 -103.91337 19.14099 06 ColimaLa función apply() ayuda a calcular los promedios por
estado de los 7 Censos de población del INEGI:
dat_mean<-data.frame(CVE_ENT=dat_1990_2020$CVE_ENT,
Val.mean=apply(dat_1990_2020[,3:9], 1, mean))
dat_mean_lab<-merge(dat_mean,lab_mexico,by="CVE_ENT")
head(dat_mean_lab)-> CVE_ENT Val.mean long lat lab_NOMGEO
-> 1 01 1073603.9 -102.36194 22.00644 Aguascalientes
-> 2 02 2763526.7 -115.11002 30.58375 Baja California
-> 3 03 539567.3 -112.07892 25.92532 Baja California Sur
-> 4 04 753407.9 -90.35982 18.84017 Campeche
-> 5 05 2541351.7 -102.04404 27.29544 Coahuila de Zaragoza
-> 6 06 588620.3 -103.91337 19.14099 ColimaCálculo de la desviación estándar
dat_sd<-data.frame(CVE_ENT=dat_1990_2020$CVE_ENT,
Val.Sd=apply(dat_1990_2020[,3:9], 1, sd))
dat_sd_lab<-merge(dat_sd,lab_mexico,by="CVE_ENT")
head(dat_sd_lab)-> CVE_ENT Val.Sd long lat lab_NOMGEO
-> 1 01 251643.3 -102.36194 22.00644 Aguascalientes
-> 2 02 729417.2 -115.11002 30.58375 Baja California
-> 3 03 180957.5 -112.07892 25.92532 Baja California Sur
-> 4 04 141831.8 -90.35982 18.84017 Campeche
-> 5 05 428298.0 -102.04404 27.29544 Coahuila de Zaragoza
-> 6 06 113687.2 -103.91337 19.14099 ColimaDeterminación del valor máximo de los 8 años de registro
dat_max<-data.frame(CVE_ENT=dat_1990_2020$CVE_ENT,
Val.max=apply(dat_1990_2020[,3:9], 1, max))
dat_max_lab<-merge(dat_max,lab_mexico,by="CVE_ENT")
head(dat_max_lab)-> CVE_ENT Val.max long lat lab_NOMGEO
-> 1 01 1425607 -102.36194 22.00644 Aguascalientes
-> 2 02 3769020 -115.11002 30.58375 Baja California
-> 3 03 798447 -112.07892 25.92532 Baja California Sur
-> 4 04 928363 -90.35982 18.84017 Campeche
-> 5 05 3146771 -102.04404 27.29544 Coahuila de Zaragoza
-> 6 06 731391 -103.91337 19.14099 Colima3.3 Creación: SpatialPointsDataFrame
Creación de un objeto espacial de puntos; SpatialPointsDataFrame
dat_mean_SPointsDF<-dat_mean_lab
coordinates(dat_mean_SPointsDF) = ~long+lat # Creación de SpatialPointsDataFrame Características de objeto de puntos espacial:
class(dat_mean_SPointsDF)-> [1] "SpatialPointsDataFrame"
-> attr(,"package")
-> [1] "sp"crs(dat_mean_SPointsDF)-> Coordinate Reference System:
-> Deprecated Proj.4 representation: NA3.4 Creación: SpatialPolygonsDataFrame
Para la creación de un objeto SpatialPolygonsDataFrame, se utiliza la función merge, y dos objetos:
SpatialPolygonsDataFrame: mexico.shp_wgs84
data.frame: dat_mean.df, el data.frame puede ser sustituido por un SpatialPointsDataFrame (dat_mean_SPointsDF).
mexico.shp_wgs84_dat<-merge(mexico.shp_wgs84, dat_mean_lab, by="CVE_ENT")
class(mexico.shp_wgs84_dat) # SpatialPolygonsDataFrame-> [1] "SpatialPolygonsDataFrame"
-> attr(,"package")
-> [1] "sp"head(mexico.shp_wgs84_dat@data,4) # Muestra los primeros 6 elementos-> CVE_ENT CVEGEO NOMGEO Val.mean long lat
-> 1 01 01 Aguascalientes 1073603.9 -102.36194 22.00644
-> 2 02 02 Baja California 2763526.7 -115.11002 30.58375
-> 3 03 03 Baja California Sur 539567.3 -112.07892 25.92532
-> 4 04 04 Campeche 753407.9 -90.35982 18.84017
-> lab_NOMGEO
-> 1 Aguascalientes
-> 2 Baja California
-> 3 Baja California Sur
-> 4 Campeche3.5 Guardando un shapefile
writeOGR(mexico.shp_wgs84_dat, "C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Investigacion_MGM/TRABAJOS de Tesis/Jose Ernesto", layer="mexico.shp_wgs84_dat", driver="ESRI Shapefile")3.6 Identificación de puntos dentro un polígono
Creación de un SpatialPointsDataFrame:
head(loc_inside)-> long lat
-> 1 -98.88529 17.33003
-> 2 -94.01534 17.33003
-> 3 -91.58036 17.33003
-> 4 -106.19023 18.72900
-> 5 -94.01534 18.72900
-> 6 -91.58036 18.72900head(y)-> rnorm.dim.loc_inside..1...0..1.
-> 1 -1.4200188
-> 2 1.8781689
-> 3 1.0594234
-> 4 0.8164765
-> 5 -1.4649089
-> 6 0.1701553dat_spdf<-SpatialPointsDataFrame(coords=loc_inside, # SpatialPointsDataFrame
data=y,
proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs"))
class(dat_spdf)-> [1] "SpatialPointsDataFrame"
-> attr(,"package")
-> [1] "sp"plot(dat_spdf)
pts.poly <- point.in.poly(dat_spdf, mexico.shp_wgs84 ) # intersecta polígono y datos: SpatialPointsDataFrame
kable(head(pts.poly,10))| rnorm.dim.loc_inside..1…0..1. | CVEGEO | CVE_ENT | NOMGEO |
|---|---|---|---|
| -1.4200188 | 12 | 12 | Guerrero |
| 1.8781689 | 30 | 30 | Veracruz de Ignacio de la Llave |
| 1.0594234 | 07 | 07 | Chiapas |
| 0.8164765 | NA | NA | NA |
| -1.4649089 | NA | NA | NA |
| 0.1701553 | 04 | 04 | Campeche |
| -1.6547056 | 23 | 23 | Quintana Roo |
| -0.1923944 | 31 | 31 | Yucatán |
| -0.5320741 | 24 | 24 | San Luis Potosí |
| 1.2803444 | NA | NA | NA |
pts.poly.df<-as.data.frame(pts.poly)
kable(head(pts.poly.df, 10))| rnorm.dim.loc_inside..1…0..1. | CVEGEO | CVE_ENT | NOMGEO | coords.x1 | coords.x2 |
|---|---|---|---|---|---|
| -1.4200188 | 12 | 12 | Guerrero | -98.88529 | 17.33003 |
| 1.8781689 | 30 | 30 | Veracruz de Ignacio de la Llave | -94.01534 | 17.33003 |
| 1.0594234 | 07 | 07 | Chiapas | -91.58036 | 17.33003 |
| 0.8164765 | NA | NA | NA | -106.19023 | 18.72900 |
| -1.4649089 | NA | NA | NA | -94.01534 | 18.72900 |
| 0.1701553 | 04 | 04 | Campeche | -91.58036 | 18.72900 |
| -1.6547056 | 23 | 23 | Quintana Roo | -89.14538 | 18.72900 |
| -0.1923944 | 31 | 31 | Yucatán | -89.14538 | 20.12796 |
| -0.5320741 | 24 | 24 | San Luis Potosí | -98.88529 | 21.52693 |
| 1.2803444 | NA | NA | NA | -96.45032 | 21.52693 |
pts.poly.df_ent<-pts.poly.df%>%
group_by(CVE_ENT)%>%
summarise(Total=sum(rnorm.dim.loc_inside..1...0..1.))
kable(head(pts.poly.df_ent, 10))| CVE_ENT | Total |
|---|---|
| 02 | 0.4447366 |
| 03 | -0.5418224 |
| 04 | 0.1701553 |
| 07 | 1.0594234 |
| 08 | -3.2288913 |
| 10 | 0.9093124 |
| 12 | -1.4200188 |
| 19 | -0.8050720 |
| 23 | -1.6547056 |
| 24 | -0.5320741 |
#residuales_SPDF<-merge(mexico.shp_wgs84, pts.poly.df_ent, by="CVE_MUN")4 Gráficos
4.1 Polígono vs valores
Para la generación del siguiente gráfico se ocupa:
shapefile: mexico.shp_wgs84
SpatialPointsDataFrame: dat_mean_SPointsDF
plot(mexico.shp_wgs84, axes=T, xlab="X Coord", ylab="Y Coord", main=" ",lty=1)
plot(dat_mean_SPointsDF, add=T, pch=21,col="blue")
Para el siguiente gráfico se ocupa:
shapefile: mexico.shp_wgs84
data.frame: dat_mean_lab
mexico.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group),
fill="white",color="grey30")+
geom_point(data=dat_mean_lab, # datos
aes(x=long,y=lat, fill=Val.mean),
shape = 21,colour="black",fill = "blue", size = 2)+
labs(x="Longitud",y="Latitud")+
theme_dark() 
4.2 Puntos de colores
Para generar el siguiente gráfico se ocupa:
shapefile: mexico.shp_wgs84
data.frame: dat_mean_lab
mexico.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white",color="grey30")+ # Shapefile
geom_point(data=dat_mean_lab, aes(x=long,y=lat, color=Val.mean), size = 3)+ # data.frame
coord_fixed(ratio = 1) + # Colores
scale_color_gradient(low = "blue", high = "orange") +
labs(x="Longitud",y="Latitud",color="Problación promedio", title = "Población promedio en cada estado, 1990-2020.")+
theme_minimal()
De acuerdo con el gráfico anterior, los estados con más población son
los de la región Centro.
4.3 Polígonos de colores
mexico.shp_wgs84_d.f<-ggplot2::fortify(mexico.shp_wgs84, region="CVE_ENT") # Conversión de SpatialPolygonsDataFrame a data.frame
dim(mexico.shp_wgs84_d.f)-> [1] 750172 7head(mexico.shp_wgs84_d.f)-> long lat order hole piece id group
-> 1 -102.2879 22.41649 1 FALSE 1 01 01.1
-> 2 -102.2875 22.41608 2 FALSE 1 01 01.1
-> 3 -102.2870 22.41613 3 FALSE 1 01 01.1
-> 4 -102.2867 22.41639 4 FALSE 1 01 01.1
-> 5 -102.2865 22.41666 5 FALSE 1 01 01.1
-> 6 -102.2863 22.41690 6 FALSE 1 01 01.1mexico.shp_wgs84_d.f$CVE_ENT<-mexico.shp_wgs84_d.f$id # Creación de la variable CVE_ENT
head(mexico.shp_wgs84_d.f)-> long lat order hole piece id group CVE_ENT
-> 1 -102.2879 22.41649 1 FALSE 1 01 01.1 01
-> 2 -102.2875 22.41608 2 FALSE 1 01 01.1 01
-> 3 -102.2870 22.41613 3 FALSE 1 01 01.1 01
-> 4 -102.2867 22.41639 4 FALSE 1 01 01.1 01
-> 5 -102.2865 22.41666 5 FALSE 1 01 01.1 01
-> 6 -102.2863 22.41690 6 FALSE 1 01 01.1 01Para la generación del gráfico se ocupa:
data.frame: dat_mean_lab
El data.frame de un shapefile: mexico.shp_wgs84_d.f
# Generación del gráfico
ggplot(dat_mean_lab, aes(map_id = CVE_ENT)) + # map_id es la columna con la clave de unión. Nos ahorra el join previo.
geom_map(aes(fill =Val.mean),map = mexico.shp_wgs84_d.f) + #map= objeto SPDF convertido a data.frame.
expand_limits(x =mexico.shp_wgs84_d.f$long, y = mexico.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="Pob. promedio")+
theme_bw()
Gráfico coloreado por rangos:
Val.intervalo: variable creada con las función
classIntervals()data.frame: dat_mean_lab
Shapefile convertido a data.frame: mexico.shp_wgs84_d.f
#--- Creación de los intervalos
ngrupos<- 3
intervalo <- classIntervals(dat_mean_lab$Val.mean, 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_lab <- dat_mean_lab%>%
mutate(Val.intervalo = cut(Val.mean, punt.cor,dig.lab=8)) # con dig.lab se indican el número se cifras
#--- Generación del gráfico
ggplot(dat_mean_lab, # datos,
aes(map_id = CVE_ENT)) +# 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 = mexico.shp_wgs84_d.f) + # map= objeto SPDF convertido a data.frame.
scale_fill_brewer(palette = "red") +
expand_limits(x =mexico.shp_wgs84_d.f$long, y = mexico.shp_wgs84_d.f$lat)+ #Especificamos el tamaño del mapa al máximo de lat y long.
labs(x="Longitud",y="Latitud", fill="Población promedio", title = "Población promedio en cada estado, 1990-2020")+
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)))
4.4 Grafico con etiquetas
Para el siguiente gráfico se ocupa:
shapefile: mexico.shp_wgs84
data.frame de etiquetas: lab_mexico
mexico.shp_wgs84%>%
ggplot()+
geom_polygon(aes(x=long,y=lat,group=group), fill="white", color="grey30")+ # Shapefile
geom_text(data=lab_mexico, aes(long, lat, label =lab_NOMGEO), size=1.5)+ # data.frame etiq
labs(x="Longitud",y="Latitud", size=0.5)+ # Nombres ejes
theme_minimal()
Para crear el siguiente gráfico se ocupa:
El data.frame de los datos: dat_mean_lab
El data.frame de las etiquetas: lab_mexico
El shapefile convertido a data.frame: mexico.shp_wgs84_d.f
#--- Creación de los intervalos: Se divide la muestra en 4 grupos
ngrupos<- 4
intervalo <- classIntervals(dat_mean_lab$Val.mean, 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_lab <- dat_mean_lab%>%
mutate(Val.intervalo = cut(Val.mean, punt.cor,dig.lab=8)) # con dig.lab se indican el número de cifras
#--- Generación del gráfico
ggplot(dat_mean_lab, # datos,
aes(map_id = CVE_ENT)) +# 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 = mexico.shp_wgs84_d.f) + # map= objeto SPDF convertido a data.frame.
scale_fill_brewer(palette = "red") +
expand_limits(x =mexico.shp_wgs84_d.f$long, y = mexico.shp_wgs84_d.f$lat)+ # Especificamos el tamaño del mapa al máximo de lat y long.
geom_text(data=lab_mexico, aes(long, lat, label =lab_NOMGEO), size=1.5)+ # data.frame label
labs(x="Longitud",y="Latitud", fill="Población promedio", title = "Población promedio en cada estado, 1990-2020")+
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)))
#--- Creación de los intervalos: Se fijan los valores del intervalo
intervalo <- classIntervals(dat_mean_lab$Val.mean, fixedBreaks=c(500000,5000000,10000000,15000000),style="fixed")
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_lab <- dat_mean_lab%>%
mutate(Val.intervalo = cut(Val.mean, punt.cor,dig.lab=8)) # con dig.lab se indican el número de cifras
#--- Generación del gráfico
ggplot(dat_mean_lab, # datos,
aes(map_id = CVE_ENT)) +# 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 = mexico.shp_wgs84_d.f) + # map= objeto SPDF convertido a data.frame.
scale_fill_brewer(palette = "red") +
expand_limits(x =mexico.shp_wgs84_d.f$long, y = mexico.shp_wgs84_d.f$lat)+ # Especificamos el tamaño del mapa al máximo de lat y long.
geom_text(data=lab_mexico, aes(long, lat, label =lab_NOMGEO), size=1.5)+ # data.frame label
labs(x="Longitud",y="Latitud", fill="Población promedio", title = "Población promedio en cada estado, 1990-2020")+
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)))
4.5 Gráfico de regiones
Creación de una variable Region8 que contiene las 8 regiones de México:
Region8<-mexico.shp@data$CVE_ENT %>%
fct_collapse(Noroeste=Nor.Oeste_cla,
Noreste=Nor.Este_cla,
Oeste=Oeste_cla,
Oriente=Oriente_cla,
Centronorte=Centro.Norte_cla,
Centrosur=Centro.Sur_cla,
Suroeste=Sur.Oeste_cla,
Sureste=Sur.Este_cla) # Clave de la entidad
dat_1990_2020$Region8<-Region8 # Asignación de la variable región a la base de datos.
head(dat_1990_2020)-> CVE_ENT state_name pop.1990 pop.1995 pop.2000 pop.2005 pop.2010
-> 1 01 Aguascalientes 719659 862720 944285 1065416 1184996
-> 2 02 Baja California 1660855 2112140 2487367 2844469 3155070
-> 3 03 Baja California Sur 317764 375494 424041 512170 637026
-> 4 04 Campeche 535185 642516 690689 754730 822441
-> 5 05 Coahuila de Zaragoza 1972340 2173775 2298070 2495200 2748391
-> 6 06 Colima 428510 488028 542627 567996 650555
-> pop.2015 pop.2020 Region8
-> 1 1312544 1425607 Centronorte
-> 2 3315766 3769020 Noroeste
-> 3 712029 798447 Noroeste
-> 4 899931 928363 Sureste
-> 5 2954915 3146771 Noreste
-> 6 711235 731391 OestePara el gráfico de las regiones, se ocupa:
data.frame: dat_1990_2020, con la variable de interés: Region8.
Un shapefle convertido a data.frame: mexico.shp_wgs84_d.f
ggplot(dat_1990_2020,
aes(map_id=CVE_ENT))+
geom_map(aes(fill=Region8), # Variable a graficar
map=mexico.shp_wgs84_d.f, color="grey30")+
expand_limits(x=mexico.shp_wgs84_d.f$long,
y=mexico.shp_wgs84_d.f$lat)+
labs(x="Longitud",y="Latitud",fill="Regiones",title = "Regiones del 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)
)
5 Apéndice
5.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 mouse, 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
5.2 Paquetes de R
El paquete rgdal, se puede emplear para para leer y escribir información vectorial; maneja los objetos:
SpatialPointsDataFrame
SpatialLinesDataFrame
SpacialPolygonsDataFrame
5.3 Funciones de R
- La función
coordinates()transforma un data.frame en un objeto espacial.
5.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
6 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 de datos: http://datos.imss.gob.mx/dataset/informacion-en-salud/resource/5c865fa6-f632-48e4-8858-54d30350f107
Cancer http://datos.imss.gob.mx/dataset/informacion-en-salud/resource/6e5388e0-a19a-4893-ba30-576e91f21aea
Enlaces para análisis espaciales:
Geospatial Data Science in R: https://zia207.github.io/geospatial-r-github.io/spatial-autocorrelation.html
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