El purín de ortigas no mejora la cosecha

Este año hemos publicado un artículo en PeerJ sobre el purín de ortigas (Urtica dioica). El resultado principal del artículo es que el tratamiento con purín de ortigas no funciona para mejorar la cosecha, al menos en patatas en los suelos de huerta de Valencia. Si algún lector ha realizado algún experimento similar, tanto si le sale que funciona, como si no, por favor dejar el enlace en comentarios.

En contra de lo que mucha gente dice, no ha sido difícil publicar un “resultado negativo”. Supongo que tiene que ver con que el diseño experimental esté bien hecho y el resultado sea claro.

El resultado es importante, ya que el tratamiento con purín de ortigas es uno de los tratamientos más utilizados en agricultura ecológica y sin embargo no existía ninguna publicación en la que se evaluase su efecto de forma objetiva.

Para evitar sesgos hemos utilizado un sistema de doble ciego, como se utiliza en medicina. De esta forma ni el agricultor, ni yo mismo mientras realizaba los análisis sabía qué producto era cual. Es decir, el agricultor no podía cuidar más a unas plantas que a otras o a la hora de medir la cosecha no se sabía de qué tratamiento era cada medida. Los códigos utilizados se desvelaron al final. En realidad tampoco hubiese sido necesario, dado que al final no había diferencias significativas entre tratamientos.

 

Diseño experimental en 6 bloques. en cada bloque se han realizado los 6 tratamientos.

 

Los tratamientos que hemos comparado son:

  • Purín de ortigas comercial, a la dosis recomendada por el fabricante.
  • Purín de ortigas comercial, a la mitad de la dosis recomendada por el fabricante.
  • Purín de ortigas comercial, al doble de la dosis recomendada por el fabricante.
  • Purín de ortigas junto con purín de cola de caballo (Equisetum arvense) ambos a la dosis recomendada.
  • Control positivo, con tratamiento de abono foliar convencional.
  • Control negativo, con agua, en la misma cantidad y fechas que los anteriores.

 

Efecto de los tratamientos sobre la cosecha (kg/m) con todos los datos (102 grados de libertad). Los tratamientos son: (A) Purín de ortiga a dosis recomendada (RD); (B) Purín de ortiga a 1/2 RD; (C) Purín de ortiga 2 RD; (D) Purín de Urtica + Equisetum; (E) Abono foliar convencional; (F) control con agua. HSD = 0.6554 kg/m. Las columnas con la misma letra no difieren significativamente en p ≤ 0.05 (HSD).

 

El único efecto que se vio del purín es una pequeña tendencia (no significativa) a mejorar el aspecto de la planta. Pero sin afectar a la cosecha ni a la cantidad de plagas. Es curioso que el tratamiento con abono  foliar convencional tampoco sirva de nada. Probablemente tampoco sea efectivo y sea otra forma de tirar dinero (en este caso para agricultores convencionales).

Un efecto que no esperábamos es un efecto del riego. Debido a la inclinación de la parcela, en la zona más baja, el agua se acumulaba más y la cosecha en esa zona ha sido mayor. Por suerte el diseño en bloques ha servido para poder lidiar con ese efecto. El artículo se puede utilizar como ejemplo de cómo analizar el efecto de una variable anulando (utilizando los residuos) el efecto de otra.

 

Conclusión: Para cultivar patatas en Valencia es mejor no gastar dinero en tratamientos que no funcionan. Si se quiere gastar en algo, mejor aumentar lo que se pueda el riego.

 

Para leer el artículo completo, aquí.

Para obtener el código de R y los datos originales con el que se han realizado los análisis, aquí.

 

Mediación y moderación. Dos formas distintas de interacción.

La mediación y la moderación son dos formas diferentes de interacción de las que no es fácil encontrar buenos ejemplos en la bibliografía. Aprovechando que el tema ha salido en clase, aprovecho para aclararlo por aquí.

Una variable mediadora intenta explicar el mecanismo por el que ocurre el efecto, mientras que una variable moderadora cambia el efecto que tiene la variable independiente sobre la dependiente. Como el post me  está quedando largo, pongo un índice. Primero cuento la moderación y pongo un ejemplo con las especies de iris. Después para sólo una especie, hago un ejemplo de mediación.

