Introducción a los datos SIG con R
Datos espacio temporales: Vectores
VE-CC
DataIntelligence12-05-2021
Abstract
En este articulo realizaremos consideraciones iniciales sobre analisis espacial y analizaremos los dos tipos de datos espacio-temporales fundamentales: vectoriales y rasters. La gran ventaja de usar R en lugar de un SIG es el de la posibilidad de automatizar procesosy leer conjuntos de datos difíciles o imposibles de gestionar con un SIG tradicional.
Índice
I Introducccion
1 Proyecciones
Todos los objetos espaciales tanto en R como en cualquier GIS necesitan un sistema de coordenadas de referencia (Coordinate reference system- CRS) que defina como se relaciona la informacion que contienen con la superficie de la tierra. Esta relacion, que en R es definida por el parametro CRS, puede ser de dos tipos: geografica (utilizan sistemas de medida angulares) o proyectada (asumen una superficie bidimensional basada en coordenadas cartesianas).
2 Transformaciones
Si trabajamos con dos objetos espaciales con diferentes sistemas de coordenadas, las operaciones de analisis espacial que queramos aplicar sobre ellos seran inútiles, ya que, a efectos prácticos, no tienen coincidencia espacial. Para solucionarlo, hay que cambiar el sistema de coordenadas de uno de ellos para que tenga coincidencia con el otro, o bien cambiar el de los dos por uno comun.
II Datos espacio-temporales I: vectorial
Para importar y manipular shapefiles, geojsons, geodatabases, etc, con R podemos utilizar los siguientes paquetes:
sp: Conjunto de clases para diferentes tipos de vectores y rasters. rgdal: importar y exportar datos espaciales. sf: nuevo paquete para importar y manipular datos vectoriales.
1 Los paquetes sp, rgdal y sf
Crearemos un data.frame con 100 coordenadas aleatorias entre -180 y +180 y entre -90 y +90 simulando 100 localizaciones sobre el planeta. Hacemos lo mismo para crear los correspondientes valores aleatorios. Luego, con coordinates() transformaremos nuestro data.frame en un objeto espacial.
La forma clásica para importar, exportar y manipular datos vectoriales se hace a traves de sp y rgdal, sin embargo esta el paquete nuevo disponible sf, el que permite:
i trabajar con una única clase para todos los tipos de datos vectoriales.
ii su estrutura de datos es del tipo dataframe
iii se instala mas fácilmente
iv es rapido y compatible con tidyverse.
Construyamos un dataframe con 100 datos aleatorios en 100 coordenadas aleatorias y trasformemoslo en un objeto espacial:
library(sp)## Warning: package 'sp' was built under R version 4.0.5set.seed(12345)
data <- data.frame(x = runif(100,-180,180),
y = runif(100,-90,90),
valores = runif(100,150,250))
coordinates(data) <- c("x","y")podemos obtener informacion espacial referida a los limites externos de nuestro data:
bbox(data)## min max
## x -179.59083 176.30530
## y -88.75274 88.80794y los representamos graficamente
spplot(data)
2 Importacion
En R, rgdal basicamente lo utilizaremos para para leer y escribir informacion geoespacial vectorial.
La funcion readOGR() se utiliza para leer informacion vectorial, por ejemplo un shp de ESRI, y carga la capa asignandole el tipo de datos de objeto espacial que le corresponde segun se trate de puntos, lineas o poligonos (SpatialPointsDataFrame, SpatialLinesDataFrame o SpacialPolygonsDataFrame, respectivamente).
library(rgdal)## rgdal: version: 1.5-12, (SVN revision 1018)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/usuario/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/usuario/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.paises <- readOGR("./libro_sig_r/05_spatial_I/world_countries.shp")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\usuario\Desktop\ds\ds_rgee\estudio\libro_sig_r\05_spatial_I\world_countries.shp", layer: "world_countries"
## with 255 features
## It has 94 fields#sss <- read_sf("./libro_sig_r/05_spatial_I/world_countries.shp")plot(paises)
# sssciudades <- readOGR("./libro_sig_r/05_spatial_I/world_Cities.shp")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\usuario\Desktop\ds\ds_rgee\estudio\libro_sig_r\05_spatial_I\World_Cities.shp", layer: "World_Cities"
## with 2540 features
## It has 13 fields
## Integer64 fields read as strings: FID ObjectID POP POP_RANK PORT_ID LABEL_FLAGplot(ciudades)
class(paises)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"class(ciudades)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"3 Visualizacion
A pesar de que estamos trabajando con objetos espaciales (diferentes que los dataframe, numeric o character, por ejemplo), podemos usar las funciones plot, lines y points para representar las capas que hemos cargado. Internamente R reconoce que esos objetos son de clase espacial y llama a las funciones correspondientes para representarlos.
plot(paises, col = 'blue', ylim = c(30,70), xlim = c(-25, 30), lwd = 0.5)
points(ciudades, col = ciudades@data$POP_RANK, cex = 0.9)
legend('bottomleft', pch = 1,
col =sort(unique(ciudades@data$POP_RANK)),
legend = sort(unique(ciudades@data$POP_RANK)),
bg = 'white', bty = 'white', title = 'Ranking de poblacion')
4 Proyecciones y coordenadas
Como objetos espaciales, los puntos y polígonos que hemos cargado tienen asociados unos atributos extra, además de la información alfanumérica correspondiente. Podemos extraer la información espacial de, por ejemplo, las coordenadas de las ciudades con la función coordinates(). Si lo aplicasemos sobre los países, nos devolvería las coordenadas de los centroides de cada uno de los polígonos que están representados en el objeto.
head(coordinates(ciudades))## coords.x1 coords.x2
## [1,] 99.5290 16.472997
## [2,] 12.1900 -5.556996
## [3,] 100.3490 16.439004
## [4,] 15.7833 -17.049996
## [5,] 24.8300 -28.659996
## [6,] 105.4380 10.382004head(coordinates(paises))## [,1] [,2]
## 0 114.01881 -0.1913581
## 1 114.72126 3.6147609
## 2 -71.24568 -35.4519491
## 3 -64.68475 -16.7068768
## 4 -74.37788 -9.1543228
## 5 -65.15516 -35.1865658Estas cordenadas estan asociadas a un sistema de proyeccion concreto, que necesitamos conocer (o definir) en cada caso. La funcion proj4string() nos devuelve la proyeccion de un objeto espacial.
proj4string(paises)## Warning in proj4string(paises): CRS object has comment, which is lost in output## [1] "+proj=longlat +datum=WGS84 +no_defs"Vemos que lo que muestra es una cadena de texto con una secuencia de parámetros importantes para la definición de la proyección: se trata de un sistema de codificación particular denominado proj4. Podemos obtener una tabla con el código EPSG y el proj4string de todas las proyecciones disponibles con la funcion make_EPSG() del paquete rgdal.
Tabla con el código EPSG y el proj4string de todas las proyecciones disponibles:
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.1 v dplyr 1.0.2
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()epsg <- make_EPSG()
head(epsg)## code note prj4 prj_method
## 1 3819 HD1909 +proj=longlat +ellps=bessel +no_defs +type=crs (null)
## 2 3821 TWD67 +proj=longlat +ellps=aust_SA +no_defs +type=crs (null)
## 3 3822 TWD97 +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs (null)
## 4 3823 TWD97 +proj=longlat +ellps=GRS80 +no_defs +type=crs (null)
## 5 3824 TWD97 +proj=longlat +ellps=GRS80 +no_defs +type=crs (null)
## 6 3887 IGRS +proj=geocent +ellps=GRS80 +units=m +no_defs +type=crs (null)# Seleccion de las proyecciones europeas
filter(epsg, str_detect(note, "Europe"))## code note
## 1 3034 ETRS89-extended / LCC Europe
## 2 3035 ETRS89-extended / LAEA Europe
## 3 3575 WGS 84 / North Pole LAEA Europe
## 4 5632 PTRA08 / LCC Europe
## 5 5633 PTRA08 / LAEA Europe
## 6 5634 REGCAN95 / LCC Europe
## 7 5635 REGCAN95 / LAEA Europe
## 8 5636 TUREF / LAEA Europe
## 9 5637 TUREF / LCC Europe
## 10 5638 ISN2004 / LAEA Europe
## 11 5639 ISN2004 / LCC Europe
## 12 9039 ISN2016 / LAEA Europe
## 13 9040 ISN2016 / LCC Europe
## prj4
## 1 +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m +no_defs +type=crs
## 2 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs +type=crs
## 3 +proj=laea +lat_0=90 +lon_0=10 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs
## 4 +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m +no_defs +type=crs
## 5 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs +type=crs
## 6 +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m +no_defs +type=crs
## 7 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs +type=crs
## 8 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs +type=crs
## 9 +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m +no_defs +type=crs
## 10 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs +type=crs
## 11 +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m +no_defs +type=crs
## 12 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs +type=crs
## 13 +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m +no_defs +type=crs
## prj_method
## 1 Lambert Conic Conformal (2SP)
## 2 Lambert Azimuthal Equal Area
## 3 Lambert Azimuthal Equal Area
## 4 Lambert Conic Conformal (2SP)
## 5 Lambert Azimuthal Equal Area
## 6 Lambert Conic Conformal (2SP)
## 7 Lambert Azimuthal Equal Area
## 8 Lambert Azimuthal Equal Area
## 9 Lambert Conic Conformal (2SP)
## 10 Lambert Azimuthal Equal Area
## 11 Lambert Conic Conformal (2SP)
## 12 Lambert Azimuthal Equal Area
## 13 Lambert Conic Conformal (2SP)Si queremos asignar una proyeccion a un objeto espacial que no tiene, debemos hacerlo con la funcion CRS(). Esto funciona tanto para objetos espaciales de clase vectorial (puntos, lineas y poligonos) compo para rasters.
# Asignamos proyeccion WGS84 en formato proj4 string
a_crs_object <- CRS("+proj=longlat +datum=WGS84 +no_defs")
class(a_crs_object)## [1] "CRS"
## attr(,"package")
## [1] "sp"a_crs_object## CRS arguments: +proj=longlat +datum=WGS84 +no_defsLa mayor parte de la veces no vamos a poder recordar la cadena de texto completa en formato proj4. Para facilitarlo, se puede usar en su lugar el codigo EPSG, una abreviatura numerica para cada una de las proyecciones disponibles.
La cadena de texto en este caso es +init=epsg:4326, siendo 4326 el codigo de la proyeccion que queremos asignar (en este caso WGS84).
Asignamos la proyeccion WGS84 en formato EPSG:
a_crs_object_epsg <- CRS("+init=epsg:4326")
class(a_crs_object_epsg)## [1] "CRS"
## attr(,"package")
## [1] "sp"a_crs_object_epsg## CRS arguments: +proj=longlat +datum=WGS84 +no_defsAhora que sabemos estraer la proyeccion de un objeto espacial y definirla si no la tenia, como la cambiamos? Aplicando una transformacion, que con R se hace con spTransform():
# Cambiar la proyeccion de un objeto espacial a Robinson
prj <- CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
paises_sptr <- spTransform(paises, prj)
plot(paises_sptr)
5 Acceso a las variables
Las variables, o campos de informacion, estan asociadas a todos los objetos espaciales. Para acceder a ellas podemos usar @, con la que accedemos a los diferentes slots o almacenes de informacion. En nuestro caso, el objeto paises tiene cinco slots:
slotNames(paises)## [1] "data" "polygons" "plotOrder" "bbox" "proj4string"donde:
data: tabla de atributo de la capa.
polygon: coordenadas de todos los vertices que componen casa uno de los poligonos de la capa(informacion espacial).
plotOrder: orden en el que se dibujan los poligonos de la capa.
bbox: coodenadas mas externas del marco que componen todos los objetos de la capa.
proj4string: tipo de proyecion de la capa.
Si queremos acceder a los datos de la tabla de atributos, podemos hacerlo llamado a data usando la @ y despues el nombre del campo que queremos recuperar, o bien escribiendo el nombre de la capa y acceder con el simbolo $ al nombre del campo directamente:
paises@data$name[1:10]## [1] "Indonesia" "Malaysia" "Chile" "Bolivia" "Peru" "Argentina"
## [7] "Dhekelia" "Cyprus" "India" "China"paises$name[1:10]## [1] "Indonesia" "Malaysia" "Chile" "Bolivia" "Peru" "Argentina"
## [7] "Dhekelia" "Cyprus" "India" "China"5.1 Añadir y eliminar atributos
Para añadir un campo nuevo escribimos el nombre que queramos despues de $ y le asignamos valores. Si queremos eliminarlo, le asignamos el valor null.
# añadimos una columna
paises$new <- sample(rep(letters, 2), length(paises), replace = TRUE)
paises$new[1:10]## [1] "b" "s" "h" "r" "k" "b" "q" "t" "u" "j"# eliminamos una columna
paises$new <- NULL6 Funciones de análisis espacial
6.1 Unir tablas (Join)
La funcion merge() nos permite unir una tabla (en formato data.frame) a la capa que tenemos cargada; simplemente, necesitamos tener un campo en comun entre ambas. Este campo no tiene por que llamarse igual, podemos definir los diferentes nombres dentro de la funcion.
dfr <- data.frame(ID = paises$name,
Campo_nuevo = round(runif(length(paises), 100, 1000)))
pm <- merge(paises, dfr, by.x = "name", by.y = "ID")
names(pm[,90:95])## [1] "name_ru" "name_sv" "name_tr" "name_vi" "name_zh"
## [6] "Campo_nuevo"6.2 Filtrado (Select by attributes) y Union (Append)
Se pueden separar los elementos de una capa para crear otras nuevas y tambien unir varias en una sola. Se tratan como dataframes. Si es una union, todas deben ser de tipo vector y tener la misma proyeccion. Así podemos utilizar funciones básicas como rbind().
separemos capas:
spain <- paises[paises$name == "Spain",]
france <- paises[paises$name == "France",]
portugal <- paises[paises$name == "Portugal",]
plot(spain)
plot(france)
plot(portugal)
Unamos los tres objetos:
x <- rbind(spain, france, portugal)
plot(x)
head(data.frame(x)[,1:3])## featurecla scalerank labelrank
## 67 Admin-0 country 0 2
## 22 Admin-0 country 0 2
## 139 Admin-0 country 0 26.3 Borrado
Si queremos eliminar una parte de la extension espacial de nuestra capa, podemos hacerlo con la funcion erase() del paquete raster.
1 Determinamos el area a recortar 2 con extent() definimos los limites externos del area. 3 obtenemos con lo anterior un objeto de la clase extent, que tenemos que convertir a un objeto espacial utilizando la funcion as() para convertirlo en un poligono (‘SpacialPolygons’)
library(raster)##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## select## The following object is masked from 'package:tidyr':
##
## extractb <- as(extent(0, 90, 0 ,90), "SpatialPolygons")
proj4string(b) <- proj4string(paises)## Warning in proj4string(paises): CRS object has comment, which is lost in outputplot(paises, lwd = 0.5)
plot(b, add= TRUE)
y usamos erase para hacer el recorte:
e <- erase(paises, b)## Loading required namespace: rgeosplot(e, lwd = 0.5)
6.4 Interseccion
Aca seleccionamos la parte comun de dos objetos espaciales
i <- raster::intersect(paises,b)## Warning in proj4string(x): CRS object has comment, which is lost in output## Warning in proj4string(y): CRS object has comment, which is lost in output## Warning in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func): spgeom1 and spgeom2
## have different proj4 strings## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
## unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4 stringsLa funcion intersect() es equivalente a la funcion crop():
i <- crop(paises,b)
plot(i,lwd = 0.5)
6.5 Diferencia
Podemos extraer el resultado de la diferencia entre dos capas. Y sera el espacio de la primera capa que no ocupe la segunda.
1 Creamos un objeto espacial de ejemplo con la misma extension que nuestro objeto espacial con las funciones extent() y as() igual que en el ejemplo de erase() y despues definimos la proyeccion de ese objeto (la misma que paises para que coincidan espacialmente)
b <- as(extent(-180, 180, -90, 90), "SpatialPolygons")
projection(b) <- projection(paises)
# comprobamos que coincidan espacialmente
plot(paises, lwd = 0.5)
plot(b,border = "red", add = TRUE)
Y ya podemos calcular la diferencia con la funcion symdif(), a la que pasamos los objetos espaciales que queremos restar:
# dif <- symdif(b,paises)
# plot(dif, col = "lightblue", lwd = 0.5)7 El paquete sf
Creemos un objeto espacial a partir de un conjunto de pares de coordenadas. Es tan fácil como crear o leer un data.frame con las coordenadas de las localizaciones a representar y los campos adicionales que sean necesarios.
datos <- data.frame(lon = c(-3.703889, -8.533333,
-5.9925, 2.166667,
-3.8, -1.644167, -0.9),
lat = c(40.4125, 42.883333,
37.3925, 41.4, 43.466667,
42.818333, 41.65),
nombre = c("Madrid", "Santiago",
"Sevilla", "Barcelona",
"Santander", "Pamplona", "Zaragoza"))Se convierte entonces en un objeto espacial con la funcion st_as_sf() identificando que campos contienen las coordenadas y que proyeccion queremos asignarle (definiendo con el codigo EPSG o el proj4string). Ademas de consultar el tipo de objeto con class(), sf permite consultar el tipo de geometría con st_geometry().
library(sf)## Warning: package 'sf' was built under R version 4.0.5## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1points_sf <- st_as_sf(datos, coords = c("lon", "lat"), crs = 4326)class(points_sf)## [1] "sf" "data.frame"st_geometry(points_sf)## Geometry set for 7 features
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -8.533333 ymin: 37.3925 xmax: 2.166667 ymax: 43.46667
## Geodetic CRS: WGS 84
## First 5 geometries:## POINT (-3.703889 40.4125)## POINT (-8.533333 42.88333)## POINT (-5.9925 37.3925)## POINT (2.166667 41.4)## POINT (-3.8 43.46667)8 Importación
La lectura de archivos, igual que con rgdal, es muy sencilla. Podemos importar un shp de ESRI o cualquier otro archivo de datos vectoriales con la funcion st_read(). El objeto cargado en memoria tendra una estructura de datos en formato tabla con una columna indicando el tipo de geometria espacial. Asi que es posible trabajar con el como si fuese un data.frame, pero con la ventaja de tener toda la informacion espacial necesaria.
paises <- st_read("./libro_sig_r/05_spatial_I/world_countries.shp")## Reading layer `world_countries' from data source `C:\Users\usuario\Desktop\ds\ds_rgee\estudio\libro_sig_r\05_spatial_I\world_countries.shp' using driver `ESRI Shapefile'
## Simple feature collection with 255 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84ciudades <- st_read("./libro_sig_r/05_spatial_I/world_Cities.shp")## Reading layer `World_Cities' from data source `C:\Users\usuario\Desktop\ds\ds_rgee\estudio\libro_sig_r\05_spatial_I\World_Cities.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2540 features and 13 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -176.1516 ymin: -54.792 xmax: 179.2219 ymax: 78.2
## Geodetic CRS: WGS 84print(paises)## Simple feature collection with 255 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84
## First 10 features:
## featurecla scalerank labelrank sovereignt sov_a3 adm0_dif level
## 1 Admin-0 country 5 2 Indonesia IDN 0 2
## 2 Admin-0 country 5 3 Malaysia MYS 0 2
## 3 Admin-0 country 6 2 Chile CHL 0 2
## 4 Admin-0 country 0 3 Bolivia BOL 0 2
## 5 Admin-0 country 0 2 Peru PER 0 2
## 6 Admin-0 country 0 2 Argentina ARG 0 2
## 7 Admin-0 country 3 3 United Kingdom GB1 1 2
## 8 Admin-0 country 6 5 Cyprus CYP 0 2
## 9 Admin-0 country 0 2 India IND 0 2
## 10 Admin-0 country 0 2 China CH1 1 2
## type admin adm0_a3 geou_dif
## 1 Sovereign country Indonesia IDN 0
## 2 Sovereign country Malaysia MYS 0
## 3 Sovereign country Chile CHL 0
## 4 Sovereign country Bolivia BOL 0
## 5 Sovereign country Peru PER 0
## 6 Sovereign country Argentina ARG 0
## 7 Dependency Dhekelia Sovereign Base Area ESB 0
## 8 Sovereign country Cyprus CYP 0
## 9 Sovereign country India IND 0
## 10 Country China CHN 0
## geounit gu_a3 su_dif subunit su_a3
## 1 Indonesia IDN 0 Indonesia IDN
## 2 Malaysia MYS 0 Malaysia MYS
## 3 Chile CHL 0 Chile CHL
## 4 Bolivia BOL 0 Bolivia BOL
## 5 Peru PER 0 Peru PER
## 6 Argentina ARG 0 Argentina ARG
## 7 Dhekelia Sovereign Base Area ESB 0 Dhekelia Sovereign Base Area ESB
## 8 Cyprus CYP 0 Cyprus CYP
## 9 India IND 0 India IND
## 10 China CHN 0 China CHN
## brk_diff name name_long brk_a3 brk_name brk_group abbrev postal
## 1 0 Indonesia Indonesia IDN Indonesia <NA> Indo. INDO
## 2 0 Malaysia Malaysia MYS Malaysia <NA> Malay. MY
## 3 0 Chile Chile CHL Chile <NA> Chile CL
## 4 0 Bolivia Bolivia BOL Bolivia <NA> Bolivia BO
## 5 0 Peru Peru PER Peru <NA> Peru PE
## 6 0 Argentina Argentina ARG Argentina <NA> Arg. AR
## 7 0 Dhekelia Dhekelia ESB Dhekelia <NA> Dhek. DH
## 8 0 Cyprus Cyprus CYP Cyprus <NA> Cyp. CY
## 9 0 India India IND India <NA> India IND
## 10 0 China China CHN China <NA> China CN
## formal_en formal_fr name_ciawf note_adm0 note_brk
## 1 Republic of Indonesia <NA> Indonesia <NA> <NA>
## 2 Malaysia <NA> Malaysia <NA> <NA>
## 3 Republic of Chile <NA> Chile <NA> <NA>
## 4 Plurinational State of Bolivia <NA> Bolivia <NA> <NA>
## 5 Republic of Peru <NA> Peru <NA> <NA>
## 6 Argentine Republic <NA> Argentina <NA> <NA>
## 7 <NA> <NA> <NA> U.K. Base <NA>
## 8 Republic of Cyprus <NA> Cyprus <NA> <NA>
## 9 Republic of India <NA> India <NA> <NA>
## 10 People's Republic of China <NA> China <NA> <NA>
## name_sort name_alt mapcolor7 mapcolor8 mapcolor9
## 1 Indonesia <NA> 6 6 6
## 2 Malaysia <NA> 2 4 3
## 3 Chile <NA> 5 1 5
## 4 Bolivia <NA> 1 5 2
## 5 Peru <NA> 4 4 4
## 6 Argentina <NA> 3 1 3
## 7 Dhekelia Sovereign Base Area <NA> 6 6 6
## 8 Cyprus <NA> 1 2 3
## 9 India <NA> 1 3 2
## 10 China <NA> 4 4 4
## mapcolor13 pop_est pop_rank gdp_md_est pop_year lastcensus gdp_year
## 1 11 260580739 17 3028000 2017 2010 2016
## 2 6 31381992 15 863000 2017 2010 2016
## 3 9 17789267 14 436100 2017 2002 2016
## 4 3 11138234 14 78350 2017 2001 2016
## 5 11 31036656 15 410400 2017 2007 2016
## 6 13 44293293 15 879400 2017 2010 2016
## 7 3 7850 5 314 2013 NA 2013
## 8 7 1221549 12 29260 2017 2001 2016
## 9 2 1281935911 18 8721000 2017 2011 2016
## 10 3 1379302771 18 21140000 2017 2010 2016
## economy income_grp wikipedia fips_10_ iso_a2
## 1 4. Emerging region: MIKT 4. Lower middle income NA ID ID
## 2 6. Developing region 3. Upper middle income NA MY MY
## 3 5. Emerging region: G20 3. Upper middle income NA CI CL
## 4 5. Emerging region: G20 4. Lower middle income NA BL BO
## 5 5. Emerging region: G20 3. Upper middle income NA PE PE
## 6 5. Emerging region: G20 3. Upper middle income NA AR AR
## 7 2. Developed region: nonG7 2. High income: nonOECD NA <NA> <NA>
## 8 6. Developing region 2. High income: nonOECD NA CY CY
## 9 3. Emerging region: BRIC 4. Lower middle income NA IN IN
## 10 3. Emerging region: BRIC 3. Upper middle income NA CH CN
## iso_a3 iso_a3_eh iso_n3 un_a3 wb_a2 wb_a3 woe_id woe_id_eh
## 1 IDN IDN 360 360 ID IDN 23424846 23424846
## 2 MYS MYS 458 458 MY MYS 23424901 23424901
## 3 CHL CHL 152 152 CL CHL 23424782 23424782
## 4 BOL BOL 068 068 BO BOL 23424762 23424762
## 5 PER PER 604 604 PE PER 23424919 23424919
## 6 ARG ARG 032 032 AR ARG 23424747 23424747
## 7 <NA> <NA> <NA> <NA> <NA> <NA> NA NA
## 8 CYP CYP 196 196 CY CYP -90 23424994
## 9 IND IND 356 356 IN IND 23424848 23424848
## 10 CHN CHN 156 156 CN CHN 23424781 23424781
## woe_note adm0_a3_is adm0_a3_us adm0_a3_un
## 1 Exact WOE match as country IDN IDN NA
## 2 Exact WOE match as country MYS MYS NA
## 3 Exact WOE match as country CHL CHL NA
## 4 Exact WOE match as country BOL BOL NA
## 5 Exact WOE match as country PER PER NA
## 6 Exact WOE match as country ARG ARG NA
## 7 No WOE equivalent. GBR ESB NA
## 8 WOE lists as subunit of united Cyprus CYP CYP NA
## 9 Exact WOE match as country IND IND NA
## 10 Exact WOE match as country CHN CHN NA
## adm0_a3_wb continent region_un subregion
## 1 NA Asia Asia South-Eastern Asia
## 2 NA Asia Asia South-Eastern Asia
## 3 NA South America Americas South America
## 4 NA South America Americas South America
## 5 NA South America Americas South America
## 6 NA South America Americas South America
## 7 NA Asia Asia Western Asia
## 8 NA Asia Asia Western Asia
## 9 NA Asia Asia Southern Asia
## 10 NA Asia Asia Eastern Asia
## region_wb name_len long_len abbrev_len tiny homepart
## 1 East Asia & Pacific 9 9 5 NA 1
## 2 East Asia & Pacific 8 8 6 NA 1
## 3 Latin America & Caribbean 5 5 5 NA 1
## 4 Latin America & Caribbean 7 7 7 NA 1
## 5 Latin America & Caribbean 4 4 4 NA 1
## 6 Latin America & Caribbean 9 9 4 NA 1
## 7 Europe & Central Asia 8 8 5 3 NA
## 8 Europe & Central Asia 6 6 4 NA 1
## 9 South Asia 5 5 5 NA 1
## 10 East Asia & Pacific 5 5 5 NA 1
## min_zoom min_label max_label ne_id wikidataid name_ar name_bn
## 1 0 1.7 6.7 1159320845 Q252 <NA> <NA>
## 2 0 3.0 8.0 1159321083 Q833 <NA> <NA>
## 3 0 1.7 6.7 1159320493 Q298 <NA> <NA>
## 4 0 3.0 7.5 1159320439 Q750 <NA> <NA>
## 5 0 2.0 7.0 1159321163 Q419 <NA> <NA>
## 6 0 2.0 7.0 1159320331 Q414 <NA> <NA>
## 7 0 6.5 11.0 1159320709 Q9206745 <NA> <NA>
## 8 0 4.5 9.5 1159320533 Q229 <NA> <NA>
## 9 0 1.7 6.7 1159320847 Q668 <NA> <NA>
## 10 0 1.7 5.7 1159320471 Q148 <NA> <NA>
## name_de name_en name_es
## 1 Indonesien Indonesia Indonesia
## 2 Malaysia Malaysia Malasia
## 3 Chile Chile Chile
## 4 Bolivien Bolivia Bolivia
## 5 Peru Peru Perú
## 6 Argentinien Argentina Argentina
## 7 Dekelia Dhekelia Cantonment Dekelia
## 8 Republik Zypern Cyprus Chipre
## 9 Indien India India
## 10 Volksrepublik China People's Republic of China República Popular China
## name_fr name_el name_hi name_hu
## 1 Indonésie <NA> <NA> Indonézia
## 2 Malaisie <NA> <NA> Malajzia
## 3 Chili <NA> <NA> Chile
## 4 Bolivie <NA> <NA> Bolívia
## 5 Pérou <NA> <NA> Peru
## 6 Argentine <NA> <NA> Argentína
## 7 Dhekelia <NA> <NA> Dekélia
## 8 Chypre <NA> <NA> Ciprus
## 9 Inde <NA> <NA> India
## 10 République populaire de Chine <NA> <NA> Kína
## name_id name_it name_ja name_ko
## 1 Indonesia Indonesia <NA> <NA>
## 2 Malaysia Malesia <NA> <NA>
## 3 Chili Cile <NA> <NA>
## 4 Bolivia Bolivia <NA> <NA>
## 5 Peru Perù <NA> <NA>
## 6 Argentina Argentina <NA> <NA>
## 7 Dhekelia Cantonment Base di Dhekelia <NA> <NA>
## 8 Siprus Cipro <NA> <NA>
## 9 India India <NA> <NA>
## 10 Republik Rakyat Tiongkok Cina <NA> <NA>
## name_nl name_pl name_pt name_ru name_sv
## 1 Indonesië Indonezja Indonésia <NA> Indonesien
## 2 Maleisië Malezja Malásia <NA> Malaysia
## 3 Chili Chile Chile <NA> Chile
## 4 Bolivia Boliwia Bolívia <NA> Bolivia
## 5 Peru Peru Peru <NA> Peru
## 6 Argentinië Argentyna Argentina <NA> Argentina
## 7 Dhekelia Cantonment Dhekelia Dekelia <NA> Dhekelia
## 8 Cyprus Cypr Chipre <NA> Cypern
## 9 India Indie Índia <NA> Indien
## 10 Volksrepubliek China Chinska Republika Ludowa China <NA> Kina
## name_tr name_vi name_zh geometry
## 1 Endonezya Indonesia <NA> MULTIPOLYGON (((117.7036 4....
## 2 Malezya Malaysia <NA> MULTIPOLYGON (((117.7036 4....
## 3 Sili Chile <NA> MULTIPOLYGON (((-69.51009 -...
## 4 Bolivya Bolivia <NA> MULTIPOLYGON (((-69.51009 -...
## 5 Peru Peru <NA> MULTIPOLYGON (((-69.51009 -...
## 6 Arjantin Argentina <NA> MULTIPOLYGON (((-67.28475 -...
## 7 Dhekelia Kantonu <NA> <NA> MULTIPOLYGON (((33.78094 34...
## 8 Kibris Cumhuriyeti <NA> <NA> MULTIPOLYGON (((33.78183 34...
## 9 Hindistan <NA> <NA> MULTIPOLYGON (((77.80035 35...
## 10 Çin Halk Cumhuriyeti <NA> <NA> MULTIPOLYGON (((78.91595 33...En el caso de una geodatabase es necesario indicar el nombre del layer como argumento. Por defecto, importa el primer layer y da aviso en caso de que exista más de uno. La función st_layers() nos permite conocer los layers de una geodatabase.
9 Visualizacion
Para representar espacialmente el objeto que hemos importado, podemos usar directamente la funcion plot(), que, al detectar que se trata de un objeto sf lo tratará de una maera diferente a la que conocemos de la clase sp; si el objeto tiene mas si el objeto tiene más de una variable (mas de un campo en la tabla de atributos) y no especificamos nada mas en la definicion de la funcion, plot() creara una figura con tantos mapas como variables (por defecto las 10 primeras), representando una a una y asignadole los intervalos de clase y color que considere mas adecuados en funcion del tipo de variable(discreta o continua) y de su rango.
plot(ciudades, max.plot =13)
Si queremos representar solo una variable del objeto, la seleccionamos mediante su nombre dentro de la funcion:
plot(ciudades["POP_RANK"])
10 Proyecciones y coordenadas
Las funciones fundamentales para proyectar, extraer las coordenadas o extraer el sistema de coordenadas son:
st_coordinates(): proyectar
st_crs(): extraer las coordenadas
st_transform():extraer el sistema de coordenadas
para asignar una proyecion es suficiente con indicar el codigo EPSG.
st_coordinates(ciudades) %>% head()## X Y
## 1 99.5290 16.472997
## 2 12.1900 -5.556996
## 3 100.3490 16.439004
## 4 15.7833 -17.049996
## 5 24.8300 -28.659996
## 6 105.4380 10.382004# extraer CRS (proyeccion)
st_crs(ciudades)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]# transformar proyeccion
# paises_rob <- st_transform(paises, 54030)
# plot(paises_rob["gdp_md_est"])11 Acceso, manipulacion y seleccion de variables
sf interpreta los objetos espaciales como un dataframe, Es por ello que el accesos a las variables se hace exactamente igual:
paises$name[1:10]## [1] "Indonesia" "Malaysia" "Chile" "Bolivia" "Peru" "Argentina"
## [7] "Dhekelia" "Cyprus" "India" "China"La manipulacion de los atributos de la clase sf es igual a la que conocimos con la coleccion de paquetes tidyverse.
En el siguiente ejemplo modificamos los atributos (variables) del objeto paises para crear una columna con el área, otra con la población y otra con la densidad de población:
library(tidyverse)
library(units)## udunits system database from C:/Users/usuario/Documents/R/win-library/4.0/units/share/udunitspaises2 <- mutate(paises,
area = st_area(paises) %>% set_units(km2),
pop_est = as.numeric(as.character(pop_est)),
densidad = pop_est/area)La funcion select() nos permite seleccionar y extraer las variables que queramos y, en combinacion con slice(), podemos extraer tambien filas (unidades espaciales) concretas.
dplyr::select(paises2, name, area:densidad)## Simple feature collection with 255 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84
## First 10 features:
## name area densidad
## 1 Indonesia 1879827.3749 [km2] 138.61950 [1/km2]
## 2 Malaysia 327884.8887 [km2] 95.71039 [1/km2]
## 3 Chile 736592.7499 [km2] 24.15075 [1/km2]
## 4 Bolivia 1086810.6457 [km2] 10.24855 [1/km2]
## 5 Peru 1289866.4011 [km2] 24.06192 [1/km2]
## 6 Argentina 2784305.9189 [km2] 15.90820 [1/km2]
## 7 Dhekelia 127.4228 [km2] 61.60593 [1/km2]
## 8 Cyprus 5395.0210 [km2] 226.42155 [1/km2]
## 9 India 3151452.4920 [km2] 406.77621 [1/km2]
## 10 China 9374591.2170 [km2] 147.13204 [1/km2]
## geometry
## 1 MULTIPOLYGON (((117.7036 4....
## 2 MULTIPOLYGON (((117.7036 4....
## 3 MULTIPOLYGON (((-69.51009 -...
## 4 MULTIPOLYGON (((-69.51009 -...
## 5 MULTIPOLYGON (((-69.51009 -...
## 6 MULTIPOLYGON (((-67.28475 -...
## 7 MULTIPOLYGON (((33.78094 34...
## 8 MULTIPOLYGON (((33.78183 34...
## 9 MULTIPOLYGON (((77.80035 35...
## 10 MULTIPOLYGON (((78.91595 33...Lo mismo pero seleccionandosolo el pais de la fila 5:
dplyr::select(paises2, name, area:densidad) %>% slice(5)## Simple feature collection with 1 feature and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.33756 ymin: -18.33775 xmax: -68.68425 ymax: -0.02909271
## Geodetic CRS: WGS 84
## name area densidad geometry
## 1 Peru 1289866 [km2] 24.06192 [1/km2] MULTIPOLYGON (((-69.51009 -...plot(paises2["area"])
cuando usamos select() para seleccionar columnas, la columna que contiene la geometria se mantiene sin necesidad de especificarla.
12 Funciones de análisis espacial
Para obtener una lista de todas las funciones disponibles del paquete sf podemos usar la siguiente funcion. Esta forma tambien es util para otros paquetes con clases propias.
methods(class = "sf")## [1] $<- [ [[<-
## [4] aggregate anti_join arrange
## [7] as.data.frame cbind coerce
## [10] dbDataType dbWriteTable distinct
## [13] dplyr_reconstruct extent extract
## [16] filter full_join gather
## [19] group_by group_split identify
## [22] initialize inner_join left_join
## [25] mask merge mutate
## [28] nest plot print
## [31] raster rasterize rbind
## [34] rename right_join rowwise
## [37] sample_frac sample_n select
## [40] semi_join separate separate_rows
## [43] show slice slotsFromS3
## [46] spread st_agr st_agr<-
## [49] st_area st_as_s2 st_as_sf
## [52] st_bbox st_boundary st_buffer
## [55] st_cast st_centroid st_collection_extract
## [58] st_convex_hull st_coordinates st_crop
## [61] st_crs st_crs<- st_difference
## [64] st_filter st_geometry st_geometry<-
## [67] st_inscribed_circle st_interpolate_aw st_intersection
## [70] st_intersects st_is st_is_valid
## [73] st_join st_line_merge st_m_range
## [76] st_make_valid st_nearest_points st_node
## [79] st_normalize st_point_on_surface st_polygonize
## [82] st_precision st_reverse st_sample
## [85] st_segmentize st_set_precision st_shift_longitude
## [88] st_simplify st_snap st_sym_difference
## [91] st_transform st_triangulate st_union
## [94] st_voronoi st_wrap_dateline st_write
## [97] st_z_range st_zm summarise
## [100] transform transmute ungroup
## [103] unite unnest
## see '?methods' for accessing help and source code13 Unir tablas (Join)
Unir un objeto sf con otra tabla con nuevos atributos es facil. Se usan las mismas funciones que aplicamos en el capitulo sobre tidyverse.
df <- data.frame(name = paises$name, var_dummy = rnorm(255))
head(df)## name var_dummy
## 1 Indonesia 0.1706780
## 2 Malaysia -1.2250713
## 3 Chile 0.3873329
## 4 Bolivia 0.4906678
## 5 Peru 2.2391261
## 6 Argentina 0.8922123paises_df <- left_join(paises, df)## Joining, by = "name"paises_df <- left_join(paises, df, by = "name")
names(paises_df)[90:96]## [1] "name_ru" "name_sv" "name_tr" "name_vi" "name_zh" "var_dummy"
## [7] "geometry"Si queremos hacer un spacial join, podemos usar la funcion st_join() o st_intersection()
14 Filtrado (Select by attributes)
La separacion de elementos se realiza con la funcion filter(), mediante la cual podemos hacer una (o varias) seleccion(es) por atributos.
paises_subset <- filter(paises2, area < set_units(1000, km2))
plot(paises_subset["area"])
En ocasiones nos puede interesar ordenar los atributos en funcion de una variable, para lo que utilizamos, en combinacion con select(), la funcion arrange().
dplyr::select(paises2, densidad, name) %>% arrange(desc(densidad))## Simple feature collection with 255 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84
## First 10 features:
## densidad name geometry
## 1 81941.480 [1/km2] Vatican MULTIPOLYGON (((12.45314 41...
## 2 20014.474 [1/km2] Macao MULTIPOLYGON (((113.5586 22...
## 3 11535.533 [1/km2] Singapore MULTIPOLYGON (((103.9608 1....
## 4 7961.280 [1/km2] Gibraltar MULTIPOLYGON (((-5.358387 3...
## 5 6940.234 [1/km2] Hong Kong MULTIPOLYGON (((114.2298 22...
## 6 3611.373 [1/km2] Maldives MULTIPOLYGON (((73.1631 -0....
## 7 2412.115 [1/km2] Bahrain MULTIPOLYGON (((50.55161 26...
## 8 1801.613 [1/km2] Sint Maarten MULTIPOLYGON (((-63.107 18....
## 9 1626.620 [1/km2] Monaco MULTIPOLYGON (((7.437454 43...
## 10 1278.334 [1/km2] Malta MULTIPOLYGON (((14.54802 35...dplyr::select(paises_subset, area, name) %>% arrange(desc(area))## Simple feature collection with 67 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -178.1857 ymin: -53.19256 xmax: 179.9067 ymax: 60.48078
## Geodetic CRS: WGS 84
## First 10 features:
## area name geometry
## 1 949.8635 [km2] Kiribati MULTIPOLYGON (((173.0383 1....
## 2 942.0536 [km2] Åland MULTIPOLYGON (((20.2776 60....
## 3 730.3126 [km2] Dominica MULTIPOLYGON (((-61.36286 1...
## 4 633.8454 [km2] Micronesia MULTIPOLYGON (((163.0261 5....
## 5 605.1307 [km2] Saint Lucia MULTIPOLYGON (((-60.88679 1...
## 6 603.1131 [km2] Tonga MULTIPOLYGON (((-173.9564 -...
## 7 584.9398 [km2] Bahrain MULTIPOLYGON (((50.55161 26...
## 8 579.3044 [km2] N. Mariana Is. MULTIPOLYGON (((145.2057 14...
## 9 574.8174 [km2] Isle of Man MULTIPOLYGON (((-4.612131 5...
## 10 564.2034 [km2] Guam MULTIPOLYGON (((144.8864 13...15 Union (Dissolve y Append)
En algunos objetos espaciales, por ejemplo los formados por multiples poligonos, se permiten operaciones adicionales sobre los atributos a través de la función group?by_(), equivalente a disolver poligonos segun una variable. En este ejemplo, las funciones st_area() y set_units() sirven para calcular el area en las unidades definidas, respectivamente.
# para establecer unidades
library(units)# anadimos un campo con el area de los paises
paises <- mutate(paises, area = st_area(paises) %>% set_units(km2))La disolucion (dissolve), a diferencia de la union (append), crea un nuevo y unico poligono a partir de todos los poligonos que lo componen. La funcion summarise() define que variables y con que datos agregadors se obtendran como resultado.
# disolver segun continentes
paises_group <- group_by(paises, continent) %>% summarise(area_continent = sum(area))## `summarise()` ungrouping output (override with `.groups` argument)## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planarhead(paises_group)## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84
## # A tibble: 6 x 3
## continent area_continent geometry
## <chr> [km2] <MULTIPOLYGON [°]>
## 1 Africa 29890688 (((-25.23131 16.923, -25.2692 16.923, -25.28083 16~
## 2 Antarctica 12358174 (((-51.73064 -82.06256, -51.2469 -82.01067, -50.76~
## 3 Asia 31182935 (((25.70436 40.1153, 25.66326 40.12617, 25.66326 4~
## 4 Europe 22979397 (((-180 68.55956, -180 68.98237, -179.7127 68.9179~
## 5 North Amer~ 24269146 (((-179.119 51.2154, -179.1367 51.22724, -179.1297~
## 6 Oceania 8513057 (((-180 -16.14911, -179.9875 -16.13812, -179.973 -~plot(paises_group["area_continent"])
Tambien se puede usar la combinacion entre la funcion group_by() y summarise() sin que esta ultima tenga argumentos, lo que provocaria la disolucion sin la creacion de nuevas variables.
Si queremos disolver todos los poligonos de nuestro objeto espacial, podemos usar la funcion st_union(). No obstante, esta disolucion provoca que perdamos el formato data.frame y pasemos a una clase sfc, que contiene unicamente la geometria.
disolv <- paises %>% st_union()## although coordinates are longitude/latitude, st_union assumes that they are planardisolv## Geometry set for 1 feature
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84## MULTIPOLYGON (((-179.9386 -16.12127, -179.9485 ...plot(disolv)
Patra volver a anadir variables tendriamos que aplicar la funcion st_sf().
# añadimos un data.frame
disolv_df <- disolv %>% st_sf()
disolv_df## Simple feature collection with 1 feature and 0 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84
## geometry
## 1 MULTIPOLYGON (((-179.9386 -...En el mismosentido, pero sin disolver limites internos, se puede aplicar la funcion st_combine().
disolv2 <- paises %>% st_combine()
class(disolv2)## [1] "sfc_MULTIPOLYGON" "sfc"plot(disolv2)
Para unir diferentes objetos espaciales manteniendo los poligonos originales (append) usamos simplemente la funcion rbind(). Si quisieramos unir diferentes geometrias o atributos, tendriamos que usar la funcion st_union(), lo que provocaria en caso de distintos tipos de geometrias una GEOMETRYCOLLECTION.
En el siguiente ejemplo generamos un objeto con todos los paises de Europa excepto Rusia y otro con todos los paises de Africa.
paises_Eu <- filter(paises2, continent == "Europe", name != "Russia")
plot(paises_Eu["area"])
paises_Afr <- filter(paises2, continent == "Africa")Con rbind(), simplemente los unimos en un objeto nuevo. Hay que tener en cuenta que deben tener la misma estructura y geometria.
EuroAsia <- rbind(paises_Eu, paises_Afr)
plot(EuroAsia["area"])
16 Interseccion (Intersect, Clip)
Algunas de las funciones aplicables con sp y rgdal tambien tambien estan disponibles en sf. Por ejemplo, podemos aplicar funciones de interseccion con st_intersection():
int <- st_intersection(ciudades, paises)## although coordinates are longitude/latitude, st_intersection assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all
## geometrieshead(int[,1:5])## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 100.352 ymin: -7.796997 xmax: 110.426 ymax: -0.02300395
## Geodetic CRS: WGS 84
## FID ObjectID CITY_NAME GMI_ADMIN ADMIN_NAME
## 162 2162 2458 Padang IDN-SBA Sumatera Barat
## 168 2168 2459 Bandung IDN-JBA Jawa Barat
## 173 2173 2460 Semarang IDN-JTN Jawa Tengah
## 179 2179 2461 Yogyakarta IDN-YGY Yogyakarta
## 238 2238 2467 Pontianak IDN-KBA Kalimantan Barat
## 240 2240 2468 Jambi IDN-JAM Jambi
## geometry
## 162 POINT (100.352 -0.954)
## 168 POINT (107.607 -6.912)
## 173 POINT (110.426 -6.971004)
## 179 POINT (110.369 -7.796997)
## 238 POINT (109.339 -0.02300395)
## 240 POINT (103.6167 -1.600002)En este ejemplo se han añadido al objeto ciudades los atributos del objeto paises (el nombre del pais) en aquellas localizaciones en las que había coincidencia espacial. Como funcionalidad extra, tambien es posible crear una sparse matrix, que no es mas que un listado indicando que indice de la primera capa coincide con que indice(s) de la segunda. En este ejemplo creamos una lista con los indices de los paises, conteniendo casa uno de ellos los indices de las ciudades con las que se intersecan.
int_spr <- st_intersects(paises, ciudades)## although coordinates are longitude/latitude, st_intersects assumes that they are planar# head(int_spr)O, por ejemplo, podemos extraer los indices que se intersecan con un solo polígono. Aqui extraemos todas las ciudades que intersectan con el poligono correspondiente a Indonesia (todas las que están dentro de su teritorio).
# primer elemento de la lista
head(paises$name_es)## [1] "Indonesia" "Malasia" "Chile" "Bolivia" "Perú" "Argentina"# ciudades de indonesia
dplyr::select(ciudades, CITY_NAME, CNTRY_NAME) %>% slice(int_spr[[1]])## Simple feature collection with 26 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 95.33334 ymin: -10.17247 xmax: 140.719 ymax: 5.566662
## Geodetic CRS: WGS 84
## First 10 features:
## CITY_NAME CNTRY_NAME geometry
## 1 Padang Indonesia POINT (100.352 -0.954)
## 2 Bandung Indonesia POINT (107.607 -6.912)
## 3 Semarang Indonesia POINT (110.426 -6.971004)
## 4 Yogyakarta Indonesia POINT (110.369 -7.796997)
## 5 Pontianak Indonesia POINT (109.339 -0.02300395)
## 6 Jambi Indonesia POINT (103.6167 -1.600002)
## 7 Palembang Indonesia POINT (104.757 -2.990997)
## 8 Tanjungkarang-Telukbetung Indonesia POINT (105.254 -5.409)
## 9 Jakarta Indonesia POINT (106.8 -6.166665)
## 10 Palu Indonesia POINT (119.851 -0.899001)17 Diferencia (Difference)
Otras funciones, como la de diferencia (st_difference()), son muy versatiles en sf. Por ejemplo, si queremos obtener la diferencia entre dos geometrias igual que haciamos con la funcion symdif() del paquete raster, podemos hacerlo de la siguiente manera: primero creamos un poligono con la dimension de, por ejemplo, Francia:
# creamos un poligono con la dimension de Francia
france <- filter(paises, name == "France")
poly_fr <- st_bbox(france) %>% st_as_sfc() %>% st_transform(4326)Y ya podemos calcular la diferencia entre las dos geometrias:
diff_pol <- st_difference(paises, poly_fr)## although coordinates are longitude/latitude, st_difference assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesdiff_pol[,1:5]## Simple feature collection with 150 features and 5 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS: WGS 84
## First 10 features:
## featurecla scalerank labelrank sovereignt sov_a3
## 1 Admin-0 country 5 2 Indonesia IDN
## 2 Admin-0 country 5 3 Malaysia MYS
## 3 Admin-0 country 6 2 Chile CHL
## 4 Admin-0 country 0 3 Bolivia BOL
## 5 Admin-0 country 0 2 Peru PER
## 6 Admin-0 country 0 2 Argentina ARG
## 9 Admin-0 country 0 2 India IND
## 10 Admin-0 country 0 2 China CH1
## 18 Admin-0 country 0 2 Pakistan PAK
## 26 Admin-0 country 0 2 South Korea KOR
## geometry
## 1 MULTIPOLYGON (((95.12696 5....
## 2 MULTIPOLYGON (((99.80044 6....
## 3 MULTIPOLYGON (((-109.22 -27...
## 4 POLYGON ((-69.51009 -17.505...
## 5 MULTIPOLYGON (((-80.84867 -...
## 6 MULTIPOLYGON (((-68.64234 -...
## 9 MULTIPOLYGON (((72.17262 10...
## 10 MULTIPOLYGON (((78.82426 33...
## 18 MULTIPOLYGON (((67.39039 24...
## 26 MULTIPOLYGON (((124.7245 37...plot(diff_pol, xlim = c(-20,10),
ylim = c(30,55))## Warning: plotting the first 9 out of 95 attributes; use max.plot = 95 to plot
## all
De manera inversa, se pede extraer la geometria de una capa a partir de la otra. Es lo mismo que la funcion intersect() o crop().
recort_fr <- st_crop(paises, poly_fr)## although coordinates are longitude/latitude, st_intersection assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesplot(recort_fr)## Warning: plotting the first 9 out of 95 attributes; use max.plot = 95 to plot
## all
18 Calculo de metricas espaciales
Una de las ventajas de trabajar con sf es que, en combinacion con el paquete units, se pueden obtener y manipular las metricas espaciales sin tener que hacer calculos adicionales. De manera muy intuitiva, la funcion st_area() calcula el area de un poligono o conjunto de poligonos de un objeto spacial en las unidades de la proyeccion, la cual podemos reconvertir con la funcion set_units():
library(units)
library(tidyverse)
st_area(paises) %>% head()## Units: [m^2]
## [1] 1.879827e+12 3.278849e+11 7.365927e+11 1.086811e+12 1.289866e+12
## [6] 2.784306e+12area <- st_area(paises) %>% set_units(km2)head(area)## Units: [km2]
## [1] 1879827.4 327884.9 736592.7 1086810.6 1289866.4 2784305.9La distancia entre puntos de un objeto espacial se calcula con st_distance(), dando como resultado un objeto de la clase units que contiene una matriz de distancias con el mismo número de filas y columnas que el número de puntos y una serie de atributos espaciales que identifican las unidades de medida.
# matrix de distancias
dist_ciudades <- st_distance(ciudades) %>% set_units(km)
dim(dist_ciudades)## [1] 2540 2540head(dist_ciudades,1)## Units: [km]
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 9907.961 87.64697 9903.808 9454.311 929.0536 4425.124 107.1802
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## [1,] 9172.47 10002.86 934.8077 993.4761 9766.836 262.9059 4192.065 937.1693
## [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 10068.86 309.9487 9347.197 3935.604 1003.045 10214.97 292.4192 9211.629
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## [1,] 978.9832 9969.104 309.3539 3777.282 9241.714 9931.194 1012.73 3678.925
## [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40]
## [1,] 319.4535 983.2194 10032.77 3570.849 9285.811 1032.829 351.0452 10059.17
## [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## [1,] 8087.867 1032.022 3483.045 9279.833 7950.376 1041.834 339.5038 3351.472
## [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## [1,] 9269.996 7798.772 334.7047 8355.575 1692.679 8491.791 3794.295 377.6704
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64]
## [1,] 8450.619 3756.077 8073.346 3998.342 1010.819 8528.004 9020.485 517.5213
## [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72]
## [1,] 3988.216 8595.021 269.9585 3605.598 1088.94 8825.836 8088.664 3440.755
## [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80]
## [1,] 299.3496 922.9937 8809.996 3317.568 1067.016 9149.458 331.3692 8511.162
## [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## [1,] 9099.661 3307.586 342.3692 1333.514 8647.706 4059.486 375.9524 9096.315
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96]
## [1,] 8960.294 1337.025 8997.484 9170.011 8921.622 5587.115 1096.213 1162.708
## [,97] [,98] [,99] [,100] [,101] [,102] [,103] [,104]
## [1,] 4134.8 9245.199 937.7043 8998.42 1281.965 4511.14 8894.961 1185.577
## [,105] [,106] [,107] [,108] [,109] [,110] [,111] [,112]
## [1,] 557.1164 3761.26 6788.928 9021.002 500.1536 1245.038 4093.285 742.6732
## [,113] [,114] [,115] [,116] [,117] [,118] [,119] [,120]
## [1,] 1181.839 3582.7 6343.148 3027.549 7662.491 3302.766 586.7669 4108.354
## [,121] [,122] [,123] [,124] [,125] [,126] [,127] [,128]
## [1,] 7688.265 536.2115 5382.852 2939.631 7615.832 701.7793 6375.319 2524.714
## [,129] [,130] [,131] [,132] [,133] [,134] [,135] [,136]
## [1,] 4073.937 7601.272 805.0967 6843.74 2246.005 7713.324 333.7403 6698.645
## [,137] [,138] [,139] [,140] [,141] [,142] [,143] [,144]
## [1,] 4689.992 2213.11 579.1532 8714.016 6908.526 544.6792 6314.613 2387.948
## [,145] [,146] [,147] [,148] [,149] [,150] [,151] [,152]
## [1,] 7888.79 460.6486 7970.924 8011.718 5929.882 2181.623 510.7153 6126.917
## [,153] [,154] [,155] [,156] [,157] [,158] [,159] [,160]
## [1,] 7998.402 8043.919 6517.023 2858.058 643.9625 7096.459 7898.002 4796.16
## [,161] [,162] [,163] [,164] [,165] [,166] [,167] [,168]
## [1,] 568.2698 1929.593 1638.809 7979.654 1614.32 657.9922 5207.165 2735.087
## [,169] [,170] [,171] [,172] [,173] [,174] [,175] [,176]
## [1,] 7952.813 353.2842 1341.714 6774.732 2857.18 8007.134 425.0513 7789.822
## [,177] [,178] [,179] [,180] [,181] [,182] [,183] [,184]
## [1,] 7923.73 404.0086 2937.808 7987.459 443.5001 1005.408 8018.33 7862.2
## [,185] [,186] [,187] [,188] [,189] [,190] [,191] [,192]
## [1,] 278.7714 7959.636 417.3243 7928.59 1157.767 5803.327 6462.62 7786.657
## [,193] [,194] [,195] [,196] [,197] [,198] [,199] [,200]
## [1,] 8202.909 903.6056 322.095 460.1847 8016.729 5635.125 984.3662 7924.899
## [,201] [,202] [,203] [,204] [,205] [,206] [,207] [,208]
## [1,] 4925.642 7970.987 654.1083 4783.783 812.5639 7973.248 766.0723 4694.901
## [,209] [,210] [,211] [,212] [,213] [,214] [,215] [,216]
## [1,] 5646.979 7965.224 711.7383 8037.205 828.44 4941.678 3295.157 7881.742
## [,217] [,218] [,219] [,220] [,221] [,222] [,223] [,224]
## [1,] 7894.827 724.0071 3180.444 7930.984 4865.65 809.9774 7852.155 3004.251
## [,225] [,226] [,227] [,228] [,229] [,230] [,231] [,232]
## [1,] 7951.796 4513.556 795.572 6522.699 3051.012 834.5007 4500.567 774.0058
## [,233] [,234] [,235] [,236] [,237] [,238] [,239] [,240]
## [1,] 7178.598 2306.785 4489.044 7249.54 914.2495 2118.635 515.3021 2048.772
## [,241] [,242] [,243] [,244] [,245] [,246] [,247] [,248]
## [1,] 4246.214 6965.177 2472.829 365.6481 2228.208 4140.707 6999.947 2362.107
## [,249] [,250] [,251] [,252] [,253] [,254] [,255] [,256]
## [1,] 2262.323 375.4066 4082.623 7063.057 376.2643 4106.146 260.4949 2500.855
## [,257] [,258] [,259] [,260] [,261] [,262] [,263] [,264]
## [1,] 2282.619 7015.68 362.4048 4242.103 842.5981 7208.248 2628.819 2279.151
## [,265] [,266] [,267] [,268] [,269] [,270] [,271] [,272]
## [1,] 840.8221 2597.941 7172.478 807.2181 8302.222 2474.213 2945.26 871.9124
## [,273] [,274] [,275] [,276] [,277] [,278] [,279] [,280]
## [1,] 819.7731 7079.028 2760.461 8605.951 3864.845 8504.564 2165.118 833.147
## [,281] [,282] [,283] [,284] [,285] [,286] [,287] [,288]
## [1,] 7086.985 3397.603 779.3029 799.3014 2596.128 6934.648 7096.579 3268.563
## [,289] [,290] [,291] [,292] [,293] [,294] [,295] [,296]
## [1,] 757.1469 7454.531 964.9735 7264.254 7150.611 2153.239 3988.313 880.0873
## [,297] [,298] [,299] [,300] [,301] [,302] [,303] [,304]
## [1,] 1890.425 871.2084 7345.726 7118.174 4698.406 897.0077 1616.685 8412.41
## [,305] [,306] [,307] [,308] [,309] [,310] [,311] [,312]
## [1,] 5137.581 7281.866 812.2109 2696.457 7522.839 3616.738 7646.649 8231.136
## [,313] [,314] [,315] [,316] [,317] [,318] [,319] [,320]
## [1,] 860.1769 2729.644 3539.521 7335.168 8737.656 929.5727 2608.467 7201.743
## [,321] [,322] [,323] [,324] [,325] [,326] [,327] [,328]
## [1,] 3067.429 1005.621 7113.383 9966.553 3138.355 2746.115 940.59 3041.822
## [,329] [,330] [,331] [,332] [,333] [,334] [,335] [,336]
## [1,] 953.5826 7064.614 10119.71 1085.059 3011.803 9918.9 865.955 2810.015
## [,337] [,338] [,339] [,340] [,341] [,342] [,343] [,344]
## [1,] 7159.136 1118.233 2695.003 7483.515 9946.894 3321.192 7505.625 1009.986
## [,345] [,346] [,347] [,348] [,349] [,350] [,351] [,352]
## [1,] 2821.386 7844.074 9302.227 3277.437 939.0106 2466.179 9629.518 1143.307
## [,353] [,354] [,355] [,356] [,357] [,358] [,359] [,360]
## [1,] 7800.643 3962.94 1084.345 4259.107 9615.596 7505.292 4992.655 7975.507
## [,361] [,362] [,363] [,364] [,365] [,366] [,367] [,368]
## [1,] 1418.872 4026.922 7487.007 8275.889 7934.213 5058.225 1462.22 3524.65
## [,369] [,370] [,371] [,372] [,373] [,374] [,375] [,376]
## [1,] 9783.731 1112.143 7846.841 9834.011 3472.409 1294.631 5331.99 9483.058
## [,377] [,378] [,379] [,380] [,381] [,382] [,383] [,384]
## [1,] 8133.809 1386.806 9956.323 3588.886 8752.887 9771.964 5628.583 1492.841
## [,385] [,386] [,387] [,388] [,389] [,390] [,391] [,392]
## [1,] 10213.67 3417.368 8632.389 5437.116 3589.074 10218.74 1499.29 8910.645
## [,393] [,394] [,395] [,396] [,397] [,398] [,399] [,400]
## [1,] 3697.834 2032.303 8773.728 10281.56 10118.64 3570.613 1545.45 10224.05
## [,401] [,402] [,403] [,404] [,405] [,406] [,407] [,408]
## [1,] 8711.862 9985.835 3536.338 9959.864 1606.717 8764.748 9985.374 10025.62
## [,409] [,410] [,411] [,412] [,413] [,414] [,415] [,416]
## [1,] 3502.486 9945.218 1719.849 8785.607 10047.63 7085.02 3559.033 1773.408
## [,417] [,418] [,419] [,420] [,421] [,422] [,423] [,424]
## [1,] 9613.453 7545.202 9574.426 3637.236 10086.22 1746.612 7793.474 8754.128
## [,425] [,426] [,427] [,428] [,429] [,430] [,431] [,432]
## [1,] 9996.662 9155.981 1079.36 3533.833 9990.144 9216.57 7800.835 1114.312
## [,433] [,434] [,435] [,436] [,437] [,438] [,439] [,440]
## [1,] 9284.288 1138.999 3547.519 9185.737 7014.903 9256.674 9538.741 1182.975
## [,441] [,442] [,443] [,444] [,445] [,446] [,447] [,448]
## [1,] 9235.572 9596.597 3559.493 903.1049 7756.633 831.6011 9059.568 8838.202
## [,449] [,450] [,451] [,452] [,453] [,454] [,455] [,456]
## [1,] 3890.145 7669.282 981.804 3962.334 9369.824 997.0791 3652.635 9245.08
## [,457] [,458] [,459] [,460] [,461] [,462] [,463] [,464]
## [1,] 8182.795 9381.644 997.4258 3892.106 9099.971 5505.566 3758.921 821.4296
## [,465] [,466] [,467] [,468] [,469] [,470] [,471] [,472]
## [1,] 9178.238 9263.82 3553.646 3753.658 5586.402 9538.415 3581.591 5631.762
## [,473] [,474] [,475] [,476] [,477] [,478] [,479] [,480]
## [1,] 9780.235 3662.673 5464.067 9711.303 3624.536 9664.582 5769.978 9847.998
## [,481] [,482] [,483] [,484] [,485] [,486] [,487] [,488]
## [1,] 3298.161 9900.564 5583.343 9739.55 3542.589 6265.762 9597.773 3105.283
## [,489] [,490] [,491] [,492] [,493] [,494] [,495] [,496]
## [1,] 6515.546 9912.822 1675.384 9350.647 5631.894 1666.194 6016.437 9735.066
## [,497] [,498] [,499] [,500] [,501] [,502] [,503] [,504]
## [1,] 1602.274 9633.679 6226.11 9539.12 2314.357 9724.904 6080.322 9718.924
## [,505] [,506] [,507] [,508] [,509] [,510] [,511] [,512]
## [1,] 2307.779 6660.079 9813.703 2146.95 5821.849 9893.036 6047.474 2114.341
## [,513] [,514] [,515] [,516] [,517] [,518] [,519] [,520]
## [1,] 6341.112 6430.919 6173.263 5979.179 3888.226 9197.864 3553.483 7243.967
## [,521] [,522] [,523] [,524] [,525] [,526] [,527] [,528]
## [1,] 3530.534 8436.244 3571.932 9454.588 5315.302 8332.474 3324.522 7625.471
## [,529] [,530] [,531] [,532] [,533] [,534] [,535] [,536]
## [1,] 9753.378 9839.388 9986.653 9917.218 10176.54 9873.814 9903.234 9799.942
## [,537] [,538] [,539] [,540] [,541] [,542] [,543] [,544]
## [1,] 9953.794 9903.887 9782.852 7511.739 17425.93 13935.1 14142.78 16553.51
## [,545] [,546] [,547] [,548] [,549] [,550] [,551] [,552]
## [1,] 16706.9 13933.8 17270.05 14284.39 17355 14416.46 17174.53 14425.86
## [,553] [,554] [,555] [,556] [,557] [,558] [,559] [,560]
## [1,] 17167.6 13814.83 17219.17 13827.8 13854.13 17223.42 13864.7 18743.28
## [,561] [,562] [,563] [,564] [,565] [,566] [,567] [,568]
## [1,] 14014.79 19062.68 14082.51 19157.24 14216.42 18916.51 13643.78 19079.22
## [,569] [,570] [,571] [,572] [,573] [,574] [,575] [,576]
## [1,] 13744.69 13389.54 18719.06 13421.44 9838.224 13639.61 10313.55 13637.8
## [,577] [,578] [,579] [,580] [,581] [,582] [,583] [,584]
## [1,] 13712.22 10434.87 13700.02 12718.66 13736.91 10656.78 13185.81 10186.07
## [,585] [,586] [,587] [,588] [,589] [,590] [,591] [,592]
## [1,] 13231.66 11745.8 13334.1 14804.85 13418.64 18282.6 16716.41 13484.5
## [,593] [,594] [,595] [,596] [,597] [,598] [,599] [,600]
## [1,] 17503.16 16534.96 13507.97 15980.1 17642.94 16674.13 15617.33 17449.34
## [,601] [,602] [,603] [,604] [,605] [,606] [,607] [,608]
## [1,] 16696.27 15849.12 17463.11 16594.88 15626.01 17136.59 16755.22 17789.83
## [,609] [,610] [,611] [,612] [,613] [,614] [,615] [,616]
## [1,] 17409.6 16716.53 17389.5 17733.08 16958.96 17318.43 17635.55 17126.78
## [,617] [,618] [,619] [,620] [,621] [,622] [,623] [,624]
## [1,] 16985.63 17422.27 16971.26 16930.06 17249.69 16582.12 17006.58 17141.81
## [,625] [,626] [,627] [,628] [,629] [,630] [,631] [,632]
## [1,] 16464.15 16886.21 18704.34 17267.34 16359.29 18420.57 16974.78 16535.92
## [,633] [,634] [,635] [,636] [,637] [,638] [,639] [,640]
## [1,] 16544.95 18251.92 16831.29 17952.9 16537.96 17848.11 16697.2 16464.04
## [,641] [,642] [,643] [,644] [,645] [,646] [,647] [,648]
## [1,] 17968.49 17041.81 16439.25 17883.92 18228.28 16389.87 18251.16 18207.73
## [,649] [,650] [,651] [,652] [,653] [,654] [,655] [,656]
## [1,] 15778.3 18298.71 18073.98 15841.09 18351.97 15797.97 17936.81 18377.28
## [,657] [,658] [,659] [,660] [,661] [,662] [,663] [,664]
## [1,] 18010.99 15725.43 18424.39 18019.99 15903.84 17688 18476.04 17670.4
## [,665] [,666] [,667] [,668] [,669] [,670] [,671] [,672]
## [1,] 15901.8 18233.75 17414.2 15799.59 17343.77 18480.74 17406.95 15770.38
## [,673] [,674] [,675] [,676] [,677] [,678] [,679] [,680]
## [1,] 18631.48 17334.7 15832.02 18019.27 17482.54 16018.68 17465.16 18494.3
## [,681] [,682] [,683] [,684] [,685] [,686] [,687] [,688]
## [1,] 15871.62 18542.04 17501.36 15920.37 17484.81 18576.48 15967.86 17427.2
## [,689] [,690] [,691] [,692] [,693] [,694] [,695] [,696]
## [1,] 18616.52 17312.07 16019.24 18622.38 16112.22 18756.18 17367.15 16019.62
## [,697] [,698] [,699] [,700] [,701] [,702] [,703] [,704]
## [1,] 17206.99 18927.7 17160.92 16243.89 19066.66 17128.2 16320.71 18789.08
## [,705] [,706] [,707] [,708] [,709] [,710] [,711] [,712]
## [1,] 17097.27 16972.18 16287.17 18838.22 16994.92 16301.44 18953.31 16328.98
## [,713] [,714] [,715] [,716] [,717] [,718] [,719] [,720]
## [1,] 16413.16 16234.89 18911.27 16340.68 16376.59 19143.13 16125.51 16648.25
## [,721] [,722] [,723] [,724] [,725] [,726] [,727] [,728]
## [1,] 19177.32 16318.79 16530.91 19154.81 9552.796 16690.25 19228.54 12493.11
## [,729] [,730] [,731] [,732] [,733] [,734] [,735] [,736]
## [1,] 17031.92 18195.1 12394.67 17027.32 18273.31 12936.31 16842.59 18267
## [,737] [,738] [,739] [,740] [,741] [,742] [,743] [,744]
## [1,] 12590.63 18307.33 16812.63 13102.69 16908.01 18324.54 9345.632 16814.96
## [,745] [,746] [,747] [,748] [,749] [,750] [,751] [,752]
## [1,] 13138.25 18354.99 16737.51 16761.22 9489.227 18417.41 10261.71 16853.66
## [,753] [,754] [,755] [,756] [,757] [,758] [,759] [,760]
## [1,] 9262.611 18394.95 19408.01 15504.56 9562.766 11050.09 19273.96 12091.03
## [,761] [,762] [,763] [,764] [,765] [,766] [,767] [,768]
## [1,] 9737.507 13015.48 19306.92 16154.42 19262.37 9319.785 10717.96 15789.47
## [,769] [,770] [,771] [,772] [,773] [,774] [,775] [,776]
## [1,] 19454.4 9765.261 12482.75 16098.35 16038.7 12845.15 9245.397 19014.7
## [,777] [,778] [,779] [,780] [,781] [,782] [,783] [,784]
## [1,] 18933.71 15782.54 12954.19 9379.742 18863.33 15246.53 12001.89 9200.955
## [,785] [,786] [,787] [,788] [,789] [,790] [,791] [,792]
## [1,] 19404.33 15735.8 13015.29 9278.114 16151.74 14964.12 12036.78 9340.353
## [,793] [,794] [,795] [,796] [,797] [,798] [,799] [,800]
## [1,] 16800.77 14970 9408.118 9420.442 17344.79 15003.65 9653.606 9249.313
## [,801] [,802] [,803] [,804] [,805] [,806] [,807] [,808]
## [1,] 18039.8 15136.38 9309.009 9408.5 15307.38 16491.82 9445.078 10225.57
## [,809] [,810] [,811] [,812] [,813] [,814] [,815] [,816]
## [1,] 9355.52 14711.73 17913.04 9476.933 11725.01 10767.57 17615.35 9457.272
## [,817] [,818] [,819] [,820] [,821] [,822] [,823] [,824]
## [1,] 17856.71 10374.55 9335.966 9958.899 18356.24 14506.97 9467.272 10181.66
## [,825] [,826] [,827] [,828] [,829] [,830] [,831] [,832]
## [1,] 18507.77 14687 10258.77 18191.14 9492.236 14764.97 18603.81 10263.38
## [,833] [,834] [,835] [,836] [,837] [,838] [,839] [,840]
## [1,] 9241.24 14644.44 10292.46 18395.61 9225.014 10855.92 10345.72 18437.72
## [,841] [,842] [,843] [,844] [,845] [,846] [,847] [,848]
## [1,] 9783.441 18053.39 9389.487 12872.08 9347.083 11512.29 8610.226 11455.7
## [,849] [,850] [,851] [,852] [,853] [,854] [,855] [,856]
## [1,] 10176.78 13077.1 9645.12 9544.17 12973.25 13083.01 9683.507 10608.59
## [,857] [,858] [,859] [,860] [,861] [,862] [,863] [,864]
## [1,] 10944.18 13434.92 11342.4 11995.04 9469.233 13446.89 11313.96 12112.05
## [,865] [,866] [,867] [,868] [,869] [,870] [,871] [,872]
## [1,] 11630.35 13658.82 11472.74 12049.81 11845.82 12123.91 10094.78 13980.3
## [,873] [,874] [,875] [,876] [,877] [,878] [,879] [,880]
## [1,] 11610.79 12176.12 12207.86 10234.03 12185.16 10192.65 11659.24 13313.16
## [,881] [,882] [,883] [,884] [,885] [,886] [,887] [,888]
## [1,] 12213.66 10234.51 11782.55 13422.74 10189.43 11798.59 10293.79 12264.82
## [,889] [,890] [,891] [,892] [,893] [,894] [,895] [,896]
## [1,] 12425.12 11943.4 12584.13 12264.45 10239.37 9687.494 11964.7 12177.33
## [,897] [,898] [,899] [,900] [,901] [,902] [,903] [,904]
## [1,] 10439.23 10257.32 12709.01 12621.36 10359.33 12278.57 12272.4 11574.61
## [,905] [,906] [,907] [,908] [,909] [,910] [,911] [,912]
## [1,] 12574.17 12273.97 10187.48 12235.08 11998.13 12578.67 10417.67 12234.43
## [,913] [,914] [,915] [,916] [,917] [,918] [,919] [,920]
## [1,] 10328.16 12225.44 12650.31 10407.62 12470.89 13140.91 12078.97 12340.11
## [,921] [,922] [,923] [,924] [,925] [,926] [,927] [,928]
## [1,] 10349.11 13153.27 12034.96 10223.53 16173.78 12124.59 10398.76 12896.54
## [,929] [,930] [,931] [,932] [,933] [,934] [,935] [,936]
## [1,] 16176.45 12170.34 9508.141 13297.23 15898.38 12079.61 9801.138 13691
## [,937] [,938] [,939] [,940] [,941] [,942] [,943] [,944]
## [1,] 12130.21 13817.98 15907.53 9662.479 12197.73 13790.08 15941.62 9736.764
## [,945] [,946] [,947] [,948] [,949] [,950] [,951] [,952]
## [1,] 11960.82 13642.38 15958.54 9945.036 11976.8 14207.34 15979.75 9913.377
## [,953] [,954] [,955] [,956] [,957] [,958] [,959] [,960]
## [1,] 14198.48 11892.61 16033.19 9963.566 11998.42 13540.15 11854.6 9745.626
## [,961] [,962] [,963] [,964] [,965] [,966] [,967] [,968]
## [1,] 11923.27 16024.56 9739.566 11436.1 13993.41 10056.55 11203.88 16069.28
## [,969] [,970] [,971] [,972] [,973] [,974] [,975] [,976]
## [1,] 9801.109 13999.14 11225.86 10224.75 11432.04 16070.54 11601.41 14310.19
## [,977] [,978] [,979] [,980] [,981] [,982] [,983] [,984]
## [1,] 11245.73 10281.08 15890.66 14592.21 11764.06 10334.97 10772.98 10232.04
## [,985] [,986] [,987] [,988] [,989] [,990] [,991] [,992]
## [1,] 14472.76 15902.87 14816.95 10491.98 9831.631 10975.95 15890.81 9862.087
## [,993] [,994] [,995] [,996] [,997] [,998] [,999] [,1000]
## [1,] 10574.61 14788.26 15891.9 9931.776 10729.52 15048.99 15890.34 9990.564
## [,1001] [,1002] [,1003] [,1004] [,1005] [,1006] [,1007] [,1008]
## [1,] 10914.24 14291.81 15908.33 10525.11 11007.57 10644.54 10976.9 14513.1
## [,1009] [,1010] [,1011] [,1012] [,1013] [,1014] [,1015] [,1016]
## [1,] 15937.03 10135.15 10838.42 10829.31 10098.94 14571.03 10835.59 15974.52
## [,1017] [,1018] [,1019] [,1020] [,1021] [,1022] [,1023] [,1024]
## [1,] 10829.29 11614.51 14874.39 16014.82 10749.21 11840.55 14869.02 10791.47
## [,1025] [,1026] [,1027] [,1028] [,1029] [,1030] [,1031] [,1032]
## [1,] 11752.45 15998.25 15119.64 11707.74 10878.74 16005.06 15146.35 12218.3
## [,1033] [,1034] [,1035] [,1036] [,1037] [,1038] [,1039] [,1040]
## [1,] 16034.99 12141.8 10871.15 15037.46 16004.83 15283.02 11997.17 10771.49
## [,1041] [,1042] [,1043] [,1044] [,1045] [,1046] [,1047] [,1048]
## [1,] 15203.26 15874.22 12005.15 10568.13 15312.46 12015.82 12089 10710.12
## [,1049] [,1050] [,1051] [,1052] [,1053] [,1054] [,1055] [,1056]
## [1,] 15894.33 15326.37 15840.05 12146.87 10616.52 15875.54 12031.99 15639.61
## [,1057] [,1058] [,1059] [,1060] [,1061] [,1062] [,1063] [,1064]
## [1,] 10671.8 12134.63 15600.3 15917.02 12080.64 10690.57 10666.52 11993.95
## [,1065] [,1066] [,1067] [,1068] [,1069] [,1070] [,1071] [,1072]
## [1,] 9908.853 15628.18 15932.53 10680.81 12132.41 10031.2 11943.43 10757.77
## [,1073] [,1074] [,1075] [,1076] [,1077] [,1078] [,1079] [,1080]
## [1,] 9974.645 11947.31 11865.56 11940.17 15959.65 15969.74 10069 11825.19
## [,1081] [,1082] [,1083] [,1084] [,1085] [,1086] [,1087] [,1088]
## [1,] 15894.9 10032.06 11964.5 15983.93 11516.12 10233.37 16080.89 11603.75
## [,1089] [,1090] [,1091] [,1092] [,1093] [,1094] [,1095] [,1096]
## [1,] 15920.35 10119.95 11741.11 15913.08 11766.35 10084.48 11581.45 9650.228
## [,1097] [,1098] [,1099] [,1100] [,1101] [,1102] [,1103] [,1104]
## [1,] 15708.78 15917.65 9995.572 11534.2 11669.32 11707.97 11841.12 10365.16
## [,1105] [,1106] [,1107] [,1108] [,1109] [,1110] [,1111] [,1112]
## [1,] 11352.39 15188.38 15954.44 12237.19 11395.97 10057.03 15370.11 10268.87
## [,1113] [,1114] [,1115] [,1116] [,1117] [,1118] [,1119] [,1120]
## [1,] 15946.36 11509.07 15495.74 12294.39 11311.24 10185.5 10447.2 17758.46
## [,1121] [,1122] [,1123] [,1124] [,1125] [,1126] [,1127] [,1128]
## [1,] 15413.68 11774.79 11943.81 10596.59 17872.79 15545.55 11702.71 11755.77
## [,1129] [,1130] [,1131] [,1132] [,1133] [,1134] [,1135] [,1136]
## [1,] 8479.803 15532.54 17936.47 8976.581 11702.98 9057.951 15590.21 9281.776
## [,1137] [,1138] [,1139] [,1140] [,1141] [,1142] [,1143] [,1144]
## [1,] 18017.39 11799.77 18010.9 9042.075 15708.22 11907.47 18068.3 9267.962
## [,1145] [,1146] [,1147] [,1148] [,1149] [,1150] [,1151] [,1152]
## [1,] 15740.52 11937.07 18124.75 16309.96 9396.437 17448.86 11899.19 9091.849
## [,1153] [,1154] [,1155] [,1156] [,1157] [,1158] [,1159] [,1160]
## [1,] 16289.48 11962.5 17561.63 9186.993 11824.9 8917.799 17595.84 11704.06
## [,1161] [,1162] [,1163] [,1164] [,1165] [,1166] [,1167] [,1168]
## [1,] 16330.8 8639.758 11654.3 17623.8 11815.17 8993.688 16314.18 11666.74
## [,1169] [,1170] [,1171] [,1172] [,1173] [,1174] [,1175] [,1176]
## [1,] 17622.25 9062.963 16335.45 17605.48 9084.301 16338.84 17785.62 11783.15
## [,1177] [,1178] [,1179] [,1180] [,1181] [,1182] [,1183] [,1184]
## [1,] 16356.34 9431.698 11651.45 11811.12 17568.9 16363.31 9217.728 11741.19
## [,1185] [,1186] [,1187] [,1188] [,1189] [,1190] [,1191] [,1192]
## [1,] 11565.5 16370.03 11578.47 16905.14 9685.056 16350.7 11582.7 16386.86
## [,1193] [,1194] [,1195] [,1196] [,1197] [,1198] [,1199] [,1200]
## [1,] 9212.83 11651.65 16380.47 16922.3 11590.49 16982.04 11584.96 16394.99
## [,1201] [,1202] [,1203] [,1204] [,1205] [,1206] [,1207] [,1208]
## [1,] 9757.903 11607.25 17099.56 16393.42 16419.45 11692.61 9776.45 17169.17
## [,1209] [,1210] [,1211] [,1212] [,1213] [,1214] [,1215] [,1216]
## [1,] 16438.17 11417.99 9270.479 17241.58 11471.19 15687.03 9748.618 11425.75
## [,1217] [,1218] [,1219] [,1220] [,1221] [,1222] [,1223] [,1224]
## [1,] 17436.4 15701.96 11473.2 17517.79 9644.488 15731.53 11410.39 16853.97
## [,1225] [,1226] [,1227] [,1228] [,1229] [,1230] [,1231] [,1232]
## [1,] 9182.206 15753.78 11486.85 11440.29 16942.99 9081.453 15753.4 11515.83
## [,1233] [,1234] [,1235] [,1236] [,1237] [,1238] [,1239] [,1240]
## [1,] 11452.69 9008.335 17050.69 15761.72 11443.27 16045.83 11542.23 9319.22
## [,1241] [,1242] [,1243] [,1244] [,1245] [,1246] [,1247] [,1248]
## [1,] 17173.08 13022.41 11108.43 9419.329 11179.41 9545.972 17176.98 14357.19
## [,1249] [,1250] [,1251] [,1252] [,1253] [,1254] [,1255] [,1256]
## [1,] 11008.69 11176.54 8972.259 17275.43 11023.35 17045.97 17192.64 17050.11
## [,1257] [,1258] [,1259] [,1260] [,1261] [,1262] [,1263] [,1264]
## [1,] 9055.867 11259.92 11041.72 17184.78 17063.95 16775.08 9358.626 11111.29
## [,1265] [,1266] [,1267] [,1268] [,1269] [,1270] [,1271] [,1272]
## [1,] 17057.45 8540.488 16732.33 11138.28 17067.76 16757.99 8793.852 11367.32
## [,1273] [,1274] [,1275] [,1276] [,1277] [,1278] [,1279] [,1280]
## [1,] 17143.2 16827.8 8575.505 11312.8 16836.93 8774.514 17241.7 11348.96
## [,1281] [,1282] [,1283] [,1284] [,1285] [,1286] [,1287] [,1288]
## [1,] 16949.5 11397.61 8486.768 16600.66 11351.75 17010.93 8971.913 17436.31
## [,1289] [,1290] [,1291] [,1292] [,1293] [,1294] [,1295] [,1296]
## [1,] 11409.64 8818.224 17368.46 16600.14 11222.78 8650.894 17729.5 16796.41
## [,1297] [,1298] [,1299] [,1300] [,1301] [,1302] [,1303] [,1304]
## [1,] 11302.13 11310.26 9064.518 16703.98 17749.6 11367.89 11258.32 8941.846
## [,1305] [,1306] [,1307] [,1308] [,1309] [,1310] [,1311] [,1312]
## [1,] 17285.59 11346.15 13338.28 11220.33 13428.4 8947.677 11182.36 17299.5
## [,1313] [,1314] [,1315] [,1316] [,1317] [,1318] [,1319] [,1320]
## [1,] 9442.691 17323.09 13481.06 11259.58 11243.35 13596.68 9083.702 16906.67
## [,1321] [,1322] [,1323] [,1324] [,1325] [,1326] [,1327] [,1328]
## [1,] 16958.4 11286.91 9858.3 13697.02 10782.21 17040.99 10884.62 13742
## [,1329] [,1330] [,1331] [,1332] [,1333] [,1334] [,1335] [,1336]
## [1,] 9711.172 9519.474 17121.91 10768.35 15396.1 9431.545 17146.73 10962.48
## [,1337] [,1338] [,1339] [,1340] [,1341] [,1342] [,1343] [,1344]
## [1,] 9306.664 14722.93 17184.46 10812.09 9378.479 11064.15 14564.95 17240.53
## [,1345] [,1346] [,1347] [,1348] [,1349] [,1350] [,1351] [,1352]
## [1,] 9397.04 18074.47 9542.582 11040.23 14795.94 9604.127 16789.93 9406.561
## [,1353] [,1354] [,1355] [,1356] [,1357] [,1358] [,1359] [,1360]
## [1,] 10987.33 14823.7 16316.56 9491.1 15078.55 10864.18 16372.25 9318.369
## [,1361] [,1362] [,1363] [,1364] [,1365] [,1366] [,1367] [,1368]
## [1,] 9531.504 15323.66 10875.33 16371.94 15630.27 9772.661 11195.02 16476.06
## [,1369] [,1370] [,1371] [,1372] [,1373] [,1374] [,1375] [,1376]
## [1,] 9787.763 11286.5 9686.378 15626.93 16503.47 9807.055 11398.66 13367.29
## [,1377] [,1378] [,1379] [,1380] [,1381] [,1382] [,1383] [,1384]
## [1,] 11001 16833.25 9208.3 11056 13349.32 16564.82 13436.4 9396.533
## [,1385] [,1386] [,1387] [,1388] [,1389] [,1390] [,1391] [,1392]
## [1,] 8379.067 16598.85 13507.87 9099.606 8406.84 16587.12 9644.663 16593.77
## [,1393] [,1394] [,1395] [,1396] [,1397] [,1398] [,1399] [,1400]
## [1,] 16611.14 9850.75 9427.324 13745.62 8685.474 16630.22 9813.937 9414.515
## [,1401] [,1402] [,1403] [,1404] [,1405] [,1406] [,1407] [,1408]
## [1,] 13703.35 8675.91 16539.96 9433.587 9751.757 13833.02 16491.75 16490.79
## [,1409] [,1410] [,1411] [,1412] [,1413] [,1414] [,1415] [,1416]
## [1,] 8796.663 9669.042 16532.49 8628.702 5926.036 8927.235 10061.69 16515.72
## [,1417] [,1418] [,1419] [,1420] [,1421] [,1422] [,1423] [,1424]
## [1,] 10111.77 8874.36 16533.72 9124.112 10039.51 16533.94 9198.125 9991.776
## [,1425] [,1426] [,1427] [,1428] [,1429] [,1430] [,1431] [,1432]
## [1,] 9165.729 16552.31 9259.661 9093.062 10136.38 16564.1 9319.609 9133.396
## [,1433] [,1434] [,1435] [,1436] [,1437] [,1438] [,1439] [,1440]
## [1,] 10191.9 16287.18 8906.651 10305.93 16334.26 8891.328 16388.82 8872.866
## [,1441] [,1442] [,1443] [,1444] [,1445] [,1446] [,1447] [,1448]
## [1,] 5957.923 16406.97 9083.567 9157.315 8877.013 9257.309 9752.031 16432.33
## [,1449] [,1450] [,1451] [,1452] [,1453] [,1454] [,1455] [,1456]
## [1,] 8964.031 16410.7 9919.3 16455.98 9869.241 16449.89 9981.005 16497.42
## [,1457] [,1458] [,1459] [,1460] [,1461] [,1462] [,1463] [,1464]
## [1,] 10013.92 16372.94 9753.508 16438.13 9663.008 16430.4 16457.21 9775.127
## [,1465] [,1466] [,1467] [,1468] [,1469] [,1470] [,1471] [,1472]
## [1,] 16495.28 9882.128 16495.13 16508.73 9765.522 16544.03 9928.909 16705.03
## [,1473] [,1474] [,1475] [,1476] [,1477] [,1478] [,1479] [,1480]
## [1,] 16753.2 9861.308 10013.2 16731.54 10549.02 16775.64 10452.46 16798.59
## [,1481] [,1482] [,1483] [,1484] [,1485] [,1486] [,1487] [,1488]
## [1,] 10647.17 16799.15 10650.3 16806.38 10300.1 16810.39 10382.89 16862.87
## [,1489] [,1490] [,1491] [,1492] [,1493] [,1494] [,1495] [,1496]
## [1,] 10249.99 10447.59 16568.74 10475.95 16610.93 10278.29 16625.3 16679.13
## [,1497] [,1498] [,1499] [,1500] [,1501] [,1502] [,1503] [,1504]
## [1,] 10590.34 16672.32 16699.26 15743.5 10742.62 16043.6 16052.11 10638
## [,1505] [,1506] [,1507] [,1508] [,1509] [,1510] [,1511] [,1512]
## [1,] 16082 16154.29 10541.57 16168.78 10769.75 16179.47 10569.89 10709.11
## [,1513] [,1514] [,1515] [,1516] [,1517] [,1518] [,1519] [,1520]
## [1,] 16210.79 16170.36 15686.99 15892.63 15874.34 15854.68 15897.62 16862.14
## [,1521] [,1522] [,1523] [,1524] [,1525] [,1526] [,1527] [,1528]
## [1,] 15806.98 15845.09 15845.83 15915.24 15913.71 15935.08 16059.19 16100.92
## [,1529] [,1530] [,1531] [,1532] [,1533] [,1534] [,1535] [,1536]
## [1,] 16125.8 16129.13 16133.91 16163.01 16173.67 16176.93 16166.48 16126.75
## [,1537] [,1538] [,1539] [,1540] [,1541] [,1542] [,1543] [,1544]
## [1,] 16129.75 16147.65 16172.38 16178.27 8108.037 6155.91 9008.376 8224.367
## [,1545] [,1546] [,1547] [,1548] [,1549] [,1550] [,1551] [,1552]
## [1,] 7571.197 6947.599 9024.55 6305.776 7019.352 8080.427 7658.726 8964.112
## [,1553] [,1554] [,1555] [,1556] [,1557] [,1558] [,1559] [,1560]
## [1,] 8185.577 7556.199 9045.759 6056.234 7052.192 8177.24 7553.122 8981.785
## [,1561] [,1562] [,1563] [,1564] [,1565] [,1566] [,1567] [,1568]
## [1,] 8104.684 5936.744 7048.971 7393.875 8139.468 8876.776 7420.002 5931.787
## [,1569] [,1570] [,1571] [,1572] [,1573] [,1574] [,1575] [,1576]
## [1,] 9061.077 7052.066 7323.343 8172.984 5897.592 8957.344 6923.946 7414.143
## [,1577] [,1578] [,1579] [,1580] [,1581] [,1582] [,1583] [,1584]
## [1,] 8135.431 6016.413 8836.42 7310.506 7072.908 8291.343 8889.2 7372.305
## [,1585] [,1586] [,1587] [,1588] [,1589] [,1590] [,1591] [,1592]
## [1,] 8860.607 7952.145 5887.566 7675.605 8851.481 7479.848 7825.435 8050.724
## [,1593] [,1594] [,1595] [,1596] [,1597] [,1598] [,1599] [,1600]
## [1,] 6003.047 7344.285 8918.401 7948.147 7459.433 6296.407 8833.403 7898.317
## [,1601] [,1602] [,1603] [,1604] [,1605] [,1606] [,1607] [,1608]
## [1,] 7473.802 7828.29 7350.856 7369.957 6210.664 7994.953 8218.894 7283.246
## [,1609] [,1610] [,1611] [,1612] [,1613] [,1614] [,1615] [,1616]
## [1,] 7289.328 7956.653 8394.043 7379.208 7458.06 7928.077 6122.087 8077.56
## [,1617] [,1618] [,1619] [,1620] [,1621] [,1622] [,1623] [,1624]
## [1,] 7311.863 6063.954 8441.078 7517.771 7986.268 8315.409 7572.491 7667.348
## [,1625] [,1626] [,1627] [,1628] [,1629] [,1630] [,1631] [,1632]
## [1,] 8001.117 6074.734 7607.257 8079.041 8278.077 7651.641 6223.259 8570.698
## [,1633] [,1634] [,1635] [,1636] [,1637] [,1638] [,1639] [,1640]
## [1,] 7613.676 7678.821 6171.523 8020.735 7653.136 6162.296 7503.881 7691.571
## [,1641] [,1642] [,1643] [,1644] [,1645] [,1646] [,1647] [,1648]
## [1,] 7627.456 8516.878 6208.457 7487.49 7904.993 7644.162 7350.481 7957.419
## [,1649] [,1650] [,1651] [,1652] [,1653] [,1654] [,1655] [,1656]
## [1,] 7405.893 5953.943 7130.634 8603.593 7328.044 5910.683 5967.578 7259.428
## [,1657] [,1658] [,1659] [,1660] [,1661] [,1662] [,1663] [,1664]
## [1,] 7727.479 7935.794 8526.442 7674.385 7006.897 7424.4 8448.832 7723.771
## [,1665] [,1666] [,1667] [,1668] [,1669] [,1670] [,1671] [,1672]
## [1,] 5949.022 7265.747 7300.891 8548.904 7758.205 5868.625 6982.9 7532.566
## [,1673] [,1674] [,1675] [,1676] [,1677] [,1678] [,1679] [,1680]
## [1,] 8481.011 7822.177 7228.156 7709.743 5826.209 7886.296 8562.993 7574.573
## [,1681] [,1682] [,1683] [,1684] [,1685] [,1686] [,1687] [,1688]
## [1,] 7079.868 8053.028 7619.479 7223.364 5859.477 7797.644 7936.791 7711.187
## [,1689] [,1690] [,1691] [,1692] [,1693] [,1694] [,1695] [,1696]
## [1,] 7503.22 7358.198 7985.352 6460.679 7471.032 7858.26 7978.33 6627.304
## [,1697] [,1698] [,1699] [,1700] [,1701] [,1702] [,1703] [,1704]
## [1,] 7670.661 8000.476 7612.425 7892.145 6541.86 7551.953 7490.411 7997.959
## [,1705] [,1706] [,1707] [,1708] [,1709] [,1710] [,1711] [,1712]
## [1,] 6038.878 7790.196 8008.21 7667.426 7577.372 7829.008 5804.509 7486.875
## [,1713] [,1714] [,1715] [,1716] [,1717] [,1718] [,1719] [,1720]
## [1,] 7885.888 7948.006 7583.043 6531.807 7990.859 7441.26 6616.817 7909.204
## [,1721] [,1722] [,1723] [,1724] [,1725] [,1726] [,1727] [,1728]
## [1,] 7334.924 6822.359 7947.494 7606.633 6259.635 7829.252 8004.285 7435.135
## [,1729] [,1730] [,1731] [,1732] [,1733] [,1734] [,1735] [,1736]
## [1,] 7255.933 7904.423 6678.494 6693.301 7332.55 7898.96 6702.028 7948.979
## [,1737] [,1738] [,1739] [,1740] [,1741] [,1742] [,1743] [,1744]
## [1,] 7620.914 7458.717 7902.622 7403.448 7614.732 8103.377 6797.79 7292.642
## [,1745] [,1746] [,1747] [,1748] [,1749] [,1750] [,1751] [,1752]
## [1,] 7468.79 7769.047 7329.36 6702.03 8070.793 7558.577 7983.368 7234.921
## [,1753] [,1754] [,1755] [,1756] [,1757] [,1758] [,1759] [,1760]
## [1,] 6734.576 8012.6 7282.112 6638.45 7943.704 7202.668 7125.488 7198.716
## [,1761] [,1762] [,1763] [,1764] [,1765] [,1766] [,1767] [,1768]
## [1,] 6635.528 7963.331 7085.213 7964.702 7197.719 7357.665 6652.078 7945.037
## [,1769] [,1770] [,1771] [,1772] [,1773] [,1774] [,1775] [,1776]
## [1,] 7120.418 8017.888 6897.564 7074.672 6556.441 7964.737 8056.615 7121.138
## [,1777] [,1778] [,1779] [,1780] [,1781] [,1782] [,1783] [,1784]
## [1,] 7938.383 8073.609 7073.314 6735.852 6600.825 7957.452 8003.39 7171.879
## [,1785] [,1786] [,1787] [,1788] [,1789] [,1790] [,1791] [,1792]
## [1,] 6575.119 7972.61 8043.72 6903.66 6598.113 6625.822 7202.745 7756.148
## [,1793] [,1794] [,1795] [,1796] [,1797] [,1798] [,1799] [,1800]
## [1,] 6563.492 7073.216 7982.688 7097.431 7881.96 6682.387 7268.635 6579.522
## [,1801] [,1802] [,1803] [,1804] [,1805] [,1806] [,1807] [,1808]
## [1,] 7733.848 7061.657 6595.664 7915.025 6612.053 7807.25 7910.664 6330.826
## [,1809] [,1810] [,1811] [,1812] [,1813] [,1814] [,1815] [,1816]
## [1,] 6890.624 6230.442 7753.764 7956.57 6318.468 6592.422 6400.083 6395.233
## [,1817] [,1818] [,1819] [,1820] [,1821] [,1822] [,1823] [,1824]
## [1,] 6616.279 7830.035 7964.935 6393.226 7107.305 6659.824 6604.845 6510.238
## [,1825] [,1826] [,1827] [,1828] [,1829] [,1830] [,1831] [,1832]
## [1,] 6597.959 5899.94 7903.826 7950.249 6637.47 6432.995 6219.141 6343.003
## [,1833] [,1834] [,1835] [,1836] [,1837] [,1838] [,1839] [,1840]
## [1,] 7969.657 6544.887 6107.836 6095.53 7831.424 6467.263 7968.959 6193.509
## [,1841] [,1842] [,1843] [,1844] [,1845] [,1846] [,1847] [,1848]
## [1,] 6102.786 7863.848 5518.833 7925.926 5784.435 7794.831 6169.772 6413.157
## [,1849] [,1850] [,1851] [,1852] [,1853] [,1854] [,1855] [,1856]
## [1,] 7949.729 5770.774 5904.876 5209.664 6759.18 7763.558 5803.297 6644.009
## [,1857] [,1858] [,1859] [,1860] [,1861] [,1862] [,1863] [,1864]
## [1,] 7775.266 6194.168 7717.066 5380.639 5832.238 7771.796 6709.162 5133.275
## [,1865] [,1866] [,1867] [,1868] [,1869] [,1870] [,1871] [,1872]
## [1,] 5791.911 6709.403 7831.376 6050.973 7883.729 5316.747 5742.743 6703.835
## [,1873] [,1874] [,1875] [,1876] [,1877] [,1878] [,1879] [,1880]
## [1,] 5798.806 7751.553 6650.645 7642.101 4898.886 7796.103 7502.093 8705.064
## [,1881] [,1882] [,1883] [,1884] [,1885] [,1886] [,1887] [,1888]
## [1,] 6720.315 7923.37 6710.561 4606.796 8747.408 6831.433 7783.961 8608.591
## [,1889] [,1890] [,1891] [,1892] [,1893] [,1894] [,1895] [,1896]
## [1,] 8538.729 5509.388 6475.341 7817.881 6577.457 7546.939 8735.094 7893.331
## [,1897] [,1898] [,1899] [,1900] [,1901] [,1902] [,1903] [,1904]
## [1,] 6605.033 5743.802 6770.822 7480.719 8082.325 7803.838 7551.196 4989.684
## [,1905] [,1906] [,1907] [,1908] [,1909] [,1910] [,1911] [,1912]
## [1,] 6772.3 8091.582 6594.412 6573.498 8724.456 8295.362 4836.176 7813.642
## [,1913] [,1914] [,1915] [,1916] [,1917] [,1918] [,1919] [,1920]
## [1,] 6760.885 8803.804 7372.856 4852.025 6583.868 4441.5 6969.146 7765.823
## [,1921] [,1922] [,1923] [,1924] [,1925] [,1926] [,1927] [,1928]
## [1,] 6426.546 6982.829 7813.696 6812.588 4588.63 6754.067 8040.238 6811.297
## [,1929] [,1930] [,1931] [,1932] [,1933] [,1934] [,1935] [,1936]
## [1,] 6630.563 8080.551 6875.198 4396.424 6719.245 6828.273 8451.443 7933.996
## [,1937] [,1938] [,1939] [,1940] [,1941] [,1942] [,1943] [,1944]
## [1,] 6902.323 4528.131 6825.183 6591.798 8329.418 6801.142 6912.988 6788.793
## [,1945] [,1946] [,1947] [,1948] [,1949] [,1950] [,1951] [,1952]
## [1,] 3731.748 7995.316 7149.43 6756.944 6627.5 6855.062 7530.402 4018.47
## [,1953] [,1954] [,1955] [,1956] [,1957] [,1958] [,1959] [,1960]
## [1,] 6547.103 7833.031 4468.985 7278.107 6833.976 4508.874 7419.504 7139.655
## [,1961] [,1962] [,1963] [,1964] [,1965] [,1966] [,1967] [,1968]
## [1,] 6483.705 6299.72 7372.697 7115.918 4268.048 6188.904 7084.106 7135.879
## [,1969] [,1970] [,1971] [,1972] [,1973] [,1974] [,1975] [,1976]
## [1,] 4316.812 7040.996 6313.371 5823.514 7765.183 7936.878 6529.3 4114.814
## [,1977] [,1978] [,1979] [,1980] [,1981] [,1982] [,1983] [,1984]
## [1,] 6942.562 5909.123 7439.722 7207.202 6516.936 6190.573 7168.05 4142.921
## [,1985] [,1986] [,1987] [,1988] [,1989] [,1990] [,1991] [,1992]
## [1,] 5879.095 7085.788 6948.88 6519.124 5961.926 6985.036 7181.904 6277.365
## [,1993] [,1994] [,1995] [,1996] [,1997] [,1998] [,1999] [,2000]
## [1,] 3706.558 6693.835 6570.336 7120.512 6803.13 6570.427 6940.613 5077.198
## [,2001] [,2002] [,2003] [,2004] [,2005] [,2006] [,2007] [,2008]
## [1,] 6212.337 6589.636 7085.983 5926.169 5022.439 6625.844 6987.223 6471.779
## [,2009] [,2010] [,2011] [,2012] [,2013] [,2014] [,2015] [,2016]
## [1,] 7167.885 6460.91 5785.533 6942.061 6698.731 5297.466 6677.914 7070.654
## [,2017] [,2018] [,2019] [,2020] [,2021] [,2022] [,2023] [,2024]
## [1,] 6061.061 7362.721 6706.999 5330.528 5645.513 6342.138 7497.792 7028.919
## [,2025] [,2026] [,2027] [,2028] [,2029] [,2030] [,2031] [,2032]
## [1,] 6805.522 5535.117 5786.666 5868.092 7438.73 6457.361 5461.418 7115.456
## [,2033] [,2034] [,2035] [,2036] [,2037] [,2038] [,2039] [,2040]
## [1,] 6413.508 7589.29 6039.029 5639.688 5528.285 5199.614 6561.515 7679.699
## [,2041] [,2042] [,2043] [,2044] [,2045] [,2046] [,2047] [,2048]
## [1,] 7078.27 5730.488 6053.559 7096.141 7475.797 6518.168 5964.118 5734.655
## [,2049] [,2050] [,2051] [,2052] [,2053] [,2054] [,2055] [,2056]
## [1,] 5673.205 6154.195 6868.596 5891.667 5905.757 6834.923 6203.695 7374.357
## [,2057] [,2058] [,2059] [,2060] [,2061] [,2062] [,2063] [,2064]
## [1,] 6131.401 6050.13 5606.947 5901.211 7067.56 6076.631 6236.585 5766.457
## [,2065] [,2066] [,2067] [,2068] [,2069] [,2070] [,2071] [,2072]
## [1,] 7539.75 6444.916 6303.501 6057.439 7448.656 5895.221 6246.079 5935.215
## [,2073] [,2074] [,2075] [,2076] [,2077] [,2078] [,2079] [,2080]
## [1,] 7677.301 4474.01 6103.307 6065.953 5960.646 6316.211 7624.755 4209.122
## [,2081] [,2082] [,2083] [,2084] [,2085] [,2086] [,2087] [,2088]
## [1,] 7666.032 5877.525 4262.518 4066.339 6072.23 6477.554 6098.405 6118.099
## [,2089] [,2090] [,2091] [,2092] [,2093] [,2094] [,2095] [,2096]
## [1,] 7626.522 7630.438 4294.777 6260.753 6141.322 6152.871 7682.618 6190.286
## [,2097] [,2098] [,2099] [,2100] [,2101] [,2102] [,2103] [,2104]
## [1,] 5118.753 4115.807 7587.78 6048.271 5190.332 3915.545 6553.951 4856.047
## [,2105] [,2106] [,2107] [,2108] [,2109] [,2110] [,2111] [,2112]
## [1,] 3899.356 5097.574 7258.559 4994.833 4007.447 7356.964 4053.645 5433.202
## [,2113] [,2114] [,2115] [,2116] [,2117] [,2118] [,2119] [,2120]
## [1,] 7322.868 5430.344 5422.251 4154.245 6704.323 5152.169 3954.723 6863.488
## [,2121] [,2122] [,2123] [,2124] [,2125] [,2126] [,2127] [,2128]
## [1,] 5210.564 5157.146 6634.414 4212.689 5160.246 4571.413 5164.785 6995.835
## [,2129] [,2130] [,2131] [,2132] [,2133] [,2134] [,2135] [,2136]
## [1,] 4194.186 5163.35 6616.589 5156.798 6869.926 4602.297 6712.317 5543.006
## [,2137] [,2138] [,2139] [,2140] [,2141] [,2142] [,2143] [,2144]
## [1,] 4589.594 6840.994 5619.047 6225.243 4576.344 6427.219 5439.048 3502.357
## [,2145] [,2146] [,2147] [,2148] [,2149] [,2150] [,2151] [,2152]
## [1,] 6298.694 5663.297 5531.832 4340.61 6203.7 3827.449 6711.158 5381.654
## [,2153] [,2154] [,2155] [,2156] [,2157] [,2158] [,2159] [,2160]
## [1,] 4154.212 3716.099 6466.279 5704.46 3805.593 8337.878 5609.871 3925.942
## [,2161] [,2162] [,2163] [,2164] [,2165] [,2166] [,2167] [,2168]
## [1,] 5461.961 8511.294 5459.813 8648.797 4138.67 8581.468 8604.597 4620.204
## [,2169] [,2170] [,2171] [,2172] [,2173] [,2174] [,2175] [,2176]
## [1,] 3638.257 8560.729 4658.574 3147.567 4669.592 8498.502 3005.539 4672.419
## [,2177] [,2178] [,2179] [,2180] [,2181] [,2182] [,2183] [,2184]
## [1,] 8680.644 3045.666 9130.023 5062.206 2786.574 9229.313 4680.042 2798.104
## [,2185] [,2186] [,2187] [,2188] [,2189] [,2190] [,2191] [,2192]
## [1,] 9293.213 4773.004 2654.972 9152.19 5862.335 2706.891 8673.664 5709.125
## [,2193] [,2194] [,2195] [,2196] [,2197] [,2198] [,2199] [,2200]
## [1,] 2431.034 9220.203 5450.66 2650.264 9377.502 5183.737 9336.006 5444.426
## [,2201] [,2202] [,2203] [,2204] [,2205] [,2206] [,2207] [,2208]
## [1,] 3613.146 9090.845 5324.554 3265.749 8883.854 5448.409 8936.316 3261.042
## [,2209] [,2210] [,2211] [,2212] [,2213] [,2214] [,2215] [,2216]
## [1,] 9000.696 5794.357 9029.124 3577.252 3643.353 5697.944 9011.59 3145.951
## [,2217] [,2218] [,2219] [,2220] [,2221] [,2222] [,2223] [,2224]
## [1,] 5794.571 3378.728 9200.692 5395.19 5479.812 2905.723 9069.623 5731.566
## [,2225] [,2226] [,2227] [,2228] [,2229] [,2230] [,2231] [,2232]
## [1,] 2912.805 8944.997 5589.474 3660.517 3753.368 8873.095 5569.287 8992.673
## [,2233] [,2234] [,2235] [,2236] [,2237] [,2238] [,2239] [,2240]
## [1,] 3515.713 4756.651 4701.673 3568.313 9027.399 4635.12 3813.002 8946.04
## [,2241] [,2242] [,2243] [,2244] [,2245] [,2246] [,2247] [,2248]
## [1,] 4665.338 3754.703 8997.252 4366.828 3705.989 8880.701 4504.444 9055.173
## [,2249] [,2250] [,2251] [,2252] [,2253] [,2254] [,2255] [,2256]
## [1,] 3764.926 8960.283 3421.84 8899.913 3514.245 3863.061 3755.959 8882.921
## [,2257] [,2258] [,2259] [,2260] [,2261] [,2262] [,2263] [,2264]
## [1,] 3826.001 8875.779 8931.674 3650.269 8922.237 3708.619 8751.131 3761.692
## [,2265] [,2266] [,2267] [,2268] [,2269] [,2270] [,2271] [,2272]
## [1,] 8802.974 3481.272 8758.028 3404.119 3802.827 8799.521 3743.231 8824.57
## [,2273] [,2274] [,2275] [,2276] [,2277] [,2278] [,2279] [,2280]
## [1,] 8779.564 3636.099 3649.3 8784.065 3752.062 3549.711 8870.417 8732.103
## [,2281] [,2282] [,2283] [,2284] [,2285] [,2286] [,2287] [,2288]
## [1,] 3628.469 8759.248 3652.529 3615.955 8781.972 3920.862 8775.66 3881.49
## [,2289] [,2290] [,2291] [,2292] [,2293] [,2294] [,2295] [,2296]
## [1,] 4040.251 8612.824 8460.05 3987.763 8509.827 3982.354 8461.995 4025.4
## [,2297] [,2298] [,2299] [,2300] [,2301] [,2302] [,2303] [,2304]
## [1,] 8563.792 3962.706 8529.768 3230.964 8498.258 2243.443 8518.203 2064.964
## [,2305] [,2306] [,2307] [,2308] [,2309] [,2310] [,2311] [,2312]
## [1,] 8941.512 8881.303 2238.431 8869.548 2273.141 8915.876 1978.37 8868.54
## [,2313] [,2314] [,2315] [,2316] [,2317] [,2318] [,2319] [,2320]
## [1,] 8960.996 2266.084 8884.251 8874.303 2400.176 8819.403 2365.739 8809.142
## [,2321] [,2322] [,2323] [,2324] [,2325] [,2326] [,2327] [,2328]
## [1,] 2213.279 2094.241 8789.527 2063.047 1920.719 8840.769 1924.394 8830.541
## [,2329] [,2330] [,2331] [,2332] [,2333] [,2334] [,2335] [,2336]
## [1,] 1834.075 8777.57 1729.497 8774.17 1803.968 1431.964 8843.842 8785.823
## [,2337] [,2338] [,2339] [,2340] [,2341] [,2342] [,2343] [,2344]
## [1,] 1260.483 8711.832 1352.882 1588.659 8712.91 8861.658 8854.685 1583.726
## [,2345] [,2346] [,2347] [,2348] [,2349] [,2350] [,2351] [,2352]
## [1,] 8891.366 1566.539 1588.826 8872.036 1646.296 8815.209 1664.414 8891.887
## [,2353] [,2354] [,2355] [,2356] [,2357] [,2358] [,2359] [,2360]
## [1,] 8780.284 1680.64 8690.726 8697.431 2214.051 8688.671 8685.151 1506.157
## [,2361] [,2362] [,2363] [,2364] [,2365] [,2366] [,2367] [,2368]
## [1,] 8690.451 1729.179 8744.109 2108.938 2182.808 2981.559 8587.794 2429.273
## [,2369] [,2370] [,2371] [,2372] [,2373] [,2374] [,2375] [,2376]
## [1,] 8747.975 8738.014 2604.876 8706.228 2618.151 8658.761 2848.568 8752.388
## [,2377] [,2378] [,2379] [,2380] [,2381] [,2382] [,2383] [,2384]
## [1,] 8761.065 2830.953 8721.018 8752.842 2843.434 8756.882 2734.043 2242.88
## [,2385] [,2386] [,2387] [,2388] [,2389] [,2390] [,2391] [,2392]
## [1,] 8558.118 2752.517 8591.415 8669.122 2392.92 8504.931 2686.892 8558.345
## [,2393] [,2394] [,2395] [,2396] [,2397] [,2398] [,2399] [,2400]
## [1,] 2170.307 8461.912 2265.222 8503.376 2335.729 8402.681 2293.039 8366.043
## [,2401] [,2402] [,2403] [,2404] [,2405] [,2406] [,2407] [,2408]
## [1,] 8338.69 2268.236 2388.533 8370.633 2348.644 8380.483 2400.977 8397.022
## [,2409] [,2410] [,2411] [,2412] [,2413] [,2414] [,2415] [,2416]
## [1,] 561.7109 8443.465 695.8718 8121.881 709.7098 8277.165 763.5075 8128.452
## [,2417] [,2418] [,2419] [,2420] [,2421] [,2422] [,2423] [,2424]
## [1,] 570.7102 8168.674 8076.196 498.5176 8107.186 3933.577 8198.418 3022.578
## [,2425] [,2426] [,2427] [,2428] [,2429] [,2430] [,2431] [,2432]
## [1,] 8091.826 8335.738 3232.825 8046.74 3832.939 3838.294 8220.702 8124.197
## [,2433] [,2434] [,2435] [,2436] [,2437] [,2438] [,2439] [,2440]
## [1,] 8209.571 707.6884 8300.286 637.0094 8340.081 708.2794 8290.815 7766.226
## [,2441] [,2442] [,2443] [,2444] [,2445] [,2446] [,2447] [,2448]
## [1,] 534.0305 645.6118 7848.867 7777.682 467.8983 7963.564 3686.139 8007.56
## [,2449] [,2450] [,2451] [,2452] [,2453] [,2454] [,2455] [,2456]
## [1,] 708.6258 706.2914 7729.044 544.5524 7855.791 7984.605 381.3328 430.4342
## [,2457] [,2458] [,2459] [,2460] [,2461] [,2462] [,2463] [,2464]
## [1,] 8024.985 353.9668 7954.237 336.499 662.6959 7966.751 726.7708 8175.168
## [,2465] [,2466] [,2467] [,2468] [,2469] [,2470] [,2471] [,2472]
## [1,] 812.0797 894.0105 8272.637 959.036 1290.352 8129.15 1325.148 8287.21
## [,2473] [,2474] [,2475] [,2476] [,2477] [,2478] [,2479] [,2480]
## [1,] 1428.622 8310.654 1231.807 8151.383 891.9199 8162.591 933.993 8354.321
## [,2481] [,2482] [,2483] [,2484] [,2485] [,2486] [,2487] [,2488]
## [1,] 981.7504 987.2562 8515.329 1032.654 8410.149 1091.888 8546.497 1112.964
## [,2489] [,2490] [,2491] [,2492] [,2493] [,2494] [,2495] [,2496]
## [1,] 1149.832 8392.369 8504.732 1036.479 634.2341 8346.885 810.1026 365.7862
## [,2497] [,2498] [,2499] [,2500] [,2501] [,2502] [,2503] [,2504]
## [1,] 8152.128 8244.083 208.2166 8140.539 512.8914 8335.995 201.1693 300.3183
## [,2505] [,2506] [,2507] [,2508] [,2509] [,2510] [,2511] [,2512]
## [1,] 8319.147 903.8711 8221.289 299.9286 8316.76 8345.999 262.1847 8438.817
## [,2513] [,2514] [,2515] [,2516] [,2517] [,2518] [,2519] [,2520]
## [1,] 287.6317 229.8077 8397.917 201.3718 197.4567 8312.108 141.1798 8426.338
## [,2521] [,2522] [,2523] [,2524] [,2525] [,2526] [,2527] [,2528]
## [1,] 64.23811 62.54718 8333.703 8350.195 174.1757 132.2388 8294.289 156.8219
## [,2529] [,2530] [,2531] [,2532] [,2533] [,2534] [,2535] [,2536]
## [1,] 8387.855 198.8418 8378.73 8322.171 8193.219 219.2111 8146.034 231.2689
## [,2537] [,2538] [,2539] [,2540]
## [1,] 230.6937 8217.803 260.3221 87.1663Ademas, podemos calcular la distancia entre un par de coordenadas y todos los puntos de un objeto espacial. Esto puede ser util, como en el ejemplo, para calcular la distancia (en linea recta) entre un punto concreto y todas las ciudades de nuestro objeto espacial.
# creacion del punto
p <- data.frame(lon = 0, lat = 40) %>% st_as_sf(coords = c("lon", "lat")) %>% st_set_crs(4326)
#calculo de la distancia
p_dist <- st_distance(ciudades, p)
dim(p_dist)## [1] 2540 1Areas de influencia (Buffer)
Tambien es posible calcular areas de influencia (buffers) alrededor de una geometria espacial (utiliza las unidades del sistema de coordenadas del objeto espacial). Si queremos calcular un area de influencia de 10000 metros alrededor de un poligono concreto del objeto paises, usaremos la funcion st_buffer()
# buffer de 10 kilometros alrededor de islandia
buf_ic <- paises %>% dplyr::filter(name == "Iceland") %>% st_transform(3055) %>% st_buffer(10000)
plot(buf_ic)## Warning: plotting the first 9 out of 95 attributes; use max.plot = 95 to plot
## all
Si por el contrario, queremos estimar un buffer interior, simplemente debemos indicar la distancia con un numero negativo.
## buffer de 10 kilometros al interior de Islandia
buf_ic <- paises %>% dplyr::filter(name == "Iceland") %>% st_transform(3055) %>% st_buffer(-10000)plot(buf_ic)## Warning: plotting the first 9 out of 95 attributes; use max.plot = 95 to plot
## all
Cuando se usa un sistema decoordenadas decimales se obtiene un aviso de que el calculo ouede ser menos exacto o correcto. Por eso, proyectamos en este caso a un sistema UTM con la funcion st_transform(). Por defecto, la unidad usada en la distancia aqui tambien es la del sistema de corrdenadas.
Cambiar el tipo de geometria
En objetos vectoriales es muy facial cambiar entre distintos tipos de geometria (poligono a linea, linea a puntos, etc). con st_cast(). Esta funcion, combinada con las de las metricas espaciales, puede dar resultados interesantes para el analisis espacial.
# convertimos los poligonos de los paises en lineas
pol_to_line <- st_cast(paises, "MULTILINESTRING")
#perimetros de los poligonos
perim <- st_cast(paises, "MULTILINESTRING") %>% st_length()
head(perim)## Units: [m]
## [1] 54827173 7584941 37353542 5901387 8069398 14598644Transformacion entre formatos: sp, sf, data.frame
Es posible que en ocaciones nos haga falta convertir entre la clase de sp, sf o data.frame. Veamos algunos ejemplos: podemos eliminar la geometria original de nuestro objeto paises, que es de poligonos, lo cual lo convierte en un simple data.frame que ya no podamos tratar como un objeto espacial porque hemos perdido las coordenadas (lon, lat), asi que este objeto es irreversible.
# eliminamos la geometria original
paises_df <- dplyr::select(paises, area, name) %>% st_set_geometry(NULL)
# data.frame (es irreversible sin lon, lat)
head(paises_df)## area name
## 1 1879827.4 [km2] Indonesia
## 2 327884.9 [km2] Malaysia
## 3 736592.7 [km2] Chile
## 4 1086810.6 [km2] Bolivia
## 5 1289866.4 [km2] Peru
## 6 2784305.9 [km2] ArgentinaSin embarto, podemos haceresa conversion manteniendo las columnas de las coordenadas con la funcion as.data.frame().
# convertir en data.frame pero manteniendo la columna de geometria
paises_df2 <- as.data.frame(paises)
class(paises_df2)## [1] "data.frame"Ahora podemos aplicar la,conversion de formatos con st_as_sf(), siempre que primero convirtamos el data.frame que acabamos de crear en un objeto espacial con as_Spacial().
# sf a sp
paises_sp <- as_Spatial(paises)
class(paises_sp)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"# sp a sf
paises_sf <- st_as_sf(paises_sp)
class(paises_sf)## [1] "sf" "data.frame"Exportacion de datos
Tanto con rgdal como con sf, las opciones para exportar datos vectoriales son muy amplias. En rgdal se utiliza la funcion writeORG() para escribir en archivo en diferentes formatos, que vienen definidos por el parametro driver. Dependiendo de la configuracion local y de las posibilidades de ciertos formatos (algunos son de solo lectura), se pueden exportar a unos o a otros,. La funcion orgDrivers() muestra un listado con todos los tipos de archivos vectoriales que rgdal es capaz de leer y escribir.
# exportacion a shapefile de ESRI
# writeOGR(paises_sp[,-95],
# "./results/paises",
# layer = 'paises',
# driver = "ESRI Shapefile")