El término “agricultura de precisión” se está haciendo cada día más común para referirse a la creciente utilización de las nuevas Tecnologías de la Información y la Comunicación (TICs) en la gestión de regadíos, viñedos, plantaciones de frutales, etc. Estas prácticas innovadoras están basadas en el monitoreo continuo del estado de los cultivos mediante dispositivos GPS, imágenes de satélites/drones y otro tipo de aparatos (sensores de humedad o temperatura), lo cual permite reaccionar a tiempo en caso de detectarse problemas. En este artículo vamos a exponer nuestros primeros resultados aplicando estas tecnologías para conseguir lo que podría denominarse una “ganadería de precisión” en las dehesas de porcino.

La dehesa es un agroecosistema ibérico de elevado valor ambiental y, como tal, se encuentra incluido en la Directiva Hábitats de la Unión Europea. La superficie que ocupa (3,1 millones de hectáreas) lo convierte en uno de los sistemas agroforestales más importantes de Europa (Papanastasis 2004). La dehesa tiene además una gran importancia socioeconómica, sobre todo en el cuadrante sudoccidental de la Península. Se trata de un sistema multi-funcional en el que se pueden desarrollar diferentes tipos de actividades agropecuarias, destacando la ganadería extensiva (Moreno y Pulido 2009) y, especialmente, los productos de alta calidad derivados del cerdo ibérico.

La clave de la calidad de los productos del cerdo ibérico es el régimen extensivo en el que se mantiene a estos animales y su alimentación, en la que juegan un papel fundamental las bellotas de encina Quercus ilex y alcornoque Quercus suber. Por ello, la baja disponibilidad de bellotas y su altísima variabilidad espacio-temporal puede limitar el potencial productivo de las fincas ganaderas. Patógenos como el oomiceto Phythophthora cinamomi, causante de la patología conocida por el nombre de “seca” (Corcobado et al. 2013), e insectos plaga que se alimentan de sus hojas (Canelo et al. 2018) provocan el decaimiento del arbolado reduciendo la disponibilidad general de bellotas. Aparte, la producción está sujeta a variabilidad dentro de las fincas por la propia heterogeneidad de las condiciones del terreno (Pérez-Ramos et al. 2014) y a las fluctuaciones inter anuales al ser las quercíneas árboles veceros (Espelta et al. 2008). A esta variabilidad en la producción de bellotas habría que añadir que generalmente no se tiene un conocimiento fino de los patrones de movimiento del ganado ya que, a diferencia de tiempos anteriores, el número de empleados por finca no permite su seguimiento constante.

Esta ausencia de información reduce la rentabilidad de las fincas y genera incertidumbres en la planificación de las explotaciones. La falta de un conocimiento en tiempo real de la disponibilidad de alimento y de los movimientos del ganado impide que se pueda responder regulando la carga ganadera, añadiendo alimento suplementario o modificando la distribución espacial del ganado en las dehesas. Para comenzar a solucionarlo, en los 3 últimos años han desarrollado un proyecto financiado por el Ministerio de Economía y Competitividad que incluía la aplicación de nuevas tecnologías a la actividad ganadera en las dehesas: concretamente las imágenes de drones y satélites para evaluar el vigor del arbolado y dispositivos GPS para el monitoreo en tiempo real del movimiento del ganado.

Las imágenes multiespectrales y la teledetección se han utilizado desde hace años para obtener información esencial del estado de cubiertas vegetales de diferentes características, permitiendo la estimación espacial y temporal de numerosos parámetros biofísicos y bioquímicos de la vegetación (Chuvieco et al. 2004, Yebra et al. 2008, Glenn et al. 2011, Melendo-Vega et al 2017). Gracias a estas imágenes se han podido calcular un gran número de índices de vegetación, entre los que destaca por su importancia el Índice de Vegetación de Diferencia Normalizado o NDVI (Rouse et al. 1974). El NDVI está fuertemente relacionado con la actividad fotosintética de las plantas, derivada de reflejar la radiación infrarroja en la estructura interna del mesófilo de las hojas verdes sanas y absorber la radiación visible roja por la clorofila de la hoja (Sellers et al. 1992). Aquellas zonas con valores altos de NDVI serían aquellas con mayor verdor o vigor vegetativo, es decir, donde la productividad primaria es mayor (Tucker 1979, Chuvieco et al 2004, Yebra et al. 2012).