Moderación

La moderación  es lo que normalmente se entiende por “interacción” en Estadística. Es cuando una variable moderadora (Z) cambia el efecto que la variable independiente (X) tiene sobre la variable dependiente (Y).

Es decir, suponiendo que X e Y están relacionadas linearmente,

Y = a + b X

la variable moderadora puede cambiar esto ya que b depende de Z.

Y = a +b X + c Z + d XZ

Ejemplo de moderación con R e iris

Lo primero es comprobar el efecto directo de la variable dependiente (Sepal.Length) sobre la variable independiente (Petal.Length).

### direct relationship
lmC <- lm(Petal.Length ~ Sepal.Length, data = iris)
summary(lmC)

### Plot
plot(Petal.Length ~ Sepal.Length, data = iris)
abline(lmC)

Coefficients:
             Estimate Std. Error t value Pr(>|t|) 
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586  21.65 <2e-16 ***

Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583 
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16

 

Se observa una fuerte relación (luego veremos que mal entendida) entre la longitud de los sépalos y de los pétalos. Lo siguiente es comprobar el efecto moderador de la variable Species.

### Moderation effect
lmZ <- lm(Petal.Length ~ Sepal.Length +
as.numeric(Species) +
I(Sepal.Length * as.numeric(Species)),
data = iris)
summary(lmZ)

Coefficients:
                                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)                           -6.91282    0.76764  -9.005 1.08e-15 ***
Sepal.Length                           1.36959    0.14659   9.343  < 2e-16 ***
as.numeric(Species)                    3.12118    0.33548   9.304  < 2e-16 ***
I(Sepal.Length * as.numeric(Species)) -0.29266    0.05889  -4.970 1.85e-06 ***

Residual standard error: 0.3915 on 146 degrees of freedom
Multiple R-squared:  0.9518,	Adjusted R-squared:  0.9508 
F-statistic: 961.1 on 3 and 146 DF,  p-value: < 2.2e-16

 

Se observa que hay interacción y que todas las variables analizadas tienen un efecto significativo por separado y por su interacción. También se puede ver que el modelo es mejor que el anterior, ya que R2 ha aumentado bastante.

Los resultados se pueden representar utilizando base R o con ggplot.

Base R plot

### Partial linear models
set <- subset(iris, Species == "setosa")
ver <- subset(iris, Species == "versicolor")
vir <- subset(iris, Species == "virginica")

lmset <- lm(Petal.Length ~ Sepal.Length, data = set)
lmver <- lm(Petal.Length ~ Sepal.Length, data = ver)
lmvir <- lm(Petal.Length ~ Sepal.Length, data = vir)

summary(lmset)
summary(lmver)
summary(lmvir)

### Plot
plot(Petal.Length ~ Sepal.Length, col = as.numeric(Species), data = iris)
abline(lmset, col = 1)
abline(lmver, col = 2)
abline(lmvir, col = 3)

ggplot

library(ggplot2)
ggplot(data = iris, aes(x = Sepal.Length, y = Petal.Length, col = Species,
shape = Species)) +
# geom_point(size = 3) +
geom_jitter(size = 2) +
geom_smooth(method = 'lm',formula = y ~ x) +
theme_bw()

 

Mediación

Para que exista mediación es necesario que:

  • Haya un efecto C, efecto directo entre la variable independiente (X) y la dependiente (Y).
  • Haya un efecto A entre la variable independiente y la variable mediadora (M).
  • Haya un efecto B entre la variable mediadora y la variable dependiente.

En ese caso la duda es cuanto efecto C queda como C’ (efecto residual) si tenemos en cuenta el efecto que pasa a través de la variable mediadora.

El efecto C lo podemos estimar en un modelo lineal a partir de la pendiente y la probabilidad (p) de que sea al azar. Para decir que hay un efecto p tiene que ser menor de 0.05, que va a ser nuestro nivel de significación.

Modelo 1: Y = a1 + C X , p1 < 0.05

