Introducción a R

16 - 20 de julio de 2018

Acceder a datos de shapefiles

Esto se realiza mediante slots
Hemos usado $ para acceder a datos de atributos
Para acceder a datos geográficos usamos el slot @
Principalmente podremos visualizar:
- Los atributos del elemento geográfico (data)
- Los valores extremos del elemento (bbox)
- Los datos de la proyección

Acceder a datos de shapefiles

# Cargar library rgdal
library(rgdal)
# Establecer el directorio de trabajo
setwd("/home/pau1esteban/R/")
# Cargar el shapefile de parroquias rurales
capa <- 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

# Visualizar información de proyección
capa@proj4string

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

Acceder a datos de shapefiles

# Visualizar atributos
capa@data

##              DESPAR__18      AREA  Area_Km2 RH_2017
## 0                  BAOS 326808484 326.80848      38
## 1               CHAUCHA 313048717 313.04872       0
## 2                 CHECA  62560593  62.56059       2
## 3            CHIQUINTAD  92948541  92.94854       7
## 4                 CUMBE  71345316  71.34532       7
## 5              EL VALLE  43204717  43.20472      63
## 6                LLACAO  16998261  16.99826      11
## 7             MOLLETURO 977343857 977.34386      16
## 8                 NULTI  26953174  26.95317       8
## 9       OCTAVIO CORDERO  20686477  20.68648       5
## 10               PACCHA  25731189  25.73119      14
## 11              QUINGEO 117594927 117.59493       5
## 12             RICAURTE  13834616  13.83462      61
## 13          SAN JOAQUIN 188631203 188.63120      22
## 14            SANTA ANA  43109447  43.10945       7
## 15              SAYAUSI 365532286 365.53229      24
## 16               SIDCAY  16663877  16.66388       7
## 17             SININCAY  24701639  24.70164      16
## 18               TARQUI 137370019 137.37002      12
## 19                 TURI  26653082  26.65308      65
## 20 VICTORIA DEL PORTETE 201698486 201.69849       3

Acceder a datos de shapefiles

# Visualizar atributos específicos
capa@data$AREA

##  [1] 326808484 313048717  62560593  92948541  71345316  43204717  16998261
##  [8] 977343857  26953174  20686477  25731189 117594927  13834616 188631203
## [15]  43109447 365532286  16663877  24701639 137370019  26653082 201698486

# Visualizar valores extremos
capa@bbox

##         min       max
## x  652796.1  739204.5
## y 9648867.0 9723945.0

Insertar mapas base

Es posible insertar mapas base alojados en servidores
Para ellos es necesario utilizar la librería ggmap
Visualizaremos un mapa de vías y de imagen satelital
Combinaremos una capa de puntos con un mapa base

Cargar un archivo con CRS de latitud/longitud

# Cargar el archivo .csv con CRS Lat/Long
datos <- read.csv("./datos.csv", stringsAsFactors = F, sep="")
# Verificar estructura de los datos
str(datos)

## 'data.frame':    94 obs. of  10 variables:
##  $ Codigo   : int  2017 2017 2017 2017 2017 2017 2017 2017 2017 2017 ...
##  $ Anio     : chr  "Enero" "Enero" "Enero" "Enero" ...
##  $ Mes      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Semana   : chr  "Domingo" "Domingo" "Domingo" "Domingo" ...
##  $ Dia      : chr  "01/ene/2017" "01/ene/2017" "01/ene/2017" "01/ene/2017" ...
##  $ Fecha    : chr  "01/01/2017" "01/01/2017" "01/01/2017" "01/01/2017" ...
##  $ Time     : chr  "00:08:14" "00:15:03" "00:13:13" "00:19:08" ...
##  $ Longitude: chr  "-79,045527" "-79,009777" "-79,028315" "-79,064" ...
##  $ Latitude : chr  "-2,879941" "-2,960284" "-2,866686" "-2,919852" ...
##  $ Canton   : chr  "CUENCA" "CUENCA" "CUENCA" "CUENCA" ...

Transformación de tipos de datos

# Coerción de valores de coordenadas (caracter a numérico)
datos$Latitude <- as.numeric(gsub(',','\\.',datos$Latitude))
datos$Longitude <- as.numeric(gsub(',','\\.',datos$Longitude))
# Verificar estructura de los datos
str(datos)

## 'data.frame':    94 obs. of  10 variables:
##  $ Codigo   : int  2017 2017 2017 2017 2017 2017 2017 2017 2017 2017 ...
##  $ Anio     : chr  "Enero" "Enero" "Enero" "Enero" ...
##  $ Mes      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Semana   : chr  "Domingo" "Domingo" "Domingo" "Domingo" ...
##  $ Dia      : chr  "01/ene/2017" "01/ene/2017" "01/ene/2017" "01/ene/2017" ...
##  $ Fecha    : chr  "01/01/2017" "01/01/2017" "01/01/2017" "01/01/2017" ...
##  $ Time     : chr  "00:08:14" "00:15:03" "00:13:13" "00:19:08" ...
##  $ Longitude: num  -79 -79 -79 -79.1 -79 ...
##  $ Latitude : num  -2.88 -2.96 -2.87 -2.92 -2.88 ...
##  $ Canton   : chr  "CUENCA" "CUENCA" "CUENCA" "CUENCA" ...