El NDVI se ha utilizado para evaluar el estado y características de la vegetación en superficies agrícolas y forestales, principalmente en cubiertas vegetales densas de plantas herbáceas o arbustivas (Glenn et al. 2011, Yebra et al. 2012, Melendo-Vega et al. 2017) y en menor medida sobre superficies arboladas (Carreiras et al. 2006, Stagakis et al. 2012). No obstante, el empleo de técnicas de teledetección a gran escala (imágenes de satélite) se encuentra condicionado por una serie de limitaciones cuando existe una gran variabilidad espacial, espectral y temporal. Estas limitaciones se acentúan en ecosistemas altamente heterogéneos como la dehesa, que combina diferentes estratos vegetativos (árboles, matorral y pasto) mezclados con varios elementos inertes (suelo, agua, etc.). La utilización de imágenes multi-espectrales con mayor resolución en dehesas para conocer el vigor del arbolado (tomadas por drones) es pues un asunto pendiente, ya que previamente sólo se han utilizado para estimar la densidad de arbolado (Carreiras et al. 2006) y siempre han sido imágenes de satélites.

En cuanto a la utilización de dispositivos de seguimiento del ganado con GPS, diferentes estudios han demostrado su viabilidad y la utilidad de la información que ofrecen sobre todo cuando los animales se mantienen en régimen extensivo (Augustine y Derner 2013). En el caso del cerdo ibérico en la dehesa, existe un estudio previo en el que se siguió a 2 ejemplares durante el periodo de montanera con este tipo de dispositivos (Aparicio et al. 2006). Se obtuvo información básica sobre distancias recorridas, pero no se relacionó el movimiento de los animales con la productividad vegetal estimada por imágenes multiespectrales.

El objetivo general del estudio es, teniendo en cuenta las peculiaridades de la dehesa, analizar la efectividad de las imágenes multiespectrales tomadas por drones y satélites para estimar la productividad del arbolado, así como de los dispositivos GPS para conocer el movimiento de los cerdos ibéricos. La integración de ambas fuentes de información permite conocer los factores que determinan el uso del espacio de los animales y optimizar la gestión ganadera. Los objetivos concretos fueron:

  • Establecer, sobre imágenes tomadas por drones, un método automático de digitalización de copas de encina a partir de su marca espectral y averiguar cuál es la época más adecuada para hacerlo.
  • Analizar la relación entre el índice NDVI de los árboles y medidas de campo de vigor y producción de bellotas.
  • Conocer si existe relación entre el índice NDVI de los árboles y el de las herbáceas circundantes y, si es así, cómo varia estacionalmente para poder subir de escala y utilizar imágenes de satélite (menor resolución pero mayor cobertura espacial).
  • Estimar parámetros básicos del movimiento del ganado porcino (distancia recorrida por hora y por día, velocidad, etc.) a partir de los dispositivos GPS.
  • Analizar la selección de hábitat que realiza el cerdo durante la montanera en función de la estructura del hábitat y la productividad vegetal estimada a partir de imágenes de teledetección.

MATERIAL Y MÉTODOS

Área de estudio

El estudio se realizó en una dehesa ubicada en el norte de Cáceres (Extremadura), concretamente en la Finca Casablanca (término municipal de Guijo de Granadilla, coordenadas 40°09’00.5″N 6°07’03.5″O). El clima en la zona es típicamente mediterráneo, con una temperatura media anual de 16.3ºC y una marcada sequía estival en la que pueden alcanzarse los 40ºC. Las precipitaciones se concentran durante el invierno y escasean en verano, la media anual es de 700 mm.

Los datos se tomaron a lo largo de 2017 en un sector de la finca de 96ha, en el que desde mediados de octubre hasta enero del año siguiente permanecieron 31 cerdos ibéricos durante la montanera en régimen extensivo (0,32 animales/ha). En esta zona de la dehesa, la mayor parte de la superficie está ocupada por herbazal con encinas Quercus ilex dispersas (densidad media 20 árboles/ha), pero hay también partes de pradera desarbolada (Fig. 1).

Muestreo

Vuelo de Dron (RPAS)