Modelo 2: M = a2 + A X, p2 < 0.05

Modelo 3: Y = a3 + B M, p3 < 0.05

Donde C, A y B representan las pendientes de la recta, cuanto mayor sea, mayor es el efecto de X en Y, de X en M o de M en Y respectivamente.

Una vez que los tres modelos anteriores son significativos (p < 0.05), lo que nos queda por saber es cuando del efecto inicial de X en Y se mantiene al tener en cuente el efecto de M en Y.

Modelo 4: Y = a4 + C’ X + B M,       pc’ = ??

donde se pueden dar dos opciones:

  • que pc’ < 0.05, moderación parcial.
  • que pc’ ≥ 0.05, moderación total.

En ambos casos la pendiente C’ debería ser menor que la pendiente C. En el primer caso, la moderación parcial significa que parte del efecto de X sobre Y, pasa a través de M, pero no todo. En el segundo caso, la moderación total, es que todo el efecto pasa a través de la variable moderadora.

Ejemplo de moderación con R e iris.

Para el ejemplo voy a utilizar los datos “iris” de longitud y anchura de pétalos y sépalos en Iris versicolor. Ya hemos visto en el ejemplo de mediación, que no tendría sentido hacerlo para las tres especies juntas.

Voy a explorar la posibilidad de mediación de la longitud de los pétalos entre la longitud de los sépalos y la anchura de los pétalos. Es decir, viendo si el efecto que tiene la longitud de los sépalos en la anchura de los pétalos es directo o si pasa por la longitud de los pétalos.

Una cosa importante es no confundir la moderación estadística con la moderación causal. El que algo sea significativo estadísticamente no es suficiente para considerar que hay una relación causa – efecto (también pasa con la mediación, la correlación y cualquier otro test estadístico).

Modelo 1

### Iris versicolor
ver <- subset(iris, Species == "versicolor")

### Plot
plot(Petal.Width ~ Sepal.Length, data = ver)
abline(lmC)

### direct relationship
lmC <- lm(Petal.Width ~ Sepal.Length, data = ver)
summary(lmC)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.08326    0.27592   0.302    0.764    
Sepal.Length  0.20936    0.04631   4.521 4.04e-05 ***

Residual standard error: 0.1673 on 48 degrees of freedom
Multiple R-squared:  0.2986,	Adjusted R-squared:  0.284 
F-statistic: 20.44 on 1 and 48 DF,  p-value: 4.035e-05

 

Se puede ver que existe un efecto de la longitud de los sépalos sobre la anchura de los pétalos. El modelo es regularcillo, ya que sólo explica el 28 % de la varianza (R2), aunque es claramente significativo (p = 0,000004).

Modelo 1: Y = 0.083 + 0.209 X

Modelo 2

Para comprobar que hay mediación necesitamos que los otros dos modelos también sean significativos.

### Mediator Petal.Width (A)
lmA <- lm(Petal.Length ~ Sepal.Length, data = ver)
summary(lmA)

### Plot
plot(Petal.Length ~ Sepal.Length, data = ver)
abline(lmA)

Coefficients:
             Estimate Std. Error t value Pr(>|t|) 
(Intercept)  0.18512 0.51421 0.360 0.72 
Sepal.Length 0.68647 0.08631 7.954 2.59e-10 ***

Residual standard error: 0.3118 on 48 degrees of freedom
Multiple R-squared: 0.5686, Adjusted R-squared: 0.5596 
F-statistic: 63.26 on 1 and 48 DF, p-value: 2.586e-10

 

Modelo 2: M = 0.185 + 0.686 X

El modelo 2 es también muy claro (R2 = 0.5596 y p << 0.05).

Modelo 3

lmB <- lm(Petal.Width ~ Petal.Length, data = ver)
summary(lmB)

### Plot
plot(Petal.Width ~ Petal.Length, data = ver)
abline(lmB)

Coefficients:
              Estimate Std. Error t value Pr(>|t|) 
(Intercept)   -0.08429 0.16070 -0.525 0.602 
Petal.Length   0.33105 0.03750  8.828 1.27e-11 ***

