Introducción a R

16 - 20 de julio de 2018

Antes de iniciar

Es importante acotar que para trabajar en este curso:
- El principal tipo de archivo a usar con el modelo vectorial: SHAPEFILE
- El principal tipo de archivo a usar con el modelo raster: IMÁGENES
- También se usarán sets de datos en formato .csv

Shapefiles

Algunas generalidades: (Adaptado de: ESRI Course "Teaching with GIS").

Extensión .shp
Almacenan una colección de entidades con el mismo tipo de geometría
Además de las entidades también muestran sus atributos.
Un shapefile se compone de al menos tres archivos:
- .shp
- .shx
- .dbf
Si el shapefile está proyectado tendrá un archivo .prj
Se ha convertido en un estándar de facto

¿Quieres saber más? click acá

Shapefiles

Geometría de un shapefile:

Atributos

Atributos de un shapefile:

Preparando el workspace

# Definir el directorio de trabajo
setwd("/home/pau1esteban/R/")
# Verificar directorio
getwd()

## [1] "/home/pau1esteban/R"

# Instalar la librería rgdal desde el menú respectivo
# Cargar la librería para utilizarla
library(rgdal)

## Loading required package: sp

## rgdal: version: 1.3-3, (SVN revision 759)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.3-1

Cargando shapefiles en R

