Obtención de información de un shapefile tipo punto mediante maptools en lenguaje R

Hace algún tiempo escribí un post sobre un uso puntual de la librería maptools de R. Hoy vamos a retomar el tema incluyendo otros aspectos. Para tenerla disponible en R hay que instalar en GNU/Linux (Debian o Ubuntu) los archivos r-cran-maptools, r-cran-sp y libgpcl-dev (vía synaptic o cónsola con apt-get o aptitude). Después de lanzado el intérprete, se ejecutan los siguientes comandos:

library(maptools)
gpclibPermit()

antes de cargar cualquier shapefile de puntos mediante readShapePoints(). Entre los paréntesis se tiene que tener, entre comillas, la ruta al archivo shape. Supongamos que en el directorio SHAPES, de mi home de usuario, incluyo el set correspondiente a sites.shp, entonces el comando siguiente en R cargará al intérprete la información relevante para operar con él:

getinfo.shape('~/Escritorio/SHAPES/sites.shp') #verifica que es un shape de puntos
Shapefile type: Point, (1), # of Shapes: 42
Pt<-readShapePoints('~/Escritorio/SHAPES/sites.shp', proj4string=CRS("+proj=utm +zone=12 +datum=WGS84"))
Pt
> Pt
           coordinates ID  COVER
0  (455552.4, 4641822)  1 shrubs
1  (449959.8, 4633803)  2  trees
2  (441201.7, 4619030)  3  rocks
3  (468214.9, 4608689)  4  grass
4  (456713.1, 4617974)  5 shrubs
.
.
.
37 (470536.3, 4639817) 38  water
38 (470958.4, 4651635) 39  water
39 (474651.6, 4623039) 40  grass
40 (428117.1, 4641083) 41  trees
41 (443628.6, 4593494) 42  grass

La información entre paréntesis son las coordenadas x,y pero no están en su tabla atributiva. Solo está la correspondiente al ID y a COVER. La información de estas columnas puede ser accedida como si fuese un data frame pero en realidad no lo es. A continuación la manera de corroborarlo:

> is.data.frame(Pt)
[1] FALSE
> Pt$coordinates
NULL
> Pt$ID
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
[26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
> Pt$COVER
 [1] shrubs trees  rocks  grass  shrubs trees  bare   bare   trees  grass 
[11] grass  rocks  grass  trees  shrubs bare   grass  bare   trees  trees 
[21] shrubs trees  rocks  grass  rocks  trees  shrubs grass  bare   rocks 
[31] grass  rocks  shrubs trees  bare   trees  grass  water  water  grass 
[41] trees  grass 
Levels: bare grass rocks shrubs trees water

Para poder acceder y manipular la información de las coordenadas hay que convertir a Pt en un data frame. Para ello:

> Pt2<-as.data.frame(Pt)
> is.data.frame(Pt2)
[1] TRUE
> Pt2
   ID  COVER coords.x1 coords.x2
0   1 shrubs  455552.4   4641822
1   2  trees  449959.8   4633803
2   3  rocks  441201.7   4619030
3   4  grass  468214.9   4608689
4   5 shrubs  456713.1   4617974
.
.
.
37 38  water  470536.3   4639817
38 39  water  470958.4   4651635
39 40  grass  474651.6   4623039
40 41  trees  428117.1   4641083
41 42  grass  443628.6   4593494
> Pt2$coords.x1
 [1] 455552.4 449959.8 441201.7 468214.9 456713.1 472330.2 456396.6 480033.1
 [9] 477395.1 491429.3 483093.2 490690.7 488580.3 488791.3 482671.1 470114.2
[17] 467898.3 473068.8 479716.6 485098.1 478766.9 474018.5 457029.7 443312.1
[25] 435503.6 432549.0 445000.4 456607.6 451437.1 445422.5 440779.6 441096.1
[33] 448377.0 458190.4 447110.8 462938.8 466210.0 470536.3 470958.4 474651.6
[41] 428117.1 443628.6
> Pt2$coords.x2
 [1] 4641822 4633803 4619030 4608689 4617974 4619874 4627577 4628104 4633169
[10] 4635702 4642561 4646782 4629582 4616286 4603307 4603518 4597292 4592966
[19] 4591700 4596448 4601091 4603518 4603413 4599614 4606051 4612276 4611432
[28] 4611221 4619452 4628949 4638023 4647837 4653007 4652691 4646148 4634541
[37] 4625783 4639817 4651635 4623039 4641083 4593494
> mean(Pt2$coords.x1) #coordenada x del centroide de los datos
[1] 462647.4
> mean(Pt2$coords.x2) #coordenada y del centroide de los datos
[1] 4622150

Para obtener un resumen de nuestro shapefile se hace un summary(Pt) que, entre otras cosas, señala las coordenadas del cuadrilátero que contiene los datos, la proyección, el número de rasgos y un resumen de los atributos.

> summary(Pt)
Object of class SpatialPointsDataFrame
Coordinates:
                min       max
coords.x1  428117.1  491429.3
coords.x2 4591699.9 4653007.2
Is projected: TRUE 
proj4string : [+proj=utm +zone=12 +datum=WGS84]
Number of points: 42
Data attributes:
       ID           COVER   
 Min.   : 1.00   bare  : 6  
 1st Qu.:11.25   grass :11  
 Median :21.50   rocks : 6  
 Mean   :21.50   shrubs: 6  
 3rd Qu.:31.75   trees :11  
 Max.   :42.00   water : 2
> plot(Pt, col='blue')  

Un plot(Pt) permite visualizar la distribución espacial de los datos.

maptools

Existen diferentes maneras para disponer de los valores de las coordenadas para incluirlas en la tabla atributiva del shapefile. Una de las más sencillas es crear una lista con las variables a copiar, etiquetar las columnas y escribir a un archivo de texto con write.table(). El proceso sería:

> Pt3<-list(Pt2$coords.x1, Pt2$coords.x2)
names(Pt3)<-c('X','Y')
> Pt3
$X
 [1] 455552.4 449959.8 441201.7 468214.9 456713.1 472330.2 456396.6 480033.1
 [9] 477395.1 491429.3 483093.2 490690.7 488580.3 488791.3 482671.1 470114.2
[17] 467898.3 473068.8 479716.6 485098.1 478766.9 474018.5 457029.7 443312.1
[25] 435503.6 432549.0 445000.4 456607.6 451437.1 445422.5 440779.6 441096.1
[33] 448377.0 458190.4 447110.8 462938.8 466210.0 470536.3 470958.4 474651.6
[41] 428117.1 443628.6

$Y
 [1] 4641822 4633803 4619030 4608689 4617974 4619874 4627577 4628104 4633169
[10] 4635702 4642561 4646782 4629582 4616286 4603307 4603518 4597292 4592966
[19] 4591700 4596448 4601091 4603518 4603413 4599614 4606051 4612276 4611432
[28] 4611221 4619452 4628949 4638023 4647837 4653007 4652691 4646148 4634541
[37] 4625783 4639817 4651635 4623039 4641083 4593494
write.table(Pt3,'datos.out',col.names=TRUE,quote=FALSE)

Un extracto del archivo es:

X Y
1 455552.418360864 4641822.05368488
2 449959.840851334 4633802.50857687
3 441201.65343075 4619029.66232529
4 468214.858005083 4608688.66994917
.
.
.
38 470536.305273189 4639817.16740788
39 470958.386594663 4651635.44440915
40 474651.59815756 4623039.43487929
41 428117.132465057 4641083.4113723
42 443628.621029225 4593493.74237611

Otra manera equivalente de hacer esto es con cbind():

Pt4<-cbind(Pt2$coords.x1, Pt2$coords.x2)
write.table(Pt4,'datos.out',col.names=TRUE,quote=FALSE)

con la diferencia de que los nombres de las columnas son V1 y V2; respectivamente:

V1 V2
1 455552.418360864 4641822.05368488
2 449959.840851334 4633802.50857687
3 441201.65343075 4619029.66232529
4 468214.858005083 4608688.66994917
.
.
.
38 470536.305273189 4639817.16740788
39 470958.386594663 4651635.44440915
40 474651.59815756 4623039.43487929
41 428117.132465057 4641083.4113723
42 443628.621029225 4593493.74237611

Si queremos modificar los nombres se puede convertir en data frame.

Esta entrada fue publicada en Lenguaje R, SIG, Software Libre. Guarda el enlace permanente.

2 respuestas a Obtención de información de un shapefile tipo punto mediante maptools en lenguaje R

  1. Pingback: Obtención de información de un sh...

  2. Pingback: Obtención de información de un sh...

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s