Residual standard error: 0.1234 on 48 degrees of freedom
Multiple R-squared: 0.6188, Adjusted R-squared: 0.6109 
F-statistic: 77.93 on 1 and 48 DF, p-value: 1.272e-11

 

Modelo 3: Y = -0.084 + 0.331 M

El modelo 3 es todavía mejor, explicando el 61 % de la varianza.

¿hay mediación o no?

Con los tres modelos anteriores ya sabemos que hay mediación ya que los tres modelos son significativos.

Lo que vamos a ver ahora es hasta donde llega esta mediación o el valor de C’ o el efecto de X sobre Y teniendo en cuenta que parte de dicho efecto pasa a través de M.

### Mediator Petal.Width (C2)
lmC2 <- lm(Petal.Width ~ Sepal.Length + Petal.Length, data = ver)
summary(lmC2)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.01561    0.20446   0.076    0.939    
Sepal.Length -0.04149    0.05218  -0.795    0.430    
Petal.Length  0.36542    0.05731   6.376 7.25e-08 ***

Residual standard error: 0.1238 on 47 degrees of freedom
Multiple R-squared:  0.6239,	Adjusted R-squared:  0.6079 
F-statistic: 38.98 on 2 and 47 DF,  p-value: 1.046e-10

 

Se puede observar que:

  • La pendiente del efecto de la longitud de los sépalos sobre la anchura de los pétalos ha cambiado de
    • de C  = 0.209 , claramente significativa (p = 0,000004), en el modelo 1.
    • a   C’ = -0.041, no significativa (p = 0.43) en este modelo.

Es decir se ha perdido el efecto de la longitud de los sépalos y se puede decir que existe una mediación total por parte de la longitud de los pétalos. Todo el efecto de X sobre Y pasa a través de M.

O dicho de otra manera, el efecto de la longitud de los sépalos sobre la anchura de los pétalos se pierde al tener en cuenta el efecto de la longitud de los pétalos, lo que supone que la longitud de los pétalos afecta a la longitud de los pétalos y es la longitud de los pétalos la que afecta a la anchura de los mismos.

Sólo recuerda que aunque haya mediación estadística, eso no significa que haya mediación causal o que la causa real sea esa.

 

R – primeros pasos (actualizado 18-09-23)

R es el lenguaje de programación de referencia para estadística, aunque tiene el potencial de poder hacer casi cualquier cosa como ocurre con casi cualquier lenguaje de programación. Ha cambiado la forma de entender la forma de analizar los datos, pasando a algo cerrado con algunos tests que hace todo el mundo, a ser algo mucho más creativo. El problema principal es su curva de aprendizaje. No es fácil empezar.

¿Por donde empezamos?

  1. Instalar R
  2. Instalar R-studio
  3. Tutoriales para empezar con R

Instalar R

La web oficial de R no es precisamente amigable, ya que como ocurre  con muchos recursos gratuitos y abiertos están más pensados para los desarrolladores que para los usuarios. Hay que saber qué es cada cosa para no perderse.

R1R2Cuando das a “download R” te envía a una lista de mirrors, que son todos lo mismo, pero en diferentes servidores (así si cae uno, quedan los demás). Para instalar R puedes seguir las instrucciones de cualquiera de los mirrors (se supone que mejor cuanto más cercano), o ir a las FAQ (Frequently asked questions) y ver el tutorial de cómo instalar en tu sistema operativo.R3

Linux, Mac or Windows

El tutorial explica sólo para  la mejor forma de instalarlo,desde repositorios oficiales de R-cran. Si se añade el repositorio, NO OLVIDAR instalar la clave para que no de errores al actualizar.

En ocasiones lo más fácil es utilizar los repositorios oficiales de tu versión Linux, que ya puede tener versiones suficientemente actualizadas de R. Para Ubuntu o Mint  valdría con:

sudo apt-get update
sudo apt-get install r-base r-base-dev

Instalar R – StudioR4