# Asignar el shapefile a la variable "capa"
capa <- readOGR(dsn = "./shp_curso/lim_urb/GEO_LIMITE_URBANO_CUENCA.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/pau1esteban/R/shp_curso/lim_urb/GEO_LIMITE_URBANO_CUENCA.shp", layer: "GEO_LIMITE_URBANO_CUENCA"
## with 1 features
## It has 6 fields

# Visualizar los nombres de atributos
names(capa)

## [1] "FID_PARROQ" "ID"         "NOMBRE"     "X_COORD"    "Y_COORD"   
## [6] "area_cue"

Visualizando shapefiles en R

# Visualización simple del shapefile "capa"
plot(capa)

Cargando archivos .csv en R

# Asignar el archivo a la variable puntos con read.csv() 
puntos <- read.csv("./at.csv")
# Verificar la estructura de "puntos"
str(puntos)

## 'data.frame':    130 obs. of  11 variables:
##  $ CODIGO   : Factor w/ 130 levels "002_AT_AQ_18",..: 18 19 20 21 22 23 24 25 26 27 ...
##  $ FECHA    : Factor w/ 31 levels "01/12/2017","02/12/2017",..: 3 3 2 2 2 1 4 5 2 2 ...
##  $ COORD_X  : Factor w/ 128 levels "673202,392214",..: 112 13 119 69 7 45 123 63 33 55 ...
##  $ COORD_Y  : Factor w/ 129 levels "9662603,9001",..: 96 66 80 18 124 36 85 47 38 99 ...
##  $ ZONA     : int  1 2 1 1 2 1 1 1 1 1 ...
##  $ HORA     : Factor w/ 115 levels "1:00:00 AM","1:00:00 PM",..: 60 66 57 36 24 71 72 6 8 40 ...
##  $ DATE     : Factor w/ 129 levels "01/12/2017 0:30:00",..: 22 20 9 10 13 3 23 25 11 14 ...
##  $ DIA      : Factor w/ 7 levels "domingo","jueves",..: 1 1 6 6 6 7 3 4 6 6 ...
##  $ TIPO     : Factor w/ 5 levels "ATROPELLO","CAIDA",..: 3 1 2 4 3 3 1 3 3 3 ...
##  $ MES      : Factor w/ 1 level "diciembre": 1 1 1 1 1 1 1 1 1 1 ...
##  $ AFECTADOS: int  3 1 1 2 3 2 1 1 4 2 ...

¿Como evitar que los atributos sean factor?

Cargando archivos .csv en R

# Leer archivo .csv sin que los caracteres se transformen en factor 
puntos <- read.csv("./at.csv",stringsAsFactors = F)
# Verificar la nueva estructura de "puntos"
str(puntos)

## 'data.frame':    130 obs. of  11 variables:
##  $ CODIGO   : chr  "1599_AT_AQ_17" "1600_AT_AQ_17" "1601_AT_AQ_17" "1602_AT_AQ_17" ...
##  $ FECHA    : chr  "03/12/2017" "03/12/2017" "02/12/2017" "02/12/2017" ...
##  $ COORD_X  : chr  "726113,17" "716199,671209" "726499,22" "722169,19" ...
##  $ COORD_Y  : chr  "9680534,71" "9679687,22604" "9680114,57" "9677281,87" ...
##  $ ZONA     : int  1 2 1 1 2 1 1 1 1 1 ...
##  $ HORA     : chr  "4:53:00 AM" "5:51:00 PM" "4:49:00 PM" "2:01:00 AM" ...
##  $ DATE     : chr  "03/12/2017 4:53:00" "03/12/2017 17:51:00" "02/12/2017 16:49:00" "02/12/2017 2:01:00" ...
##  $ DIA      : chr  "domingo" "domingo" "sábado" "sábado" ...
##  $ TIPO     : chr  "CHOQUE" "ATROPELLO" "CAIDA" "COLISION" ...
##  $ MES      : chr  "diciembre" "diciembre" "diciembre" "diciembre" ...
##  $ AFECTADOS: int  3 1 1 2 3 2 1 1 4 2 ...

Verificando atributos

Práctica:

Del atributo AFECTADOS obtener:
- El valor máximo
- El valor mínimo
- La media

Verificando atributos

# Valor máximo
max(puntos$AFECTADOS)

## [1] 15

# Valor mínimo
min(puntos$AFECTADOS)

## [1] 1

# Media
mean(puntos$AFECTADOS)

## [1] 2.6

La forma óptima de resolver esto es …

Preparando información geográfica

# Coerción de atributos caracter a numérico (Coordenadas)
puntos$COORD_X <- as.numeric(gsub(',','\\.',puntos$COORD_X))
puntos$COORD_Y <- as.numeric(gsub(',','\\.',puntos$COORD_Y))
# Comprobar la coerción
str(puntos)

## 'data.frame':    130 obs. of  11 variables:
##  $ CODIGO   : chr  "1599_AT_AQ_17" "1600_AT_AQ_17" "1601_AT_AQ_17" "1602_AT_AQ_17" ...
##  $ FECHA    : chr  "03/12/2017" "03/12/2017" "02/12/2017" "02/12/2017" ...
##  $ COORD_X  : num  726113 716200 726499 722169 699640 ...
##  $ COORD_Y  : num  9680535 9679687 9680115 9677282 9691822 ...
##  $ ZONA     : int  1 2 1 1 2 1 1 1 1 1 ...
##  $ HORA     : chr  "4:53:00 AM" "5:51:00 PM" "4:49:00 PM" "2:01:00 AM" ...
##  $ DATE     : chr  "03/12/2017 4:53:00" "03/12/2017 17:51:00" "02/12/2017 16:49:00" "02/12/2017 2:01:00" ...
##  $ DIA      : chr  "domingo" "domingo" "sábado" "sábado" ...
##  $ TIPO     : chr  "CHOQUE" "ATROPELLO" "CAIDA" "COLISION" ...
##  $ MES      : chr  "diciembre" "diciembre" "diciembre" "diciembre" ...
##  $ AFECTADOS: int  3 1 1 2 3 2 1 1 4 2 ...

Generando Data Frames ESPACIALES

# Instalar la librería raster desde el menu Tools, luego cargarla
library(raster)

# Verificar el indice de los atributos de coordenadas
names(puntos[,3:4])

## [1] "COORD_X" "COORD_Y"

# Generar el data frame proyectado con base en el archivo .csv
puntos_coord <- SpatialPointsDataFrame(puntos[,3:4], # Columnas de coords.
                puntos,  # Variable que contiene al .csv
                proj4string = crs("+init=epsg:32717")) # Proyección

Generando Data Frames ESPACIALES

# Verificar el tipo de dato de la nueva variable
class(puntos_coord)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

# Verificar información de proyección
#Forma 1
crs(puntos_coord)

## CRS arguments:
##  +init=epsg:32717 +proj=utm +zone=17 +south +datum=WGS84 +units=m
## +no_defs +ellps=WGS84 +towgs84=0,0,0

#Forma 2
proj4string(puntos_coord)

## [1] "+init=epsg:32717 +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Graficando Data Frames ESPACIALES

# Grafica simple del data frame espacial
plot(puntos_coord)

Algunas propiedades de plot

Se aplican a la simbología:

pch: Forma del símbolo
cex: Tamaño del símbolo
col: Color del símbolo
main: Título del gráfico

Algunas propiedades de plot

La propiedad pch tiene diversas formas:

Fuente: "R plot pch symbols" (STHDA).

También permite uso de caracteres como . * y otros

Personalizando plots

pch recibe el valor numérico del símbolo o un caracter
cex es un valor numérico (no necesariamente entero)
col es el color como tal (blue, green, red …)
El título debe ser acorde a lo que presentaremos y no muy extenso

UNA SIMBOLOGIA APROPIADA EXPRESA MEJOR LOS RESULTADOS

Mapeando con plots

# Aplicación de propiedades de plot: cex como caracter
plot(puntos_coord, pch=".", cex=2, col="red", main="Accidentes registrados")

Mapeando con plots

# Aplicación de propiedades de plot: cex como símbolo
plot(puntos_coord, pch=20, cex=0.5, col="blue", main="Accidentes")

Práctica:

Graficar el conjunto de datos anterior con:
- Un cuadrado rojo como simbología base
- Un tamaño de símbolo con valor entre 1 y 2
- Un título acorde a lo realizado

Combinando capas

# Cargar un shapefile existente
rurales <- readOGR(dsn="./shp_curso/p_rur/GEO_PARROQUIAS_RURALES_CUENCA.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/pau1esteban/R/shp_curso/p_rur/GEO_PARROQUIAS_RURALES_CUENCA.shp", layer: "GEO_PARROQUIAS_RURALES_CUENCA"
## with 21 features
## It has 4 fields

Para graficar dos capas es necesario usar la propiedad ADD (de plot)

Esta propiedad acepta valores T o F
Debe aplicarse a la segunda capa que se graficará

Combinando capas

# Graficar la capa de parroquias rurales
plot(rurales)
# Graficar los puntos de accidentes
plot(puntos_coord, pch=20, cex=0.5, col="blue", main="Accidentes", add = T)

Creando nuevas capas

# Transformar el Spatial Data Frame en Shapefile
writeOGR(obj=puntos_coord, #  Nuestro Spatial Data Frame
         dsn=".",  # Ubicación del directorio de trabajo 
         layer="puntos_nuevo",  # Nombre del nuevo shapefile
         driver="ESRI Shapefile",  # Tipo de archivo
         overwrite_layer = T, check_exists = T) # Verificación

Práctica:

Graficar el shapefile exportado anteriormente
Combinarlo con la capa de parroquias rurales
Simbolizarlo adecuadamente.
¿Qué sucede si no establezco el valor de la propiedad cex?

Resultado

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/pau1esteban/R/puntos_nuevo.shp", layer: "puntos_nuevo"
## with 130 features
## It has 11 fields

Trabajando con atributos

Un aspecto esencial de S.I.G. es la posibilidad de trabajar con atributos.
Los atributos son las características de una entidad:
- nombre, edad, estatura …
¿Atributos geográficos?
Muchas veces será necesario seleccionar ciertos atributos
Para ello usaremos la librería dplyr y el mismo set de datos

Preparando los datos

# Cargar la libreria dplyr 
library(dplyr)
# Verificar los atributos del archivo
str(puntos)

## 'data.frame':    130 obs. of  11 variables:
##  $ CODIGO   : chr  "1599_AT_AQ_17" "1600_AT_AQ_17" "1601_AT_AQ_17" "1602_AT_AQ_17" ...
##  $ FECHA    : chr  "03/12/2017" "03/12/2017" "02/12/2017" "02/12/2017" ...
##  $ COORD_X  : num  726113 716200 726499 722169 699640 ...
##  $ COORD_Y  : num  9680535 9679687 9680115 9677282 9691822 ...
##  $ ZONA     : int  1 2 1 1 2 1 1 1 1 1 ...
##  $ HORA     : chr  "4:53:00 AM" "5:51:00 PM" "4:49:00 PM" "2:01:00 AM" ...
##  $ DATE     : chr  "03/12/2017 4:53:00" "03/12/2017 17:51:00" "02/12/2017 16:49:00" "02/12/2017 2:01:00" ...
##  $ DIA      : chr  "domingo" "domingo" "sábado" "sábado" ...
##  $ TIPO     : chr  "CHOQUE" "ATROPELLO" "CAIDA" "COLISION" ...
##  $ MES      : chr  "diciembre" "diciembre" "diciembre" "diciembre" ...
##  $ AFECTADOS: int  3 1 1 2 3 2 1 1 4 2 ...

Selección simple

# Seleccionar las coordenadas,tipo y cantidad de afectados
seleccion <- select(puntos, COORD_X, COORD_Y, TIPO, AFECTADOS)
# Mostrar los primeros registros de la selección
head(seleccion)

##    COORD_X COORD_Y      TIPO AFECTADOS
## 1 726113.2 9680535    CHOQUE         3
## 2 716199.7 9679687 ATROPELLO         1
## 3 726499.2 9680115     CAIDA         1
## 4 722169.2 9677282  COLISION         2
## 5 699639.8 9691822    CHOQUE         3
## 6 720471.3 9678478    CHOQUE         2

Selección basada en un criterio (Filtrar)

# Mostrar las coordenadas de los accidentes con tres o más afectados
puntos %>%
  filter(AFECTADOS >= 3) %>%
  select(COORD_X,COORD_Y,AFECTADOS)

##     COORD_X COORD_Y AFECTADOS
## 1  726113.2 9680535         3
## 2  699639.8 9691822         3
## 3  719517.6 9678733         4
## 4  720628.8 9679777        11
## 5  729507.4 9672370         4
## 6  678899.0 9693300         6
## 7  720502.0 9677929         3
## 8  723224.9 9677781         3
## 9  721096.9 9682045         4
## 10 721862.9 9681468         3
## 11 717873.8 9676660         5
## 12 714654.5 9683600         6
## 13 724350.1 9679451         3
## 14 726281.5 9680386         3
## 15 723988.6 9680961         3
## 16 718700.2 9680097         4
## 17 721890.0 9678146         5
## 18 717149.9 9675093        15
## 19 721876.1 9677906         3
## 20 718662.3 9684875         4
## 21 726140.9 9683267         4
## 22 720592.9 9678278         3
## 23 726113.1 9681436         4
## 24 719382.9 9679544         3
## 25 720344.2 9680146         6
## 26 721409.5 9679482         4
## 27 721550.9 9680257         5
## 28 728310.2 9674693         4
## 29 722813.8 9678913         3
## 30 723140.8 9680556         3
## 31 716293.3 9681034         5
## 32 725249.3 9679234        11
## 33 716723.4 9680904         9
## 34 719097.6 9680337         4
## 35 714846.9 9678193         3
## 36 726193.6 9683360         3
## 37 721328.8 9680780         3
## 38 724866.9 9679119         4
## 39 720862.9 9677796         3
## 40 727019.4 9687840         3
## 41 725352.4 9679465         4
## 42 718197.8 9676977         3
## 43 718868.6 9677440         6
## 44 719067.8 9677702         3
## 45 724019.7 9679335         4
## 46 723731.6 9681193         3
## 47 726306.1 9683574         4
## 48 722372.9 9677290         3
## 49 717681.0 9664679         4

Selección de atributos únicos

# Seleccionar los tipos de accidente registrados
tipos_accid <- distinct(puntos,TIPO)
# Mostrar la selección
tipos_accid

##          TIPO
## 1      CHOQUE
## 2   ATROPELLO
## 3       CAIDA
## 4    COLISION
## 5 VOLCAMIENTO

Selección de atributos únicos

# Seleccionar la cantidad de accidentes por tipo
tipos_accid <- distinct(count(puntos,TIPO))
# Mostrar la selección
tipos_accid

## # A tibble: 5 x 2
##   TIPO            n
##   <chr>       <int>
## 1 ATROPELLO      26
## 2 CAIDA           7
## 3 CHOQUE         84
## 4 COLISION        5
## 5 VOLCAMIENTO     8