Generación de data frame espacial

# Verificación de índice de atributos 
names(datos[,8:9])

## [1] "Longitude" "Latitude"

# Generar el data frame espacial con coordenadas Lat/Long
puntos_latlong <- SpatialPointsDataFrame(datos[,8:9], #8(Long) 9(Lat)
                  datos,  # Variable que contiene al .csv
                  proj4string = CRS("+init=epsg:4326")) #Proyección

# Generar un nuevo shapefile con base en el data frame espacial
writeOGR(obj=puntos_latlong,dsn=".",layer="capa_latlong",
         driver="ESRI Shapefile", overwrite_layer = T,check_exists = T)

Graficación de los datos

# Graficar el data frame espacial
plot(puntos_latlong, pch=20, cex=0.7, col="blue")

Reproyectar shapefiles

# Cargar librería raster
library(raster)
# Cargar shapefile en CRS UTM WGS84 Zona 17S
capa_parroquias <- readOGR(dsn="./shp_curso/p_urb/GEO_PARROQUIAS_URBANAS_CUENCA.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/pau1esteban/R/shp_curso/p_urb/GEO_PARROQUIAS_URBANAS_CUENCA.shp", layer: "GEO_PARROQUIAS_URBANAS_CUENCA"
## with 15 features
## It has 5 fields

# Verificar CRS
capa_parroquias@proj4string

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

# Verificar valores extremos
extent(capa_parroquias)

## class       : Extent 
## xmin        : 714264.1 
## xmax        : 734270.1 
## ymin        : 9675306 
## ymax        : 9687848

Reproyectar shapefiles

# Verificar valores extremos
extent(capa_parroquias)

## class       : Extent 
## xmin        : 714264.1 
## xmax        : 734270.1 
## ymin        : 9675306 
## ymax        : 9687848

Reproyectar shapefiles

# Cargar shapefile en CRS Geograficas WGS84
capa_puntos <- readOGR(dsn="./capa_latlong.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/pau1esteban/R/capa_latlong.shp", layer: "capa_latlong"
## with 94 features
## It has 10 fields

# Verificar CRS
crs(capa_puntos)

## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

# Verificar valores extremos
extent(capa_puntos)

## class       : Extent 
## xmin        : -79.39686 
## xmax        : -78.89776 
## ymin        : -3.086492 
## ymax        : -2.765934

Reproyectar shapefiles

# Verificar valores extremos
extent(capa_puntos)

## class       : Extent 
## xmin        : -79.39686 
## xmax        : -78.89776 
## ymin        : -3.086492 
## ymax        : -2.765934

Reproyectar shapefiles

# Reproyectar capa UTM en Lat/Long
parr_latlong <- spTransform(capa_parroquias, crs(capa_puntos))
# Verificar el CRS
crs(parr_latlong)

## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

Reproyectar shapefiles

# Graficar los resultados
plot(parr_latlong)
plot(capa_puntos,add=T, pch=20,col="red",cex=0.8)

Geoprocesamiento: Clasificar datos

Realizaremos un mapa de clasificación de datos de robos ocurridos en las parroquias rurales, siguiendo estos pasos:
- Instalar y cargar las librerías RColorBrewer, prettymapr, maps y mapdata
- Revisar la información de la capa respectiva
- Crear una paleta de colores e intervalos de clasificación
- Estructurar el mapa coherentemente
- Exportar el resultado como imagen PNG.

Manejo de datos Raster

Para manejar datos raster usamos la library raster
Recordar que los raster son imágenes

# Cargar la libreria
library(raster)
# Crear un elemento de tipo raster
mi_raster <- raster()
# Ver las propiedades del raster
mi_raster

## class       : RasterLayer 
## dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

Manejo de datos Raster

# Ver resolución del raster
res(mi_raster)

## [1] 1 1

# Asignar un CRS al raster
projection(mi_raster) <- "+init=epsg:32717"
# Verificar la proyección
mi_raster

## class       : RasterLayer 
## dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32717 +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

Manejo de datos Raster

Crearemos un raster con valores de celda aleatorios

# Asignar dimensiones al raster
mi_raster <- raster(ncol=10, nrow=10)
# Verificar cantidad de celdas del raster
ncell(mi_raster)

## [1] 100

# Generación de números aleatorios
set.seed(0)
# Asignar valores al raster mediante una distribución uniforme
values(mi_raster) <- runif(ncell(mi_raster))

Manejo de datos Raster

# Graficar el raster
plot(mi_raster, main="Raster")