R se puede utilizar directamente desde consola y utilizando archivos de texto plano para hacer los scripts y pegar trozos en la consola. Personalmente me gusta mucho cómo R – studio separa la pantalla en cuatro partes: editor texto, consola, variables y graficos (ayuda, etc), facilitando el correr líneas de código o ver los gráficos en el mismo programa. Lo  único malo es  que no está en los repositorios y que hay que instalarlo  por tanto desde su web oficial. Es tan fácil como dar a “Download RStudio” y luego “Download Opensource RStudio Desktop“, se abre desde el instalador y se instala.

CÓMO EMPEZAR CON R

R5Se supone que si has llegado hasta aquí, ya tienes R y RStudio instalados, por lo que lo que te encuentras al abrir RStudio es cómo  la figura de  la derecha. Cuatro ventanas. Pinchando en los símbolos de pestaña de la parte superior derecha de cada ventana, puedes maximizar o minimizar cada una. Por defecto las de la izquierda son arriba el editor de texto y abajo la consola, y las de la izquierda las variables (environment) y el historial (history) arriba y todo lo demás,  incluyendo plots y ayuda, abajo. Aunque la colocación de las ventanas se puede modificar al gusto de cada cual.  Las más importantes son:

EDITOR: donde se guarda y edita el código. Aquí se abren los scripts, se modifican y se prueban. No olvides dar a guardar de vez en cuando.

CONSOLA: donde se corre el código. Para correr una línea de código del editor en la consola, se presiona la tecla “Run” (parte superior del editor).

Una vez dicho esto, lo mejor para empezar el leer y practicar. Para ello se pueden seguir:

Tutoriales

Aconsejo echarle un vistazo a ambos.

Webs y libros interesantes

Aunque hay muchos libros interesantes, para aconsejar me centraré en los materiales gratuitos.

Para profundizar con R y/o estadística:

Para las nociones básicas sobre código reproducible:

Algunas webs interesantes donde buscar soluciones a problemas  concretos o incluso preguntar si no están ya publicadas:

Guis

Otras guis interesantes para utilizar con R, aunque no las he probado mucho son:

 

 

 

My R lecture notes / Mis apuntes de R (updated 18-09-23)

I link you here the “R for life sciences” lecture notes I am doing, so they are easily available.

(If you prefeer in Spanish / Si prefieres en español)

Lecture notes Web html pdf epub
 1. Easy R start  html  pdf  epub
 2. Operations  html  pdf  epub
 3. Basic graphics and data  html  pdf  epub
 4. Hypothesis testing  html  pdf  epub

 

Some other interesting books and webs for R and statistics are:

For reproducible code basics:

To look for a specific topic:

Automatizar el paso de GPX a SHP desde R

Desde la lógica más pura del trabajo con datos en R, lo mejor es tener los datos originales sin modificar, ya que siempre que se modifican es una posibilidad nueva de error. Para poder trabajar con esos datos originales lo que hacemos es convertirlos en una matriz de datos con los que podamos trabajar.

En nuestro caso tenemos en una carpeta los datos GPX tal cual salen del GPS y una tabla CSV con el mismo nombre que contiene los datos recogidos en la libreta (Especie, nombre parcela, otros datos) sobre los mismos puntos que contiene el GPX. Es importante que el número de registros en uno y en otro sean los mismos y que estén en el mismo orden, aunque no es necesario que se llamen igual, es interesante hacerlo para poder comprobar que se han unido bien.

Pasos a seguir:

  1. Instalar rgdal
  2. Elegir directorios de trabajo
  3. Importar datos GPX en R
  4. Importar datos CSV en R
  5. Juntar todos los datos en una única SpatialPointsDataFrame
  6. Corregir errores de transcripción
  7. Eliminar repetidos
  8. Exportar SHP

Instalar rgdal

Para importar datos GPX en R necesitaremos instalar la librería rgdal. Al intentar instalarla puede ocurrir que falte gdal-config, por lo que habrá que instalarlo primero en la consola de Linux:

sudo apt-get install libgdal1-dev libproj-dev

y en la consola de R:

install.packages("rgdal", dep=T)         #rgdal package
library(rgdal)

Una vez instalada, sólo es necesario correr la última línea de código para volver a activar la librería las siguientes veces.

