Uso del paquete sf

Se describe el uso del paquete sf para el uso de datos espaciales

Author

Affiliation

Yandira Cuvero

Escuela Politécnica Nacional

Datos espaciales se puede representar como vectores de información que pueden representar pontos, rectas y polígonos. Además, estos vectores pueden tener columnas adicionales que contenga información adicional para describir el proceso. El paquete sf permite trabajar con datos espaciales vectoriales.

Code

library(sf)

Shapefiles

Vectores con información espacial se suelen almacenar en archivos de formato shapefile. Para poder describir completamente un archivo con datos espaciales, este formato incluye varias archivos, que en conjunto los describen. Los archivos principales son:

• .shp tiene la información geométrica de los datos.

• .shx indica el índice pocisional de los datos geométricos que permite manejar los datos del archivo .shp

• .dbf contiene los atributos de cada forma.

  Otros archivos que también se pueden incluir son:

• .prj describe la proyección de los datos.

• .sbn y .sbx otros índices para describir los datos, y

• .shp.xml que tiene los meta datos en formato XML.

  Entonces, al trabajar con shapefiles, es importante recordar que no es aconsejable trabajar solo con el archivo .shp sino con los archivos antes mencionados.

Para leer estos archivos usamos la función st_read() del paquete sf se puede utilizar para cargar shapefiles.

Cargar la base de datos

Cargamos la base de datos, como resultados revisamos el tipo y la clase de los datos cargados.

Code

# Colocar la dirección del archivo
#map <- st_read('EDUCATION.shp')
class(map)

[1] "sf"         "data.frame"

Revisamos la estructura de los datos obtenidos:

Code

str(map)

Classes 'sf' and 'data.frame':  4389 obs. of  19 variables:
 $ PERIOD    : chr  "2019-2020" "2019-2020" "2019-2020" "2019-2020" ...
 $ I_CODE    : chr  "17H01486" "17H01486" "17H01486" "17H01486" ...
 $ I_NAME    : chr  "SAN ANDRES" "SAN ANDRES" "SAN ANDRES" "SAN ANDRES" ...
 $ UTM_X     : num  776430 776430 776430 776430 776430 ...
 $ UTM_Y     : num  9960077 9960077 9960077 9960077 9960077 ...
 $ TYPE_OF_ED: chr  "Regular education" "Regular education" "Regular education" "Regular education" ...
 $ SUSTENANCE: chr  "Fiscal" "Fiscal" "Fiscal" "Fiscal" ...
 $ ZONE      : chr  "Rural" "Rural" "Rural" "Rural" ...
 $ REGIME    : chr  "Sierra" "Sierra" "Sierra" "Sierra" ...
 $ JURISDICTI: chr  "Hispana" "Hispana" "Hispana" "Hispana" ...
 $ SEX       : chr  "Feminine" "Feminine" "Feminine" "Feminine" ...
 $ GRADE     : chr  "1er año Básica" "2do año Básica" "3er año Básica" "4to año Básica" ...
 $ P_NPROMOTE: num  0 0 0 0 0 0 0 0 0 0 ...
 $ P_DESERTER: num  0 0 0 0 0 0 0 0 0 0 ...
 $ P_PROMOTED: num  100 100 100 100 100 100 100 100 100 100 ...
 $ TOTAL     : num  100 100 100 100 100 100 100 100 100 100 ...
 $ FID_DPA_SE: num  4079 4079 4079 4079 4079 ...
 $ DPA_SECTOR: chr  "170150999033" "170150999033" "170150999033" "170150999033" ...
 $ geometry  :sfc_POINT of length 4389; first list element:  'XY' num  776430 9960077
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
  ..- attr(*, "names")= chr [1:18] "PERIOD" "I_CODE" "I_NAME" "UTM_X" ...

Proyecciones

Revisamos el tipo de proyección disponible, esto será fundamental en el caso que se requiera modificar la proyección, ya sea en el mapa o en los puntos considerados.

Code

st_crs(map)

Coordinate Reference System:
  User input: WGS 84 / UTM zone 17S 
  wkt:
PROJCRS["WGS 84 / UTM zone 17S",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 17S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-81,
            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",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",32717]]

Gráficos

Para gráficar todas las variables usamos:

Code

