Cartografía Básica

Aquí empezaremos con algunos ejercicios para trabajar con objetos espaciales y hacer algunos mapas b?sicos.

Vamos a generar unos datos aleatorios de coordenadas geográficas y de precipitación.

name <- LETTERS[1:10]
longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5,
-120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9,
36.2, 39, 41.6, 36.9)
stations <- cbind(longitude, latitude)
# Simulated rainfall data
set.seed(0)
# runif: genera valores aleatorios siguiendo una distribuci?n uniforme
precip <- (runif(length(latitude))*10)^3

Ahora vamos a graficarlos y el tamaño de los puntos será proporcional a la cantidad de precipitación.

psize <- 1 + precip/500
plot(stations, cex=psize, pch=20, col='red', main='Precipitation')

# add names to plot
text(stations, name, pos=4)

# add a legend
breaks <- c(100, 500, 1000, 2000)
legend("topright", legend=breaks, pch=20, pt.cex=psize, col='red', bg='gray')

Puntos espaciales:

Vamos a generar vectores con datos de longitud y latitud, y luego los combinamos.

longitude <- c(-116.7, -120.4, -116.7, -113.5, -115.5, -120.8, -119.5, -113.7, -113.7, -110.7)
latitude <- c(45.3, 42.6, 38.9, 42.1, 35.7, 38.9, 36.2, 39, 41.6, 36.9)
lonlat <- cbind(longitude, latitude)

Vamos a usar el paquete sp. Debemos instalarlo usando la función install.packages(“sp”)

library(sp)
pts <- SpatialPoints(lonlat)

¿Qué clase de objeto es?

class(pts)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"

Exploramos que tiene

showDefault(pts)
## An object of class "SpatialPoints"
## Slot "coords":
##       longitude latitude
##  [1,]    -116.7     45.3
##  [2,]    -120.4     42.6
##  [3,]    -116.7     38.9
##  [4,]    -113.5     42.1
##  [5,]    -115.5     35.7
##  [6,]    -120.8     38.9
##  [7,]    -119.5     36.2
##  [8,]    -113.7     39.0
##  [9,]    -113.7     41.6
## [10,]    -110.7     36.9
## 
## Slot "bbox":
##              min    max
## longitude -120.8 -110.7
## latitude    35.7   45.3
## 
## Slot "proj4string":
## CRS arguments: NA

El sistema de coordenadas de referencia de nuestro objeto espacial no está definido. Vamos a definirlo con la ayuda de la función CRS

crdref <- CRS('+proj=longlat +datum=WGS84')
pts <- SpatialPoints(lonlat, proj4string=crdref)

Cargamos la libreria raster

library(raster)
pts
## class       : SpatialPoints 
## features    : 10 
## extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

Vamos a generar (inventar) unos datos de precipitación asociados a esas coordenadas. Estos datos NO son reales, es solo para mostrar como podemos visualizarlos como objetos espaciales.

df <- data.frame(ID=1:nrow(lonlat), precip=(latitude-30)^3)

Combinamos los datos espaciales y los datos de precipitación

ptsdf <- SpatialPointsDataFrame(pts, data=df)
ptsdf
## class       : SpatialPointsDataFrame 
## features    : 10 
## extent      : -120.8, -110.7, 35.7, 45.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       : ID,   precip 
## min values  :  1,  185.193 
## max values  : 10, 3581.577

Vemos que tiene ese objeto espacial

str(ptsdf)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  10 obs. of  2 variables:
##   .. ..$ ID    : int [1:10] 1 2 3 4 5 6 7 8 9 10
##   .. ..$ precip: num [1:10] 3582 2000 705 1772 185 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:10, 1:2] -117 -120 -117 -114 -116 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "longitude" "latitude"
##   ..@ bbox       : num [1:2, 1:2] -120.8 35.7 -110.7 45.3
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "longitude" "latitude"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

También podemos verlos así:

showDefault(ptsdf)
## An object of class "SpatialPointsDataFrame"
## Slot "data":
##    ID   precip
## 1   1 3581.577
## 2   2 2000.376
## 3   3  704.969
## 4   4 1771.561
## 5   5  185.193
## 6   6  704.969
## 7   7  238.328
## 8   8  729.000
## 9   9 1560.896
## 10 10  328.509
## 
## Slot "coords.nrs":
## numeric(0)
## 
## Slot "coords":
##       longitude latitude
##  [1,]    -116.7     45.3
##  [2,]    -120.4     42.6
##  [3,]    -116.7     38.9
##  [4,]    -113.5     42.1
##  [5,]    -115.5     35.7
##  [6,]    -120.8     38.9
##  [7,]    -119.5     36.2
##  [8,]    -113.7     39.0
##  [9,]    -113.7     41.6
## [10,]    -110.7     36.9
## 
## Slot "bbox":
##              min    max
## longitude -120.8 -110.7
## latitude    35.7   45.3
## 
## Slot "proj4string":
## CRS arguments:
##  +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

Líneas y polígonos espaciales. Hacer puntos espaciales es bastante fácil. Hacer líneas y polígonos es un poco más complicado pero vamos a usar las funciones SpatialLines y SpatialPolygons para esto.

lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7)
lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6)
lonlat <- cbind(lon, lat)
lns <- spLines(lonlat, crs=crdref)
lns
## class       : SpatialLines 
## features    : 1 
## extent      : -117.7, -111.9, 37.6, 42.9  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