Elegir directorios de trabajo

Lo primero que vamos a hacer es fijar los directorios de trabajo:

#   Directorio de entrada de datos
DirDatosBrutos <- "Directorio donde están los datos GPX y CSV/"
#   Directorios para exportar datos
DirDatosSHPwp  <- "Directorio donde se querran guardar los archivos shp de puntos/"
DirDatosSHPtrk <- "Directorio donde se querrán guardar los shp de tracks/"

Observar que al principio no hay barra y al final hay una barra inclinada (en Mac o en Windows puede ser diferente). Para comprobar que están bien elegidos, podemos utilizar:

list.files()
list.files(DirDatosBrutos)
list.files(DirDatosBrutos, ".csv")
list.files(DirDatosSHPwp, ".shp")
list.files(DirDatosSHPtrk, ".shp")

Importar datos GPX en R

Para importar los datos GPX y convertirlos en variables espaciales de puntos (waypoints) o lineas (tracks) utilizaremos readOGR()

#=============================================================================== 
#          LEER DATOS GPX y csv y crear listas
#=============================================================================== 
GPXinNombres <- list.files(DirDatosBrutos, ".GPX")  # Lista nombres gpx importar
CSVinNombres <- list.files(DirDatosBrutos, ".csv")  # Lista nombres csv importar
#
#                         GPXwp           Lista de SpatialPointsDataFrame PUNTOS
GPXwp <- list()
for (i in GPXinNombres){
    nombre <- paste(DirDatosBrutos,i,sep="")
    GPXwp[[i]] <- readOGR(dsn = nombre, layer="waypoints")
}
#
#                       GPXtrk             Lista de SpatialLinesDataFrame TRACKS
GPXtrk <- list()
for (i in GPXinNombres){
  nombre <- paste(DirDatosBrutos,i,sep="")
  GPXtrk[[i]] <- readOGR(dsn = nombre, layer="tracks")
}
#                       COMPROBAR
for (i in 1:length(GPXwp)){
   plot(GPXwp[[i]])
   plot(GPXtrk[[i]],add=T)
}
# Debería salir un dibujo con los traks y waypoints de cada archivo gpx

Importar datos CSV en R

Lo ideal es que tengan el mismo nombre que los correspondientes GPX, pero si no, al menos deberían estar en el mismo orden alfabético. También los nombres de las variables deberían ser en todas las mismas. En el ejemplo son ID, Tipo, Lugar Parcela y Fecha, pero habrá que poner en cada caso los datos tomados en la libreta.

#                         CSVwp     Lista de matrices de datos espaciales
CSVwp <- list()
for (i in CSVinNombres){
nombre <- paste(DirDatosBrutos,i,sep="")
CSVwp[[i]] <- read.table(nombre,header=T,sep=",")      #,sep="\t")     #   Leer tabla
names(CSVwp[[i]]) <- c("ID","Tipo","Lugar","Parcela","Fecha")          #   Arreglar nombres variables
}
#

Limpiar variables del GPX

Podemos ver todas las variables introducidas en cada una de las tablas utilizando str() o head(). En el ejemplo sólo nos quedaremos con las variables 5, 1 y 2, “nombre”, “ele” y “time”, respectivamente, eliminando todas las demás ya que están vacías o no tienen información relevante. Comprobar siempre antes, ya que pueden variar en función del modelo de gps. Para eliminarlas de toda la lista de tablas a la vez, utilizaremos for().

#===============================================================================
#          LIMPIAR GPXs
#===============================================================================
#          Seleccionar las variables útiles y crear variable VarUtiles
head(GPXwp[[1]])
VarUtiles <- c(5,1,2)                           # Puede cambiar según modelo GPS
for(i in 1:length(GPXwp)){
GPXwp[[i]] <- GPXwp[[i]][,VarUtiles]
}
head(GPXwp[[1]])                                             # Comprobar cambios
#

Juntar las variables de los GPX con las de los CSV