plot(map[,8:18]) # plot first attribute

Para una sola variable

Code

plot(map['P_PROMOTED'])

Una forma mejor es usando el paquete ggplot

Code

library(ggplot2)
ggplot(map) + geom_sf(aes(color = GRADE)) + theme_bw()

Una forma mejor es usando el paquete ggplot

Code

library(plotly)
g <- ggplot(map) + geom_sf(aes(color = JURISDICTI))
ggplotly(g)

Maps

Determinamos la proyección de nuestros datos

Code

library(mapview)
mapview(map)

División política

Code

# Colocar la dirección del archivo
#map <- st_read('nxparroquias.shp')

De esto, vamos a tomar únicamente en Pichincha y de esta tomar únicamente la Parroquia Quito, para ello tomamos:

Code

p <-  dp[dp$DPA_PROVIN == 17, ]
cq <- p[p$DPA_DESCAN == "QUITO", ]
q <- cq[cq$DPA_DESPAR=="QUITO" ,]

plot(cq)

Dividimos en una grilla esta área

Code

grid <- q %>%
  st_make_grid(cellsize = c(1000,1000)) %>%
  st_intersection(q) %>% 
   st_cast("MULTIPOLYGON") %>%
  st_sf()
plot(grid)

Incluimos a la vez los puntos y la grilla generada

Code

points <- map
ggplot() + geom_sf(data = grid) + geom_sf(data = points)

Conteo

Code

inter <- st_intersects(grid, points)
grid$count <- lengths(inter)
ggplot(grid) + geom_sf(aes(fill = count))

Usando el comando scale_fill_continuous podemos obtener una mejor representación de los colores.

Code

ggplot(grid) + geom_sf(aes(fill = count))+
scale_fill_continuous(type = "viridis")

Selección de zonas

Finalmente, podemos identificar las zonas que tienen conteo no nulo. En general, para realizar cualquier análisis se requiere que el área con la que se está trabajado tengo menos del 30% de observaciones nulas.

Code

grid.selec <- grid[grid$count>0, ]
ggplot(grid.selec) + geom_sf(aes(fill = count))+
geom_sf_text(aes(label = count), size = 2, color = "white")+
scale_fill_continuous(type = "viridis") +
theme(
  axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
  axis.title.x = element_blank(),
  axis.title.y = element_blank())

Agregamos las celdas nulas faltantes en la grilla anterior:

Code

ref<- as.numeric(rownames(grid.selec))
relle <- function(x){
  if(x[2]<4){
    r <- as.character(seq(x[1],x[1]+x[2], by=1))
  }else{
    r <- as.character(x[1])
  }
  return(r)
}
ix <- data.frame(ref, lag = c(diff(ref)-1, 0))
sl <- as.character(unlist(apply(ix,1, relle)))

Code

grid2 <- grid[sl, ]
grid2$sl <- sl
ggplot(grid2) + geom_sf(aes(fill = count))+
geom_sf_text(aes(label = sl), size = 1.5, color = "white")+
theme(
  axis.text.x = element_blank(),
  axis.title.x = element_blank(),
  axis.title.y = element_blank())

Finalmente, agregamos y eliminamos manualmente las observaciones que han quedado disjuntas al área de análisis

Code

f <- c(sl[-1], as.character(c(29, 76, 79, 80, 81, 269, 270,
                              425, 426, 427, 428, 429, 430 )))
grid3 <- grid[f, ]
ggplot(grid3) + geom_sf(aes(fill = count))+
geom_sf_text(aes(label = count), size = 2, color = "white")+
scale_fill_continuous(type = "viridis") +
theme(
  axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
  axis.title.x = element_blank(),
  axis.title.y = element_blank())

Referencias

• https://www.ecuadorencifras.gob.ec/documentos/web-inec/Geografia_Estadistica/Micrositio_geoportal/index.html
• Moraga, Paula. (2023). Spatial Statistics for Data Science: Theory and Practice with R. Chapman & Hall/CRC. <https://www.paulamoraga.com/book-spatial/spatial-data-in-r.html>
• https://samanthacsik.github.io/posts/2022-10-24-quarto-blogs/
• https://quarto.org/docs/websites/website-listings.html
• https://quarto.org/docs/output-formats/html-themes.html#sass-variables