Para la obtención de las imágenes multiespectrales de alta resolución se utilizó un sensor con 4 bandas espectrales (rojo, verde, infrarrojo y borde rojo) y una resolución (tamaño de píxel) de 10cm, aeroportado por un dron (RPAS) de tipo ala fija. Las imágenes se tomaron en la primavera (25 mayo) y otoño (20 octubre) en condiciones de cielo descubierto y ausencia de viento. Se eligieron estas fechas por presentar el mayor contraste fenológico en el ecosistema, que se debería traducir en diferencias en el contraste espectral entre pasto y arbolado. En primavera, tanto herbáceas como árboles están verdes y con intensa actividad fotosintética, por el contrario, al final del verano principio de otoño el pasto está seco (Mendiguren et al. 2015, Migliavacca et al. 2017; Fig. 2).

Cámara hemisférica

De forma simultánea a la adquisición de las imágenes multiespectrales, y con objeto de obtener datos de campo sobre alguna variable biofísica que nos permita validar la metodología propuesta, se tomaron fotografías con una cámara hemisférica de las copas de 92 árboles seleccionados al azar por toda la superficie de la finca. Esta metodología permite calcular de forma indirecta el LAI (Leaf Area Index) de la cubierta, que hace referencia a la superficie foliar fotosintéticamente activa (contabilizada por una sola cara) por unidad de suelo (Watson 1947). Las fotografías se tomaron con una cámara digital Nikon Coolpix 995, a la cual se acopló una lente hemisférica “fisheye” Nikon FC-E8. Las imágenes se tomaron en agosto de 2017, a 30 cm del suelo, sobre un trípode, en posición horizontal y con la cámara orientada hacia el norte (Fig. 3).

Medidas de cosechas de bellotas

En el año anterior al estudio se eligieron aleatoriamente 32 árboles en los que se midió la cosecha de bellotas de 2017. Para ello, en todas estas encinas se colgaron trampas de semillas de la parte inferior de sus ramas. Estas trampas fueron recipientes de plástico de 0,19m2 de superficie que se colocaron antes del verano, y en las cuales fueron recogiéndose todas las bellotas que cayeron desde entonces hasta el final de la temporada. El número de trampas de semillas osciló de 4 a 8 y fue proporcional a la superficie de la copa del árbol (que se calculó a partir de ortofotografía con un Sistema de Información Geográfica). Las bellotas se recogieron cada 15 días desde principios de septiembre hasta que la caída cesó a finales de Diciembre. El contenido de las trampas de semillas se llevó al laboratorio, se contó el número total de bellotas crecidas y, a partir de éste, se calculó la producción de bellotas por metro cuadrado de cada árbol. Multiplicando ésta producción por la superficie de la copa obtuvimos la estimación del número total de bellotas producido por cada encina. Esta metodología para estimar cosechas a partir de muestras de trampas de semillas ha demostrado ser muy fiable en muchos estudios previos en diferentes encinares y dehesas (Bonal et al. 2007).

Seguimiento del movimiento del ganado por medio de dispositivos GPS

Para el seguimiento del movimiento del ganado se le colocó a un cerdo ibérico durante 11 días (del 15 al 25 de diciembre de 2017, ambos inclusive) un emisor de señales GPS fabricado por la empresa Microsensory. El emisor tenía forma de petaca (peso 500gr) y se le acopló al animal por medio de un arnés fabricado a tal efecto; contenía una batería, un emisor GPS y una tarjeta GPRS. El dispositivo estaba programado para encenderse cada 15 minutos, buscar la señal del satélite y calcular la ubicación del animal y, por último, cada 24 horas, mandar una señal de todas las ubicaciones al satélite gracias a la tarjeta GPRS. Todas esas ubicaciones (longitud y latitud) quedaban registradas en el servidor de la empresa (Micorsensory) al cual accedíamos para descargarlas diariamente.

Procesamiento de imágenes

Imágenes multiespectrales

Las imágenes adquiridas por el sensor multiespectral aeroportado (RPAS) y que contienen los datos de reflectividad, fueron corregidas radiométricamente y geométricamente por la empresa encargada de realizar los vuelos. Para la visualización y tratamiento de imágenes se utilizó la interfaz ArcMap del Software ArcGis 10.2.2., utilizando el mismo sistema de proyección de coordenadas original (WGS84) para las imágenes y para los datos GPS obtenidos en campo.