Un detalle importante en este paso es que los nombres de las variables en el CSV no deben coincidir con las variables que se hayan conservado del GPX. Si es así debería cambiarse el nombre de una de las dos. En este momento tenemos una lista de gpx con tres variables cada uno y una lista de csv con 4 variables cada uno. Podrían ser más y no pasaría nada mientras estén en el mismo orden y lógicamente el número de puntos en los gpx coincida con el número de lineas de datos de cada csv. Si es así, esto debería funcionar.

#===============================================================================
#          JUNTAR GPX Y CSV
#===============================================================================
for(i in 1:length(GPXwp)){
for (j in names(CSVwp[[i]])){
GPXwp[[i]][[j]] <- CSVwp[[i]][[j]]
}
}
#                       COMPROBAR si se han añadido
head(GPXwp[[1]])

Juntar todos los datos en una única SpatialPointsDataFrame

para juntar todos los datos en una única tabla de datos espaciales “SpatialPointsDataFrame” creamos un objeto introduciremos la primera tabla y luego con for() vamos añadiendo las demás utilizando rbind() para añadir cada vez más filas. Más tarde quitaremos los elementos duplicados por lo que no debemos  preocuparnos si los archivos solapan o si se tomaron puntos repetidos.

WP <- GPXwp[[1]] # Crea tabla GPXWP
for(i in 2:length(GPXwp)){
WP <- rbind(WP,GPXwp[[i]]) # Une resto de tablas a GPXWP
}
# COMPROBAR
WP

Corregir errores de transcripción

Además de comprobar que se han introducido todos los casos es interesante comprobar también si existen errores de transcripción que haya que corregir. Para las variables cuantitativas suele ser interesante ver los valores máximos y mínimos. Las cualitativas se pueden comprobar viendo los niveles de cada factor.

levels (WP$Tipo)                                 # Debería haber sólo 4: A,H,R,S

En nuestro ejemplo resulta que en la variable Tipo se han introducido en algunos casos “Ref” en vez de “R” y “X”  en vez de “H”. Para corregirlo

levels(WP$Tipo)[levels(WP$Tipo)=="REF"] <- "R"   # Cambia REF por R
levels(WP$Tipo)[levels(WP$Tipo)=="X"] <- "H"     # Cambia X por H
levels (WP$Tipo)                                 # Debería haber sólo 4: A,H,R,S

Eliminar repetidos

Los puntos que coinciden en coordenadas y en Tipo se consideran puntos repetidos y por tanto serán eliminados para que en cada coordenada sólo pueda haber un punto de cada tipo.

#===============================================================================
#          ELIMINAR REPETIDOS
#===============================================================================
coordTipo <- cbind(coordinates(WP),WP$Tipo)           # Variables para repetidos
dupl <- duplicated (coordTipo)                        # casos duplicados
#                       COMPROBAR
# length(subset(dupl,dupl==TRUE))                     # Numero casos duplicados
# length(subset(dupl,dupl==FALSE))                    # Numero casos sin duplicar
#
str(WP)
WP <- WP[!dupl,]                                      # ELIMINAR CASOS REPETIDOS
str(WP)
#

Exportar SHP

Ahora ya podemos exportar un shp con todos los puntos de todos los gpx y con los datos del csv incluidos en la tabla de atributos utilizando la función writeOGR().

writeOGR(WP, dsn = DirDatosSHPwp, "WP-Definitivo", driver = "ESRI Shapefile", overwrite_layer = T)  # Crea WP.shp DEFINITIVO

También se puede exportar en otros formatos (GPX, CSV, etc.). Si además queremos exportar una capa de shp para cada uno de los archivos de traks de cada gpx, se puede hacer así

for(i in 1:length(GPXtrk)){
writeOGR(GPXtrk[[i]],DirDatosSHPtrk,names(GPXtrk)[i],driver="ESRI Shapefile",overwrite_layer = T)
}

No he conseguido (tampoco lo he necesitado por ahora) juntar todas las “SpatialLinesDataFrame” de las tracks de los gpx en una sólo ni eliminar los segmentos repetidos. Si en algún momento veo cómo hacerlo lo añadiría aquí. Si sabes cómo hacerlo puedes comentarlo.