Graficamos

plot(lns)

pols <- spPolygons(lonlat, crs=crdref)
pols
## class       : SpatialPolygons 
## features    : 1 
## extent      : -117.7, -111.9, 37.6, 42.9  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

Graficamos

plot(pols)

Vamos a ver la estructura del objeto pols

str(pols)
## Formal class 'SpatialPolygons' [package "sp"] with 4 slots
##   ..@ polygons   :List of 1
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. .. .. ..@ Polygons :List of 1
##   .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
##   .. .. .. .. .. .. ..@ labpt  : num [1:2] -114.7 40.1
##   .. .. .. .. .. .. ..@ area   : num 19.7
##   .. .. .. .. .. .. ..@ hole   : logi FALSE
##   .. .. .. .. .. .. ..@ ringDir: int 1
##   .. .. .. .. .. .. ..@ coords : num [1:8, 1:2] -117 -114 -113 -112 -114 ...
##   .. .. .. ..@ plotOrder: int 1
##   .. .. .. ..@ labpt    : num [1:2] -114.7 40.1
##   .. .. .. ..@ ID       : chr "1"
##   .. .. .. ..@ area     : num 19.7
##   ..@ plotOrder  : int 1
##   ..@ bbox       : num [1:2, 1:2] -117.7 37.6 -111.9 42.9
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

Graficamos

plot(pols, axes=TRUE, las=1)
plot(pols, border='blue', col='yellow', lwd=3, add=TRUE)
points(pts, col='red', pch=20, cex=3)

Ahora vamos a manipular algunos archivos tipo raster

El paquete sp permite trabajar con archivos de clase SpatialGridDataFrame y SpatialPixelsDataFrame. Estos archivos son archivos que usamos muy a menudo en operaciones espaciales y SIG. De aquí en adelante nos vamos a empezar a usar el paqute rasterpara trabajar con archivo tipo raster (celdas). Con este paquete podemos hacer operaciones matemáticas con rasters y un sinnúmero de otros procesos que se usan en programas dedicadas a SIG (e.g., ArcGIS, QGIS, etc.).

¿Qué es un capa de tipo raster? Un RasterLayer es un objeto de R que contiene datos de una sola variables. Como vimos en la clase este objeto debe tener información del número de celdas, hileras (rows) y columnas (columns), la extensión geográfica o espacial (x, y - min-, -max-) y lo más importante en SIG: el sistema de coordenadas de referencia.

library(raster)
r <- raster(ncol=10, nrow=10, xmx=-80, xmn=-150, ymn=20, ymx=60)
r
## class       : RasterLayer 
## dimensions  : 10, 10, 100  (nrow, ncol, ncell)
## resolution  : 7, 4  (x, y)
## extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

El objeto r que acabamos de crear es solo el “esqueleto” del raster. Aún no teneoms los datos asociados (e.g., temperatura, precipitación, o cobertura de árboles). Vamos a generar un vector con números aleatorios que llenen las 100 celdas de nuestro raster. Usaremos la función runif.

values(r) <- runif(ncell(r))
r
## class       : RasterLayer 
## dimensions  : 10, 10, 100  (nrow, ncol, ncell)
## resolution  : 7, 4  (x, y)
## extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.01339033, 0.9926841  (min, max)

Ahora, podremos gráficar y ver como sale nuestro raster

plot(r)

Ahora vamos a subir unas capas climáticas que tenemos guardadas en nuestro computador.

capas <- list.files(path="/Users/juvelas/Dropbox/Cursos_Charlas/Zaragosa_2018/Scripts/Capas_clima/", pattern=".asc", full.names = T)
capas <- stack(capas)
capas
## class       : RasterStack 
## dimensions  : 2183, 3799, 8293217, 19  (nrow, ncol, ncell, nlayers)
## resolution  : 0.008333334, 0.008333334  (x, y)
## extent      : -118.3684, -86.71002, 14.5296, 32.72126  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## names       : bio1, bio10, bio11, bio12, bio13, bio14, bio15, bio16, bio17, bio18, bio19, bio2, bio3, bio4, bio5, ...

Ahora vamos a graficar una capa climática (bio1)

plot(capas[[1]])

library(rgdal)
## rgdal: version: 1.3-1, (SVN revision 747)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-7
mx <- readOGR(dsn="/Users/juvelas/Dropbox/Cursos_Charlas/Zaragosa_2018/Scripts/Poligono_Mexico/", layer="gadm36_MEX_1")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/juvelas/Dropbox/Cursos_Charlas/Zaragosa_2018/Scripts/Poligono_Mexico", layer: "gadm36_MEX_1"
## with 32 features
## It has 10 fields
plot(mx)

Vamos a extraer un estado de México y cortaremos las capas usando este estado

estado <- subset(mx, NAME_1=="Chihuahua")
plot(estado)

Cortamos las capas usando mask. También podemos usar crop

# Usando mask
chi.clima <- mask(capas, estado)
plot(chi.clima[[1]])

# Usando crop
chi.clima.v2 <- crop(capas, estado)
plot(chi.clima.v2[[1]])

Eso es todo por hoy!

Adaptado de http://www.rspatial.org/intr/index.html (R. Hijmans, 2016)