Para obtener el valor medio de NDVI de la copa de cada árbol de estudio, se realizó una clasificación no supervisada de 3 y 4 clases en función de sus valores de NDVI (métodos de digitalización automática 3 y 4), con el objeto de separar en una nueva capa ráster los elementos (suelo, agua, pasto seco y árboles) con respuestas espectrales diferentes. Esta clasificación se hizo en otoño, ya que es cuando el pasto está seco y su respuesta espectral puede diferenciarse de la de los árboles (Fig. 2). A continuación, depuramos la capa ráster hasta quedarnos únicamente con la última clase generada, que se corresponde con las copas de todos los árboles de la finca. Por último, obtuvimos un shapefile de polígonos con todas las copas de nuestros árboles de estudio. A partir de éstas copas digitalizadas pudimos extraer los valores medios de NDVI de los píxeles correspondientes a nuestros árboles de estudio.

Para analizar la exactitud de este método de clasificación seleccionamos una submuestra aleatoria de 24 árboles. En estas encinas se realizaron 2 tipos de digitalizaciones manuales de las copas en las imágenes de otoño, una digitalización bruta (polígono más o menos circular sobre la copa) y otra digitalización fina (ajustando esos polígonos a la copa real). Se extrajeron los valores medios de NDVI de cada una de las copas y se compararon con los obtenidos a partir de las dos clasificaciones supervisadas automáticas.

Con el objetivo de determinar la época de mayor contraste entre pasto y árbol, se realizó un buffer de 2m alrededor de la copa de cada árbol en primavera y otoño, obteniendo los valores medios de NDVI de todos los píxeles contenidos en el buffer. Esta información permite por un lado confirmar la época óptima de delimitación de copas y, por otro, saber si en algún periodo del año la relación entre el NDVI de la copa y el pasto es tan alta que posibilitaría utilizar imágenes con menor resolución (imágenes de satélite) (Fig. 4) para conocer el vigor del arbolado de una dehesa.

Para relacionar el uso del hábitat por parte de los cerdos en función de la productividad vegetal (estimada según el índice NDVI) se utilizaron imágenes del Satélite Sentinel con una resolución (tamaño de pixel) de 10m de lado. La razón de tomar la imagen de satélite es que el seguimiento del movimiento de los cerdos se hizo en plena montanera (diciembre), por lo que no había imagen reciente del dron pero sí del satélite (fecha de la imagen 18/12/2017).

Fotografías hemisféricas

En cada uno de los 92 árboles se tomaron 3 fotografías de la copa árbol en diferentes orientaciones, calculando una vez procesadas la media de sus resultados (Índice LAI). Para el análisis de las fotografías hemisféricas se utilizó el software especializado ImageJ (Abràmoff et al. 2004), que nos permite realizar una binarización de las imágenes (Jonckheere et al. 2004) y establecer un umbral de selección automático que clasifique los pixeles en cielo y vegetación, para obtener a continuación la “gap fraction” o frecuencia de huecos en la cubierta (Weiss et al. 2004), relacionada a su vez con el LAI (Índice de Área Foliar).

RESULTADOS Y DISCUSIÓN

Métodos de digitalización

Las imágenes tomadas por los drones permiten individualizar con precisión la copa de cada árbol de la dehesa y conocer su valor de NDVI, aunque la eficiencia de la delimitación automática de las copas difiere dependiendo de la estación. Las imágenes de la mayoría de los satélites (Landsat, Sentinel) tienen resoluciones de más de 10m de tamaño de pixel, lo que impide poder individualizar copas. En bosques de dosel continuo y monoespecíficos esto no es problema (Fernández-Martinez et al. 2015), pero en la dehesa esa resolución mezcla en el mismo pixel pasto y árbol haciendo imposible distinguir la marca espectral de uno y otro. Esta marca espectral de hecho difiere de forma marcada en otoño, cuando el pasto está seco y las encinas aún mantienen su actividad fotosintética (Fig. 2). Es en esta época cuando, a partir de la imagen NDVI del dron, es más efectivo hacer una rápida clasificación automática para delimitar las copas del pasto circundante. El resultado de la clasificación en 3 clases fue mejor que el de 4 clases. El modelo de regresión lineal robusto mostró que el NDVI medio de las copas de los árboles delimitados a partir de la clasificación automática estaba fuertemente relacionado con ese mismo valor calculado dentro de los polígonos de las copas delimitados manualmente (R2=0,95; F1,27=553; t=23; P<0,001).

Relación del índice de vegetación (NDVI) con el vigor del arbolado y la producción de bellotas

El índice NDVI se relacionó positivamente con el vigor del arbolado y la producción de bellotas. Los árboles con copas con valores más altos de NDVI alcanzaron valores más altos de Índice de Área Foliar (LAI). El modelo de regresión lineal fue altamente significativo (F1,90=124; t=11,16; P<0,001; Fig. 3). La relación positiva entre ambas variables ha sido documentada con anterioridad, aunque el valor de los coeficientes de regresión ha sido muy variable, yendo desde el 30% de la varianza explicada (Franklin et al. 1997) hasta el 84% (Melendo-Vega et al. 2017). En nuestro caso, teniendo en cuenta que hemos usado imágenes de gran resolución (10x10cm de píxel), la varianza explicada (58%) puede considerarse moderadamente alta, aunque podría haberlo sido más si el vuelo del dron se hubiese retrasado a agosto, que es cuando se tomaron las imágenes hemisféricas.

El Índice de Área Foliar es un indicador del vigor, ya que está relacionado con la cobertura de hojas de la copa y por tanto con la cantidad de luz que puede interceptar el árbol. Por esta razón se encuentra correlacionado con el NDVI, que mide la actividad fotosintética. El decaimiento y la defoliación provocan la pérdida de superficie de hojas lo cual, en último extremo, reduce la producción de bellotas (Canelo et al. 2018). De hecho, los árboles con mayores valores de NDVI en mayo fueron también los que produjeron un mayor número de bellotas (F1,30=6,87; t=2,62; P=0,01). El modelo de regresión explicó significativamente 18,6% de la varianza de la producción. La fracción no explicada puede deberse a factores accidentales que hacen que falle el cuajado de frutos en algunas encinas con potencial para engordarlos, por ejemplo, que su floración coincida con días de lluvia (Koenig et al. 2015). La relación entre producción de bellotas e índice de productividad vegetal ya se había documentado anteriormente, pero siempre en bosques de dosel continuo y con imágenes de menor resolución (tamaños de píxel de 30 a 4m) tomadas por satélites (Fernández-Martínez et al. 2015) o cámaras portadas por avionetas (Yao y Sakai 2010). Es la 1ª vez que se encuentra esta correlación a partir de imágenes de drones y en una sábana de quercíneas como la dehesa.

Relación entre el NDVI del árbol y el pasto circundante: época de digitalización y escalado a satélite

Las encinas más vigorosas (índice NDVI más alto) se encontraron en zonas de la finca donde el índice de vegetación del pasto circundante fue también más alto, aunque la relación árbol-pasto circundante varió con la fecha. En mayo, ambos (árbol y pasto), presentan valores de NDVI positivos y altos, y es cuando se da la relación positiva. Por el contrario, en el otoño, tras la sequía estival, esta relación desaparece ya que, pese a mantenerse la variabilidad de NDVI entre árboles, no la hay en el pasto al estar seco en toda la finca (Fig. 2). El modelo de ANCOVA muestra estas diferencias, ya que el factor momento del año (mayo/octubre) tiene un efecto significativo sobre la relación entre la covariable (NDVI del pasto circundante) y la variable dependiente (NDVI de la copa del árbol) (estimate=-1,13±0,25 SE; t=-4,51; P<0,001). Algunos estudios previos habían tratado el contraste fenológico entre árbol y pasto (Mendiguren et al. 2015; Migliavacca et al. 2017), pero su relación nunca había sido estimada por medio de imágenes tomadas por drones. Además, en este caso queremos destacar sus implicaciones prácticas ya que, si bien el verano es el periodo ideal para la digitalización automática de copas de árboles, la primavera lo sería para la utilización de imágenes de satélite, ya que es cuando ambos estratos presentan una respuesta espectral más parecida.

Patrones de movimiento del cerdo ibérico en montanera medidos con GPS y factores que determinan su uso del hábitat

El cerdo equipado con el GPS se movió fundamentalmente por el día (considerando como horas de luz en Diciembre desde las 7 a las 17:00h solar), por la noche permaneció normalmente reposando dentro del cobertizo que sirve de refugio a la piara (Fig. 1). No obstante, al estar el cobertizo abierto y al igual que se ha visto anteriormente (Aparicio et al. 2006), se registró algún movimiento durante la madrugada, aunque con menos frecuencia (Fig. 6). En ambos casos, la distancia recorrida por 15 minutos más frecuente fue menor a 100m. Lo que corresponde a una velocidad de 0,4km/h. También se registraron desplazamientos de hasta 700m en 15 minutos, lo cual correspondería a una velocidad máxima de 4,5km/h. Estas últimas velocidades se relacionan con movimientos direccionales dentro de diferentes zonas de la dehesa, mientras que los más frecuentes serían los desplazamientos más lentos de alimentación. A partir de los desplazamientos por 15 minutos se calculó que la distancia media recorrida por día durante las horas de luz fue de 6,51km (desviación estándar ±2.19). Este valor es algo superior al registrado por el único estudio previo disponible, en el que la media se situó ligeramente por debajo de los 5km diarios (Aparicio et al. 2006).

El cerdo no utilizó toda la finca de forma homogénea, su movimiento estuvo condicionado por la presencia de refugios, la estructura de la vegetación y su productividad. Muchas de las localizaciones se dieron cerca del cobertizo que les sirve de refugio, al igual que se observa otra pequeña agregación cerca de una pequeña charca (Fig. 1). Asimismo, la presencia de encinas es un factor determinante, ya que en la zona de pradera apenas se registraron localizaciones en comparación a la gran superficie que ocupa (Chi2=87,01; d.f.=1; p<0,001); Fig. 1). Es sabido que la hierba es una fuente de alimento importante durante la montanera (García et al. 2009), pero los resultados muestran que la consumen en zonas en las que también hay bellotas. Los cerdos seleccionaron aquellas en las que la productividad vegetal fue mayor, ya que el NDVI medio correspondiente a los puntos donde se localizó al cerdo fue mayor que el del mismo número de puntos distribuidos al azar por la finca (t-test: t=-5,38; N=286; df=548; p<0,001). El uso del hábitat del cerdo coincide pues con el de especies de fauna silvestre (Wiegand et al. 2008) pero esta no sería la razón por la que no utiliza la zona de pradera ya que, pese a no haber árboles, en algunas zonas los valores de NDVI son altos.

CONCLUSIONES

La utilización de imágenes multiespectrales de precisión tomadas por drones permite estimar con antelación (en primavera) cómo será la cosecha de bellotas de ese año, pudiendo anticipar capacidades de carga de las fincas, necesidades de alimentación suplementaria, etc. La relación entre productividad de las encinas y el pasto circundante muestra, además, que se pueden utilizar imágenes de satélite en determinadas épocas del año, lo que facilita una mayor cobertura espacial de las fincas. El seguimiento del movimiento de los cerdos corrobora la importancia de mantener una buena densidad de encinas en las fincas, ya que las zonas arboladas son seleccionadas positivamente. Más novedoso es el resultado que muestra la importancia que tiene la productividad vegetal en el uso del espacio, y que significa que sería bueno tener esta información con antelación para decidir las zonas prioritarias donde introducir los cerdos en montanera.

De cara al futuro, nuestro objetivo es continuar esta línea de investigación aumentando las zonas y las series temporales para corroborar que es una metodología generalizable. En relación a las imágenes, queremos probar con otros índices aparte del NDVI para saber si puede que predigan aún mejor las variables de interés. En cuanto al estudio del movimiento, pretendemos utilizar nuevos dispositivos de menor tamaño (como mini-emisores GPS en crotales) que reduzcan las faltas en la frecuencia de emisión. Asimismo, el objetivo es marcar un mayor número de cerdos y abordar cuestiones relacionadas con la productividad, por ejemplo, la relación entre la distancia diaria recorrida y la calidad de la carne. En resumen, los resultados de este estudio muestran la eficacia que puede tener la aplicación de las nuevas Tecnologías de la Información y la Comunicación para conseguir una verdadera “ganadería de precisión” en las dehesas de porcino.

Autores:

  • Carlos Pérez-Izquierdo, Tara Canelo, Raúl Bonal. Grupo de Investigación Forestal, INDEHESA (Instituto de Investigación de la Dehesa). Escuela de Ingeniería Forestal y del Medio Natural. Universidad de Extremadura.
  • Álvaro Gaytán. Department of Environment, Ecology and Plants, Universidad de Estocolmo, Suecia.

Deja un comentario