Modelos de distribución de especies aplicados a estudios de flora amenazada (guión de prácticas)

Blas Benito de Pando Unidad de Conservación Vegetal Departamento de Botánica Universidad de Granada

Nota del autor Este trabajo contiene una guía práctica muy completa (y compleja) para modelar la distribución de especies. Trata todo el proceso de modelado, desde la confección de las variables ambientales hasta la preparación de los resultados para publicación. El objetivo de este escrito es proporcionar al no iniciado los conocimientos teóricos y prácticos mínimos, así como mostrar las herramientas necesarias para trabajar con garantías en este área. Todo el material contenido se basa en la bibliografía existente, pero usted debe entender que no es una transcripción literal, sino una interpretación de ese material bibliográfico, basada en la experiencia del autor en el campo de los modelos de distribución. Teniendo en cuenta que es un campo en desarrollo activo, será fácil encontrar investigadores y literatura que no estén en acuerdo con el modo en que algunas cuestiones relativas a los modelos de distribución. Tanto en caso de desacuerdo, como en caso de error, le agradecería que me lo comunicara para su revisión o corrección a alguno de los siguientes correos: [email protected] [email protected] Teniendo en cuenta que todo el proceso expuesto se basa en herramientas de software gratuito y libre, este trabajo tiene una licencia libre Creative Commons.

Reconocimiento – Sin obra derivada – No comercial Según los términos de esta licencia, usted puede distribuir y copiar este trabajo, siempre que otorgue al autor el crédito que le corresponde. Usted no puede obtener beneficio comercial de los contenidos, ni puede realizar obras derivadas de los mismos.

Índice SECCIÓN 1: PREPARACIÓN DEL ENTORNO DE TRABAJO.............................................................1 1.1 Instalación de GRASS.................................................................................................................1 1.2 Instalación de Notepad++............................................................................................................1 1.3 Instalación de Open Office..........................................................................................................1 1.4 Instalación de Google Earth........................................................................................................2 1.5 Instalación de OpenModeller.......................................................................................................2 1.6 Instalación de Octave..................................................................................................................2 SECCIÓN 2: PREPARACIÓN DE LAS VARIABLES AMBIENTALES..................................................3 2.1 CREACIÓN DE LA BASE DE DATOS DE GRASS......................................................................3 2.1.1 Generar las localidades de GRASS.....................................................................................3 2.1.2 El interfaz de GRASS..........................................................................................................4 2.1.3 Preparación de la región de trabajo.....................................................................................5 2.2 PREPARACIÓN DEL MODELO DE ELEVACIONES...................................................................6 2.2.1 Descarga del modelo de elevaciones..................................................................................7 2.2.2 Importación del MDE la base de datos de GRASS.............................................................7 2.3 CREACIÓN DE UNA MÁSCARA.................................................................................................8 2.4 VARIABLES TOPOGRÁFICAS....................................................................................................9 2.4.1 Pendientes y orientaciones..................................................................................................9 2.4.2 Posición topográfica..........................................................................................................10 2.4.3 Índice topográfico de humedad..........................................................................................11 2.4.4 Radiación solar potencial...................................................................................................12 2.4.5 Diversidad topográfica.......................................................................................................13 2.5 VARIABLES DE TELEDETECCIÓN..........................................................................................14 2.5.1 Importación y procesamiento de bandas Landsat.............................................................14 2.5.2 NDVI..................................................................................................................................16 2.6 VARIABLES CLIMÁTICAS.........................................................................................................16 2.6.1 Precipitación......................................................................................................................16 2.6.2 Temperatura.......................................................................................................................18 2.6.2.1 Cálculo de variables predictoras................................................................................18 2.6.2.2 Cálculo de la ecuación de regresión..........................................................................19 2.7 LA CORRELACIÓN ENTRE VARIABLES..................................................................................22 2.8 ANÁLISIS DE COMPONENTES PRINCIPALES.......................................................................25 2.9 CAMBIO DE RESOLUCIÓN DE LAS VARIABLES....................................................................26 2.9.1 Introducción.......................................................................................................................26 2.9.2 Cambiando resoluciones...................................................................................................27 2.10 EXPORTACIÓN DE LAS VARIABLES AMBIENTALES...........................................................28 SECCIÓN 3: PREPARACIÓN DE LOS REGISTROS DE PRESENCIA..............................................29 3.1 Tipos y calidad de los datos.......................................................................................................29 3.2 Números mínimo y máximo de presencias................................................................................30

3.3 Densidad de los datos...............................................................................................................30 3.4 Partición aleatoria de datos.......................................................................................................31 3.5 Preparación de datos de presencia...........................................................................................31 SECCIÓN 4: MODELOS DE DISTRIBUCIÓN CON OPENMODELLER.............................................33 4.1 Preparación de las variables.....................................................................................................33 4.2 Ejecución de un experimento de prueba...................................................................................34 4.3 Lógica interna de los distintos algoritmos..................................................................................34 4.3.1 Bioclim...............................................................................................................................34 4.3.2 Climate Space Model.........................................................................................................35 4.3.3 Envelope Score.................................................................................................................36 4.3.4 Environmental Distance.....................................................................................................36 4.3.5 GARP................................................................................................................................37 4.3.6 Support Vector Machines...................................................................................................40 4.4 Algunas consideraciones sobre los resultados..........................................................................41 SECCIÓN 5: EVALUACIÓN DE MODELOS DE DISTRIBUCIÓN.......................................................42 5.1 Evaluación de modelos binarios (presencia/ausencia)..............................................................42 5.1.1 Cálculo de la sensibilidad, error de omisión y error de comisión.......................................42 Caso práctico: evaluación de modelos binarios midiendo la sensibilidad..............................44 5.1.2 Partición aleatoria de datos...............................................................................................44 Caso práctico: evaluación de modelos binarios utilizando partición aleatoria de datos.........45 5.1.3 registros de ausencia y matriz de confusión......................................................................49 Caso práctico: evaluación de modelos binarios utilizando partición aleatoria de datos y registros de ausencia....................................................................................................................51 5.1.4 registros aleatorios como sustitutos de las ausencias.......................................................53 Caso práctico: evaluación de modelos binarios utilizando partición aleatoria de datos, puntos aleatorios y técnicas de remuestreo......................................................................................55 5.2 Evaluación de modelos continuos.............................................................................................60 5.2.1 La curva ROC....................................................................................................................60 Caso práctico: Cálculo de la curva ROC para los modelos de Linaria nigricans...................63 SECCIÓN 6: TRANSFORMACIÓN DE MODELOS CONTINUOS EN BINARIOS..............................67 Caso práctico: transformación del mejor modelo continuo de Linaria nigricans en binario. . .69 SECCIÓN 7: APLICACIONES PRÁCTICAS DE LOS MODELOS DE DISTRIBUCIÓN.....................70 7.1 CARTOGRAFÍA DE POBLACIONES.........................................................................................70 7.1.1 Medición del área de ocupación........................................................................................70 7.1.2 Cartografía de detalle........................................................................................................71 7.1.3 Cartografía de reconocimiento..........................................................................................71 7.1.4 Comparación con una cartografía real...............................................................................72 7.1.5 Discusión...........................................................................................................................73 7.2 BÚSQUEDA DE NUEVAS POBLACIONES DE Linaria nigricans..............................................73 7.2.1 Procedimiento de búsqueda de nuevas poblaciones.........................................................74 7.2.2 Discusión...........................................................................................................................77 7.3 ENSAMBLADO DE BIODIVERSIDAD PARA ASISTIR EN EL DISEÑO DE RESERVAS DE

FLORA.............................................................................................................................................78 7.3.1 Preparación de los modelos de distribución......................................................................79 7.3.2 Discusión...........................................................................................................................83 7.4 EVALUACIÓN DEL IMPACTO POTENCIAL DEL CAMBIO CLIMÁTICO EN LA DISTRIBUCIÓN DE UNA ESPECIE AMENAZADA....................................................................................................84 7.4.1 Escenarios de Cambio Climático.......................................................................................85 7.4.2 Proyección de modelos de distribución con OpenModeller...............................................85 7.4.3 Análisis de los modelos ....................................................................................................86 7.4.3.1 Distribución potencial actual y futura.........................................................................86 7.4.3.2 Diferencial de idoneidad: frentes de avance y retroceso...........................................86 7.4.4 Discusión...........................................................................................................................87 7.5 ENSAMBLADO DE MODELOS PARA PROYECCIONES DE DISTRIBUCIÓN EN ESCENARIOS DE CAMBIO CLIMÁTICO......................................................................................................87 7.5.1 Escenarios climáticos........................................................................................................87 7.5.2 Proyección de modelos de distribución.............................................................................88 7.5.3 Análisis de los modelos entre escenarios..........................................................................88 7.4.3.2 Ensamblado de modelos...........................................................................................89 7.4.4 Discusión...........................................................................................................................89 7.6 EJERCICIOS PROPUESTOS...................................................................................................90 7.6.1 Cartografía de poblaciones................................................................................................90 7.6.2 Búsqueda de nuevas poblaciones.....................................................................................90 7.6.7 Explora libremente.............................................................................................................90 BIBLIOGRAFÍA....................................................................................................................................91

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

SECCIÓN 1: PREPARACIÓN DEL ENTORNO DE TRABAJO (Tiempo estimado: 10 minutos) Para trabajar cómodamente se ha preparado un directorio de trabajo para el curso. Puedes encontrarlo en C:, y se llama CURSO_MODELOS. Este directorio contiene los instaladores de los programas necesarios para seguir el curso (carpeta PROGRAMAS), una carpeta con el material geográfico necesario (GEODATOS), y las carpetas llamadas GUIONES, EVALUACION, y RESULTADOS, de las que hablaremos más tarde. En primer lugar vamos a preparar el entorno de trabajo para elaborar las variables ambientales que después utilizaremos para generar los modelos de distribución. Todos los programas que vamos a utilizar durante la práctica son gratuitos y de código abierto

1.1 Instalación de GRASS GRASS (http://grass.itc.it/) es un programa de procesamiento de datos geográficos libre y gratuito, muy potente, con el que vamos a hacer gran parte del trabajo. No se ha seleccionado como software para este curso por su condición de gratuito, sino porque es el más completo programa de procesamiento raster que puedes encontrar. Para descargar instaladores para distintas plataformas debes visitar: http://grass.itc.it/download/index.php. El instalador nativo para Windows está en: http://grass.itc.it/grass63/binary/mswindows/native/. En este momento no tienes que descargar nada, porque el instalador está en C:\CURSO_MODELOS\ PROGRAMAS\GRASS. Ejecuta la instalación por defecto. No es necesario que lo inicies aún, porque después lo veremos con detalle.

1.2 Instalación de Notepad++ Notepad++ (http://notepad-plus.sourceforge.net/es/site.htm) es un procesador de texto plano que nos va a resultar muy útil para determinadas operaciones. Puedes encontrar el instalador en C:\CURSO_MODELOS\PROGRAMAS\NOTEPAD_PLUS. El instalador de este programa puedes descargarlo de: http://sourceforge.net/project/showfiles.php?group_id=95717&package_id=102072.

1.3 Instalación de Open Office De esta suite ofimática (http://es.openoffice.org/) en realidad solo necesitamos el programa de hojas de cálculo Calc, por lo que durante la instalación, puedes descartar el resto de aplicaciones. Para las operaciones que nos ocupan lo preferimos frente a Microsoft Excel porque es más respetuoso con los formatos en la exportación a texto plano. Puedes encontrar el ejecutable en la carpeta OPENOFFICE del directorio del curso. Cuando termines la instalación, abre Calc, y en el menú Herramientas abre Corrección Automática, y desactiva todas las entradas de la pestaña Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 1

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Opciones; así evitaremos comportamientos molestos del programa más adelante.

1.4 Instalación de Google Earth. Sin duda ya conoces Google Earth (http://earth.google.es/); el instalador, Google Updater.exe está en la carpeta PROGRAMAS/GOOGLE_EARTH.

1.5 Instalación de OpenModeller OpenModeller (http://openmodeller.sourceforge.net/) es un programa diseñado para generar modelos de distribución de especies utilizando distintos algoritmos. Para instalar la última versión (1.0.7) de OpenModeller, necesitas actualizar la librería C++ (el lenguaje informático en el que está programado OpenModeller). La tienes disponible en CURSO_MODELOS/PROGRAMAS/OPENMODELLER/ vcredist_x86.exe: ejecuta el archivo para instalarla. Esta librería también puedes descargarla gratuitamente de: http://www.microsoft.com/downloads/thankyou.aspx?familyId=32bc1beea3f9-4c13-9c99-220b62a191ee&displayLang=en#). Una vez instalada la librería, puedes comenzar la instalación de OpenModeller ejecutando openModellerDesktopSetup1.0.7.exe., dejando las opciones por defecto. Los instaladores de OpenModeller en sus distintas versiones, disponibles en la carpeta mencionada anteriormente están disponibles en http://sourceforge.net/project/showfiles.php?group_id=101808&package_id=142057.

1.6 Instalación de Octave Octave (http://octave.sourceforge.net/) es un lenguaje de cálculo numérico, compatible con Matlab., con una sintaxis relativamente sencilla, que permite programar algoritmos rápidamente. En la carpeta CURSO_MODELOS/PROGRAMAS/OCTAVE está el instalador octave-3.0.1-setup.exe. Instala el programa teniendo en cuenta estas opciones: ●

en la ventana de selección de CPU, elige la opción Generic.

●

en la selección de gráficos, elige gnuplot.

Este instalador se descargó desde: http://sourceforge.net/project/showfiles.php? group_id=2888&package_id=40078. Si ya has instalado sin problemas estos programas, tienes el ordenador preparado para la primera fase de cualquier proceso de modelización, que es la creación de las variables ambientales necesarias.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 2

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

SECCIÓN 2: PREPARACIÓN DE LAS VARIABLES AMBIENTALES (Tiempo estimado: 2 horas) Preparar un conjunto apropiado de variables ambientales para hacer modelos de distribución es un proceso laborioso, que exige conocimientos relativamente avanzados en Sistemas de Información Geográfica, y cierta destreza manejando en cada caso el programa adecuado. Es cierto que en la red existen conjuntos de datos preparados para hacer modelos de distribución, como el de la página WorldClim (http://www.worldclim.org/), o el Atlas Climático Digital de la Península Ibérica (http://opengis.uab.es/wms/iberia/mms/index.htm) pero tienen limitaciones en su aplicación a estudios de pequeña escala. En general, para conseguir un buen conjunto de variables, será necesario buscar, procesar e integrar información espacial de distintas fuentes. A continuación vamos a desarrollar un protocolo de trabajo útil para crear un conjunto de variables ambientales de alta resolución espacial en el que se integrará información topográfica, climática y obtenida mediante teledetección.

2.1 CREACIÓN DE LA BASE DE DATOS DE GRASS GRASS va a ser nuestra herramienta principal para elaborar los mapas de variables. Para trabajar con GRASS es necesario definir: ●

Una base de datos en la que se guardarán las capas (en C:\CURSO_MODELOS\GEODATOS\GRASSDB).

●

Una localidad de trabajo (LOCATION), que define el área de trabajo y el sistema de coordenadas (NOTA: dentro de una misma localidad solo se puede utilizar un mismo sistema de coordenadas).

●

Un directorio de mapas (MAPSET), que guarda capas concretas, y comparte área de trabajo y sistema de coordenadas con la localidad.

●

Una REGION, que determina el área real de trabajo del programa, qué puede ocupar toda la extensión de la localidad, o una ventana concreta dentro de la misma, y la resolución de trabajo, que se puede ir cambiando sobre la marcha como veremos a continuación.

Para nuestro trabajo vamos a crear dos localidades definidas dentro de los límites geográficos de la provincia de Almería; una en sistema de coordenadas geográficas (latitud-longitud), con datum WGS 84 (código EPSG 4326), y otra con un sistema de coordenadas proyectadas UTM, con datum ED50 (código EPSG 23030). 2.1.1 Generar las localidades de GRASS Para generar ambas localidades, inicia GRASS. Posiblemente saldrá un mensaje de error (porque aún no hay una base de datos definida); pulsa ok para continuar, y sigue estos pasos. ●

GIS Data Directory = C:/CURSO_MODELOS/GEODATOS/GRASSDB

●

Define new location with = EPSG codes

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 3

Modelos de distribución y estudios de flora amenazada ○

Name of new location = ALMERIA_latlong

○

EPSG code number of projection = 4326

○

Pulsa el botón Define location

Guión de prácticas

Ya hemos creado la base de datos en coordenadas geográficas; ahora vamos a generar la base de datos en UTM, repitiendo algunos de los pasos anteriores. ●

Define new location with = EPSG codes ○

Name of new location = ALMERIA_utm

○

EPSG code number of projection = 32630 (utm wgs84 zona 30N)

○

Define location

○

Select datum transformation parameters = marcar la opción 5; Used in Spain (except Northwest).

○

Ok.

2.1.2 El interfaz de GRASS Ahora que tenemos ambas bases de datos configuradas, habrás observado que en ambas se crea automáticamente un Mapset denominado PERMANENT, en el que se almacenarán todos los datos geográficos. Para entrar en GRASS, en Project Location selecciona ALMERIA_latlong, y en Accesible Mapsets selecciona PERMANENT. Pulsa Enter GRASS para iniciar el programa. Habrás observado que el interfaz de GRASS tiene varias ventanas (ver Figura 1): ●

GRASS 6.3.0 GIS Manager: ○

barra de menús de esta ventana accedemos a la mayor parte de las funcionalidades de procesamiento de GRASS: menús Raster, Imagery y Vector para operaciones con cada tipo de capa. Menú File para importación, exportación y conversión de archivos. Menú Databases, para trabajar con tablas, y menú Config para configurar el entorno de trabajo.

○

Ventana de mapas activos: en la que se listan los mapas que queremos visualizar.

○

Configuración de mapas activos: se configuran opciones de visualización de los mapas.

●

Map Display: es el visualizador de mapas, y tiene las típicas opciones de cualquier navegador de mapas, incluida la exportación a distintos formatos de imagen.

●

Output GIS.m: muestra el estado de las operaciones en curso, muestra resultados (aquellos que por su naturaleza se dan en formato texto), y permite lanzar comandos.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 4

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Figura 1: Interfaz gráfico de GRASS, mostrando las ventanas principales y las herramientas más importantes.

2.1.3 Preparación de la región de trabajo Solo nos queda definir para cada localidad la región de trabajo, que indica a GRASS la resolución espacial y los límites geográficos. La resolución espacial corresponde a la del modelo de elevaciones que prepararemos en el apartado siguiente: unos 90 metros por celda, si hablamos de la proyección UTM, o 0.00083333 grados decimales por celda (equivalente a 0º 00' 03'' en notación sexagesimal), si tratamos con coordenadas geográficas. Para definir la resolución espacial, en la ventana Output, en la línea de comandos teclea y ejecuta (pulsando el botón Run): g.region res=0:00:03 Los límites geográficos de la provincia de Almería los obtendremos de un archivo de provincias de toda Europa que se encuentra disponible en la geoplataforma EDIT, del Museo Nacional de Ciencias Naturales (http://edit.csic.es/GISdownloads.html). En esta dirección, en la sección Administrative Units se ha obtenido el fichero provinBlas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 5

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

ce-europe-shape.zip que está en C:\CURSO_MODELOS\ GEODATOS\LIMITES. Descomprímelo en la misma ubicación. Para importar el archivo vectorial de la provincia de Almería (y no el de toda Europa), sigue estos pasos: ●

GIS Manager > File > Import vector map > Multiple formats using OGR, para abrir el módulo v.in.ogr, que se encarga de importar diversos formatos vectoriales. ○

Pestaña Options: ■

Override dataset projection = marcado

■

OGR datasource name = C:/CURSO_MODELOS/GEODATOS/LIMITES/province-europe-shape/province-europe.shp

■

Name for output vector map = ALMERIA

■

WHERE conditions of SQL statement... = IDR_ID=42

■

Type = boundary

○

Run

○

Cierra v.in.ogr.

Para establecer la región de trabajo dentro del área cubierta por el archivo vectorial importado: ●

GIS Manager > Config > Region > Change region settings, para abrir el módulo v.region, que tiene las utilidades para controlar los límites geográficos y resolución espacial de la región de trabajo. ○

Pestaña Existing: ■

Set current region to match this vector map = ALMERIA

○

Run

○

Cierra g.region.

2.2 PREPARACIÓN DEL MODELO DE ELEVACIONES En cualquier trabajo de modelado de distribución, es básico disponer de un modelo de elevaciones con una resolución apropiada para los objetivos del proyecto. En esta práctica vamos a utilizar los resultados del proyecto SRTM. Los datos que vamos a descargar están en coordenadas geográficas (datum= WGS84), y tienen una resolución horizontal de 90m (la resolución vertical ronda los 16m). Es cierto que el modelo de elevaciones que podemos obtener del SRTM no puede considerarse de muy alta resolución, y para trabajos en los que la precisión espacial es importante, es preferible optar por modelos de elevaciones de mayor resolución, pero no siempre están disponibles. Los resultados del proyecto SRTM están disponibles a través de la la web del CGIAR (http:\\srtm.csi.cgiar.org\). Esta institución ha preparado una aplicación para Google Earth desde la que resulta sencillo descargar cualquiera de las teselas de datos disponibles. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 6

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

2.2.1 Descarga del modelo de elevaciones Los pasos a seguir para descargar el modelo de elevaciones (MDE en adelante) son: 1. Cargar el programa de descarga de teselas SRTM en Google Earth: Tras la instalación, descargamos y abrimos la aplicación del CGIAR, que puedes descargar del siguiente enlace (http:\\www.ambiotek.com\srtm). Si no funciona la descarga, también tienes el archivo, llamado srtm3.kmz, en la carpeta C:\CURSO_MODELOS\GEODATOS\ELEVACION. Solo tienes que pulsar dos veces sobre el archivo para que Google Earth lo cargue. 2. Búscar y descargar el archivo: Ahora buscaremos la tesela que corresponde al SE de la Península Ibérica (srtm_36_05), y pulsamos sobre el triángulo verde que presenta en el centro. Pasada la mitad de la página que se abre podemos encontrar el listado de servidores desde el que podemos descargar el conjunto de datos seleccionado. Selecciona el primer servidor listado, y al pulsar sobre él ya podrías iniciar la descarga, pero ya has aprendido a buscar datos SRTM, y no hace falta que los descargues, porque ya están en C:\GEODATOS\ORIGINALES\ELEV\srtm_36_05.zip. Es un archivo comprimido; para extraerlo, solo tienes que pulsar sobre él con el botón derecho del ratón y pinchar sobre extraer todo. Si todo ha ido bien, en la carpeta C:\GEODATOS\ORIGINALES\ELEV\srtm_36_05 habrá un archivo llamado Z_36_5.ASC de unos 170 MB aproximadamente; es nuestro modelo de elevaciones, aunque ocupando un área más extensa que la apropiada para nuestro trabajo. 2.2.2 Importación del MDE la base de datos de GRASS Para empezar a generar nuestras variables, importaremos el modelo de elevaciones a la Localidad en la que estamos trabajando. ●

GIS Manager > File > Import raster map > ESRI GRID, para abrir el módulo r.in.arc, el encargado de importar ficheros en formato .asc. ○

ARC/INFO raster file to be imported = C:/CURSO_MODELOS/GEODATOS/ELEV/srtm_36_05/Z_36_5.ASC

○

Name for output raster map = ELEV

○

Storage type for resultant raster map = CELL

●

Run Tardará un poco, porque está importando el mapa entero, no solo el área de trabajo; es una costumbre de GRASS...En la ventana Output podrás observar el progreso de la operación.

●

Cierra r.in.arc.

Para visualizar el mapa de elevación, hay que seleccionarlo en la ventana GIS Manager, pulsando el icono añadir nueva capa raster. En la ventana de mapas activos aparecerá una entrada para la nueva capa raster. En el panel de configuración de mapas activos pulsa sobre el icono idéntico a abrir nueva capa raster. De la lista que se despliega, selecciona ELEV y pulsa Ok. Para visualizar el mapa, en la ventana de mapas primero pulsa sobre el icono Zoom to (el que tiene una lupa sobre un mapa), y pulsa sobre Zoom display to selected map. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 7

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

2.3 CREACIÓN DE UNA MÁSCARA Una máscara es un mapa que indica a GRASS la áreas a las que debe aplicar los cálculos y a las que no. Las áreas en las que no se aplican cálculos se denominan habitualmente NODATA, pero en GRASS la nomenclatura propia es Null. Para crear la máscara, en esta ocasión, tenemos un ligero problema debido a las distintas fuentes de datos que usamos; la línea de costa definida por el mapa ELEV, y la definida por el contorno provincial del mapa ALMERIA no coinciden. A continuación se describen los pasos para solucionarlo. 1. Vamos a darle valor 1 a todas las áreas con datos válidos del mapa ELEV, guardando el resultado con el nombre MASCARA_ELEV. ○

GIS Manager > Raster > Change category values and labels > Recode interactively, para abrir la ventana Interactive rules entry. ■

Map to recode = ELEV

■

Recoded map = MASCARA_ELEV

■

En la ventana escribir = -30:4000:1 que significa literalmente: recodificar a valor 1 todos los valores dentro del intervalo -30:4000

○

Apply. Saldrá un error, pero puedes ignorarlo.

○

Ok, para cerrar la ventana.

2. Ahora es necesario transformar el fichero vectorial ALMERIA en un raster: ○

GIS Manager > File > Map type conversions > Vector to raster, para abrir el módulo v.to.rast. ■

Name of input vector map = ALMERIA

■

Name for output raster map = ALMERIA

■

Source of raster values = val

■

Type = area

○

Run

○

Cierra v.to.rast.

3. Ahora disponemos de las capas ALMERIA, que indica los límites del área de trabajo “tierra adentro”, y MASCARA_ELEV, que indica los límites de la línea de costa que nos interesa. El área resultante de la intersección entre ambos mapas es nuestro área de trabajo. Para generar la intersección multiplicaremos ambas capas entre sí. En la ventana Output, en la línea de comandos teclea la siguiente expresión y pulsa el botón Run: r.mapcalc MASK=ALMERIA*MASCARA_ELEV Ya tenemos perfectamente delimitada nuestro área de trabajo, puedes comprobarlo visualizando el mapa MASK y usando de nuevo la opción Zoom display to selected map. Esta pequeña obsesión por tener una máscara limpia es debida a un requerimiento de algunos de los programas que usaremos para hacer modelos de distribución: exigen que todos los mapas de variables que utilicemos, independientemente de su origen, coincidan perfectamente, celda a celda. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 8

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

2.4 VARIABLES TOPOGRÁFICAS 2.4.1 Pendientes y orientaciones Para muchas especies la pendiente y la orientación son factores con gran influencia en la distribución espacial, por lo que deben incluirse en cualquier trabajo de modelización, especialmente cuando se utilizan resoluciones altas y regiones poco extensas. Para calcularla con GRASS: ●

GIS Manager > Raster > Terrain analysis > Slope and aspect, para abrir el módulo r.slope.aspect. ○

Pestaña Options: ■

Raster elevation file name = ELEV

■

Output slope filename = PEND

■

Output aspect filename = OR_360

○

Run

○

Cierra r.slope.aspect.

Los valores 0-360 de la variable ORIENTACION no son apropiados para trabajar con ellos en modelos de distribución, porque tanto el valor 0 como el 360 indican el mismo azimut, y esto crea problemas en los algoritmos de modelado. Para transformarla, vamos a crear dos gradientes: ●

OR_NS: con valores 100 en las exposiciones norte, y valores 0 en las sur, con toda la gama intermedia de valores (este y oeste presentan ambos valor 50).

OR_EO: que sigue el mismo esquema que el anterior, pero con valores 100 en la orientación oeste, y 0 en las este (norte y sur presentan ambos valor 50). Esta operación se denomina recodificación (lo hiciste antes, al crear el mapa MASCARA_ELEV), porque implica un cambio en los valores de las celdas; por ejemplo, los valores 360 y 0 de la capa ORIENTACIÓN, serán sustituidos por el valor 100 en la capa OR_NS. Previamente a la recodificación de un mapa raster es necesario crear un fichero que indique los valores que van a variar y sus reemplazos. En GRASS, este fichero debe tener la extensión “.csv” (comma separate values) y el siguiente formato: ●

0:0:100:100 1:1:100:100 2:2:100:100 3:3:99:99 ...hasta completar todos los valores a sustituir. Los dos primeros números de cada fila indican el intervalo de valores del mapa original que se quieren sustituir; los dos últimos números indican el intervalo de valores por el que se quiere cambiar el valor original. La primera línea diría textualmente Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 9

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

“cambia los valores de 0 a 0 por valores de 100 a 100”, o sustituye 0 por 100. Para elaborar estos ficheros es muy útil disponer de una hoja de cálculo (para esto funciona mejor Calc qué Excel), y el editor de texto Notepad++, Los ficheros correspondientes a las dos reclasificaciones que hay que realizar se encuentran en C:\CURSO_MODELOS\GUIONES (GRADIENTE_EO.txt y GRADIENTE_NS.txt). Puedes abrirlos con Notepad++ para observar su contenido (aparecerán modificados si los abres con el bloc de notas de windows). Para recodificar el mapa de orientaciones: ●

GIS Manager > Raster > Change category values and labels > Recode using rules file, para abrir el módulo r.recode.

○

■

Raster map to be recoded = OR_360

■

Name for output raster map = OR_NS

■

File containing recode rules = C:/CURSO_MODELOS/GUIONES/GRADIENTE_NS.csv

Run ■

Raster map to be recoded = OR_360

■

Name for output raster map = OR_EO

■

File containing recode rules = C:/CURSO_MODELOS/GUIONES/GRADIENTE_EO.csv

○

Run

○

Cierra r.recode.

2.4.2 Posición topográfica La posición topográfica de una celda es la diferencia entre la altitud real de la celda, y el promedio de las celdas circundantes en un radio dado. Esta variable contrasta áreas muy expuestas a los elementos, como cimas o crestas, con respecto a áreas resguardadas, como valles, y contrasta muy bien las zonas llanas. Para calcularla son necesarias dos capas: una de elevación, de la que ya disponemos (mapa ELEV), y la elevación media en un radio determinado, que en este caso será de 1000 metros aproximadamente (un radio de 11 celdas de 90 metros cada una). 1. Para calcular la elevación media en un radio determinado: ○

GIS Manager > Raster > Neighbourhood analysis > Moving window, para abrir el módulo r.neighbors. ●

Use circular neighbourhood = marcado

●

Name of input raster map = ELEV

●

Name for output raster map = ELEV_MEDIA_1000m

●

Neighbourhood operation = Average

●

Neighbourhood size = 11 (NOTA: el número debe ser impar)

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 10

Modelos de distribución y estudios de flora amenazada ■

Run

■

Cierra r.neighbors.

Guión de prácticas

2. Para restar la capa ELEV_MEDIA_1000m a la capa ELEV, en el cuadro de ejecución de comandos de la ventana Output teclea y ejecuta: r.mapcalc POS_TOPO=ELEV-ELEV_MEDIA_1000m 3. Para eliminar el mapa ELEV_MEDIA_1000m (no lo necesitarás más), teclea y ejecuta la siguiente expresión en la línea de comandos: g.remove ELEV_MEDIA_1000m 2.4.3 Índice topográfico de humedad La variable Índice Topográfico de Humedad refleja el comportamiento hidrológico del territorio. Esta variable topográfica destaca zonas de suelos profundos (valles, cuencas sedimentarias) de otras con suelos poco profundos, en los que la escorrentía es la norma. Para calcular esta variable GRASS tiene una limitación, y es que solo puede hacerlo con datos en proyección UTM. Por tanto es necesario reproyectar el mapa ELEV antes de realizar el cálculo siguiendo los siguientes pasos: 1. Inactivar la máscara actual, tecleando y ejecutando en la línea de comandos la siguiente expresión: g.rename MASK,MASK_inactiva 2. Cerrar GRASS (cierra GIS Manager y las demás ventanas caerán con ella), y abrirlo de nuevo en la localidad ALMERIA_utm, mapset PERMANENT 3. Es necesario configurar previamente la región de trabajo. Prepararemos la región según un rectángulo que acoge aproximadamente la provincia de Almería (las coordenadas se han sacado de Google Earth): GIS Manager > Config > Region > Change region settings, para abrir el módulo g.region. ○

○

Pestaña Bounds: ■

Value for the northern edge = 4200000

■

Value for the southern edge = 4050000

■

Value for the eastern edge = 625000

■

Value for the wetstern edge = 485000

Pestaña Resolution: ■

Grid Resolution 2D = 90

○

Run

○

Cierra g.region.

4. Para reproyectar el mapa ELEV, en GIS manager > Raster > Develop map > Reproject, para abrir el módulo r.proj. ○

Location of input raster map = ALMERIA_latlong

○

Mapset of input raster map = PERMANENT

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 11

Modelos de distribución y estudios de flora amenazada ○

Name of input raster map to re-project = ELEV

○

Name for output raster map = ELEV

○

Interpolation method = nearest

○

Run

○

Cierra r.proj.

Guión de prácticas

5. Ahora hay que sustituir los valores NULL de ELEV por valores 0 (cero): GIS Manager > Raster > Develop map > Null values, para abrir el módulo r.null. ○

Raster map for which to edit null file = ELEV

○

The value to replace the null value by = 0

○

Run

○

Cierra r.null

6. Para generar el mapa del índice topográfico de humedad, de nuevo en GIS manager > Raster > Hydrologyc modeling > Topographic index map. Se abre el diálogo r.topidx: ○

Input elevation map = ELEV

○

Output topographic index map = ITP

○

Run Tardará un par de minutos.

○

Cierra r.topidx

7. Ahora, este mapa hay que reproyectarlo a latitud-longitud, por lo que cerramos grass, y volvemos a abrirlo en la Location ALMERIA_latlong. 8. Es necesario reactivar la máscara antes de continuar; en la línea de comandos teclea y ejecuta: g.rename MASK_inactiva,MASK 9. Para reproyectar el mapa ITP al sistema de coordenadas latlon: GIS manager > Raster > Develop map > Reproject, y se abre el módulo r.proj. ○

Location of input raster map = ALMERIA_utm

○

Mapset of input raster map = PERMANENT

○

Name of input raster map to re-project = ITP

○

Interpolation method = nearest

○

Run

○

Cierra r.proj

2.4.4 Radiación solar potencial La radiación solar incidente tiene una gran importancia en la distribución de los vegetales, y teniendo un modelo del terreno es posible calcularla, utilizando ecuaciones que simulan el movimiento del sol en fechas y horas determinadas, y teniendo en cuenta el efecto de ocultación de la topografía. Si bien a primera vista estas circunstancias quedarían cubiertas por las variables de orientación, la realidad es que los modelos de radiación solar son mucho más precisos. Para calcular la raBlas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 12

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

diación solar incidente mínima (en torno al 21 de Diciembre, en el solsticio de invierno) seguiremos los siguientes pasos: 1. Crear un mapa de latitudes: es tan sencillo como ir a la línea de comandos y ejecutar la siguiente expresión: r.mapcalc LATITUD=y() 2. Para el cálculo de radiación solar: GIS Manager > Raster > Solar radiance and shadows > Solar irradiance, para abrir el módulo r.sun. ○

○

○

Ficha Options: ■

Incorporate the shadowing effect of terrain = marcado

■

No. of day of the year = 355

■

Time step when... = 1.0 (es para ahorrar tiempo; los resultados son más precisos si este parámetro se configura en 0.5)

Ficha Input options: ■

Name of the input elevation raster map... = ELEV

■

Name of the input aspect... = OR_360

■

Name of the input slope... = PEND

■

Name of the latitudes input raster map = LATITUD

Ficha Output options: ■

○

Output beam irradiance... = RAD_INV

Run, y paciencia, porque tardará un rato!!.

NOTA: Mientras se realiza este cálculo, pasa al punto 2.4.5, y después vuelve aquí. 3. Para calcular la radiación solar máxima (correspondiente al solsticio de verano, en torno al 21 de Junio), se vuelve a ejecutar el módulo, manteniendo los mismos parámetros excepto: ○

No. of day of the year = 175

○

Output beam irradiance... = RAD_VER

2.4.5 Diversidad topográfica La homogeneidad o heterogeneidad topográfica del territorio circundante también es un factor a tener en cuenta. Para calcularla, menú Raster > Neighbourhood analysis > Moving window. ●

Use circular neighbourhood = marcado

●

Name of input raster map = ELEV

●

Name for output raster map = DIV_TOPO

●

Neighbourhood size = 11 (para calcularla en unos 1000 metros a la redonda)

●

Neighbourhood operation = diversity

●

Run

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 13

Modelos de distribución y estudios de flora amenazada ●

Guión de prácticas

Cierra r.neighbors.

NOTA: Vuelve al punto 2.4.4 para comprobar si se terminó el primer cálculo de radiación solar. Si terminó, pon a funcionar el siguiente, y vuelve aquí para continuar.

2.5 VARIABLES DE TELEDETECCIÓN 2.5.1 Importación y procesamiento de bandas Landsat Las imágenes Landsat proporcionan una información multiespectral muy valiosa del territorio. Las longitudes de onda que recoge el sensor ETM+ proporcionan datos sobre la naturaleza geológica del sustrato, la potencia de la cubierta vegetal, los usos del suelo, e incluso la presencia exacta de determinadas formaciones vegetales, y todo con una resolución espacial en torno a los 30 metros. Una fuente que podemos utilizar para descargar un mosaico de imágenes Landsat ya procesado es Image 2000 (http://image2000.jrc.it/), un proyecto europeo de estudio de los usos del suelo (CORINNE Land Cover) para el que se han procesado cientos de imágenes Landsat de toda Europa (años 1998-2000). Es necesario registrarse y explicar la naturaleza del proyecto de investigación para obtener las imágenes (tardarán algunos días en proporcinaros una clave de acceso), pero mientras aquí podéis ver el interfaz de selección y descarga de datos: http://mapserver.jrc.it/website/image2000/viewer.htm. Si haces un zoom sobre el SE de la Península Ibérica, y en el panel de la derecha marcas Mosaic, podrás ver que las cuadrículas U18 y U19 contienen los datos que necesitamos. Ambas teselas están almacenados en C:\CURSO_MODELOS\GEODATOS\ORIGINALES\LANDSAT\U18_MS.bil y U19_MS.bil; descomprime ambas carpetas. Cada uno de los conjuntos de mapas está compuesto por 6 bandas, que representan la reflectancia del terreno ante distintas longitudes de onda. Las imágenes están en una proyección de coordenadas geográficas (latitud-longitud), en un datum llamado ETRS89, que a efectos prácticos, para nuestro trabajo, es equivalente a WGS84 (esto no es cierto para trabajos topográficos de muy alta resolución). Por lo tanto podemos importarlos directamente a nuestra localidad ALMERIA_latlong en GRASS: En la ventana GIS Manager > File > Import raster map > Multiple formats using GDAL, para abrir el módulo r.in.gdal. ●

Pestaña Options: ○

●

●

Override projection = marcado

Pestaña Required: ○

Raster file to be imported = C:/CURSO_MODELOS/GEODATOS/LANDSAT/U18_MS.bil/U18_MS.bil

○

Name for output raster map = U18

Run El conjunto de imágenes U19 se importa del mismo modo, nombrando los re sultados como U19.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 14

Modelos de distribución y estudios de flora amenazada ●

Guión de prácticas

Cierra r.in.gdal.

Ahora tenemos en nuestra base de datos de GRASS 6 parejas de imágenes (U18.1 y U19.1; U18.2 y U19.2, etc...), cada una de las cuales corresponde a una banda Landsat concreta. Las labores que deberíamos realizar ahora son: ●

Unir las imágenes que pertenecen a una misma banda: tenemos las imágenes U18.1 y U19.1, que corresponden a una misma banda Landsat, cada una de las cuáles cubre parcialmente la provincia de Almería. Es necesario unirlas.

●

Borrar las imágenes parciales (U18.1 y U19.1 en el ejemplo anterior).

En total hay que hacer 6 uniones y borrar 12 imágenes, una labor relativamente tediosa. Para hacer operaciones repetitivas, lo más operativo es escribir las órdenes en un fichero de texto plano (script, o guión) para que GRASS las ejecute secuencialmente. En nuestro caso, para unir las imágenes se utiliza la orden r.patch, que tiene la siguiente sintaxis: r.patch input=nombre_imagen_A,nombre_imagen_B output=nombre_imagen_AB Para borrarlas, la expresión también es bastante sencilla: g.remove nombre_imagen_A Ahora que conoces la sintaxis apropiada, puedes escribir y ejecutar las siguientes órdenes secuencialmente en la línea de comandos, cambiando los números correspondientes a cada capa según el caso. r.patch r.patch r.patch r.patch r.patch r.patch

input=U18.1,U19.1 output=LSAT_1 input=U18.2,U19.2 output=LSAT_2 input=U18.3,U19.3 output=LSAT_3 input=U18.4,U19.4 output=LSAT_4 input=U18.5,U19.5 output=LSAT_5 input=U18.6,U19.6 output=LSAT_6 g.remove U18.1 g.remove U18.2 g.remove U18.3 g.remove U18.4 g.remove U18.5 g.remove U18.6 g.remove U19.1 g.remove U19.2 g.remove U19.3 g.remove U19.4 g.remove U19.5 g.remove U19.6

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 15

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

2.5.2 NDVI El Índice Diferencial de Vegetación Normalizado estima la cantidad, calidad y desarrollo de la vegetación a partir de la intensidad de la radiación que la vegetación emite o refleja en distintas bandas del espectro electromagnético. El sensor ETM+ del satélite Landsat recoge información para componer este índice, que se computa a partir de las bandas 3 (rojo) y 4 (infrarrojo cercano) según la siguiente expresión: NDVI = (infrarrojo_cercano – rojo) / (infrarrojo_cercano + rojo) Podemos calcular esta expresión directamente desde la linea de comandos ejecutando la siguiente expresión: r.mapcalc NDVI=1.0*(LSAT_4-LSAT_3)/(LSAT_4+LSAT_3)

2.6 VARIABLES CLIMÁTICAS 2.6.1 Precipitación La precipitación es una variable con gran influencia en la distribución de las especies vegetales, porque representa la disponibilidad de un recurso vital como es el agua. Pero es una variable para la que es muy difícil obtener modelos espaciales y temporales fiables. Distintos métodos proporcionan resultados muy diferentes, a la hora de calcular modelos de precipitación. Para hacer un modelo de precipitaciones de un área concreta, en primer lugar son necesarios los datos de pluviometría y las georreferencias de las estaciones meteorológicas que la cubren. Una vez obtenidos, es necesario realizar el análisis estadístico para obtener los valores descriptivos que nos interesa expresar espacialmente, como podrían ser la precipitación acumulada anual, la precipitación durante el mes mas seco, etc. En la carpeta C:\CURSO_MODELOS\GEODATOS\CLIMA puedes encontrar una hoja de cálculo (PRECIPITACION.ods) con los datos de precipitación de la provincia de Almería. Estos datos son el resultado de un análisis estadístico sobre una base de datos proporcionada por la AEMET, con información meteorológica del periodo 1960-1990. Las coordenadas de las estaciones meteorológicas están en coordenadas geográficas, datum WGS84 (NOTA: las georreferencias de las estaciones meteorológicas proporcionadas por la AEMET están originalmente en ED50. Las que aparecen en las tablas han sido transformadas para este tutorial). Las precipitaciones se muestran en milímetros, y los campos de la tabla son: ●

cat = código identificador de la estación meteorológica

●

LONGITUD = coordenada x

●

LATITUD = coordenada y

●

P_INV = precipitación acumulada media en invierno

●

P_PRI = precipitación acumulada media en primavera

●

P_VER= precipitación acumulada media en verano

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 16

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

●

P_OTO = precipitación acumulada media en otoño

●

P_ANU = precipitación acumulada media anual (suma de las anteriores)

●

P_MAX = precipitación acumulada media del mes más húmedo

●

P_MIN = precipitación acumulada media del mes más seco

●

P_RAN = diferencial de precipitación entre el mes más húmedo y el más seco

Uno método para hacer modelos espaciales de precipitación es la interpolación espacial de los datos puntuales de todas las estaciones meteorológicas. Los pasos a seguir para generar un mapa de precipitaciones mediante interpolación espacial son: 1. Exportar los datos a formato dBASE (dbf): Abre la hoja de cálculo PRECIPITACION con Calc. ○

○

Menú Archivo > Guardar como ■

Nombre: C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\PE RMANENT\dbf\PRECIPITACION.dbf

■

Tipo: Texto dBASE (dbf)

■

Guardar, y pulsa sobre Sí en la advertencia. Marca No volver a mostrar esta advertencia.

■

Cierra la hoja de cálculo.

Ahora es necesario hacer un arreglo: abre con Calc la tabla PRECIPITACION.dbf que acabas de guardar, y en la primera celda sustituye ID,N,5,2 por ID,N,5,0. Ahora puedes guardar los cambios y cerrar Calc.

2. Creación del archivo vectorial con las localidades de las estaciones pluviométricas: GIS Manager > Vector > Generate points > Generate points from database, para abrir el módulo v.in.db. ○

Input table name = PRECIPITACION

○

Name of column containing x coordinate = LONGITUD

○

Name of column containing y coordinate = LATITUD

○

Must refer to an integer column = CAT

○

Name for output vector map = ESTACIONES_PLUVIOMETRICAS

○

Run

○

Cierra v.in.db.

3. Interpolación espacial de los datos de precipitación con el algoritmo Regularized Tension Splines: Gis Manager > Raster > Interpolate surfaces > Regularized spline tension, para abrir el módulo v.surf.rst. ■

Pestaña Options: ●

Name of input vector map = ESTACIONES_PLUVIOMETRICAS

●

Name of the attribute column... = P_ANU

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 17

Modelos de distribución y estudios de flora amenazada

■

■

●

Output surface raster map = PRE_ANU

●

Run

Guión de prácticas

Pestaña Settings: ●

Maximum number of points in a segment = 30

●

Minimum number of points for approximation = 90

Run

Ya has generado tu mapa de precipitación mediante interpolación espacial. Del mismo modo, calcula los siguientes mapas de precipitación: ●

Precipitación mes más seco (campo P_MIN), nombrando el resultado como PRE_MIN.

●

Precipitación del mes más húmedo (campo P_MAX), nombrando el resultado como PRE_MAX.

●

Cierra v.surf.rst.

●

Rango de precipitación anual, ejecutando la siguiente expresión en la línea de comandos: r.mapcalc PRE_RAN=PRE_MAX-PRE_MIN

2.6.2 Temperatura Para calcular la temperatura se podría recurrir a la interpolación espacial, como en el caso de la precipitación, pero lo cierto es que la variable temperatura está muy ligada a la topografía y a la situación geográfica. La mayoría de los autores coinciden en que la temperatura del aire en superficie depende directamente de la altitud, la distancia al mar, la latitud y la longitud. La técnica más utilizada para hacer modelos espaciales de temperatura es la regresión múltiple con corrección de residuos. La idea es calcular la ecuación de regresión múltiple que relaciona la variable respuesta, en este caso podría ser la temperatura media mínima del mes más frío, con el conjunto de variables predictoras, que son: ELEV, LATITUD, LONGITUD, y DIS_MAR. La ecuación debe ser de este tipo: T=m+a*LATITUD+b*LONGITUD+c*ELEV+d*DIS_MAR Esta ecuación proporciona un mapa de temperatura que no tiene en cuenta los valores de los residuales de cada una de las estaciones meteorológicas. Para incorporar los residuos, se interpolan espacialmente los valores de los residuales de cada estación, y el mapa resultante se suma al anterior, para obtener el mapa de temperatura final. 2.6.2.1 Cálculo de variables predictoras

De las variables predictoras de la temperatura, aún no hemos calculado la distancia al mar y la longitud. Para calcular la longitud, igual que calculamos la latitud anteriormente; en la línea de comandos ejecuta la siguiente expresión: r.mapcalc LONGITUD=x() Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 18

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Para calcular la distancia al mar sigue estos pasos: 1. Inactivar la máscara cambiándole el nombre: teclea y ejecuta la siguiente expresión en la línea de comandos. g.rename MASK,MASK_inactiva 2. Invertir el mapa MASCARA_ELEV: GIS Manager > Raster > Develop map > Null Values, para abrir el módulo r.null. ■

Raster map for wich edit null file = MASCARA_ELEV

■

List of cell values to be set to null = 1

■

The value to replace the null value by = 1

■

Run

■

Cierra r.null.

3. Creación de un mapa con valor 9 en todas las celdas (9x10 metros que tiene de lado cada celda; indica el coste en metros al atravesar una celda, pero usamos 9 en lugar de 90 para obtener un mapa menos “pesado”; los mapas, cuanto mayores son los números que almacenan en sus celdas, más espacio en disco ocupan, y más difícil es manejarlos). Lo utilizaremos para crear el mapa de distancia al mar. En la ventana Output tecleamos y ejecutamos la siguiente expresión. r.mapcalc COSTE=9 4. Cálculo del mapa de distancia al mar: GIS Manager > Raster > Terrain analysis > Cost surface, para abrir el diálogo r.cost. ■

Name of raster map containing grid cell cost information = COSTE

■

Name for output raster map = DIS_MAR

■

Starting points raster map = MASCARA_ELEV

■

Run

■

Cierra r.cost

5. Reactivar la máscara: Teclea y ejecuta la siguiente expresión en la línea de comandos: g.rename MASK_inactiva,MASK 2.6.2.2 Cálculo de la ecuación de regresión

En la carpeta C:\CURSO_MODELOS\GEODATOS\ORIGINALES\CLIMA puedes encontrar la hoja de cálculo TEMPERATURA.ods, que contiene los datos de temperatura media anual del periodo 1960-1990 para la provincia de Almería. Los campos que contiene la tabla son: ●

CAT: código de la estación meteorológica

●

LATITUD: coordenada y

●

LONGITUD: coordenada x

●

TMED_ANU: temperatura media anual

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 19

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

●

TMED_INV: temperatura media invierno

●

TMED_PRI: temperatura media primavera

●

TMED_VER: temperatura media verano

●

TMED_OTO: temperatura media otoño

●

TMED_MAX: temperatura media del mes más cálido

●

TMED_MIN: temperatura media del mes más frío

●

TMED_RAN: diferencia entre las medias del mes más cálido y el más frío (rango medio anual)

Estos son los pasos a seguir para calcular la ecuación de regresión: 1. En la tabla TEMPERATURA, entre los campos LONGITUD y TMED_ANU crearemos dos columnas; ELEV y DIS_MAR, y rellenaremos los registros con valores 0,0, cambiando el formato de las celdas para que muestren decimales. Una vez hecho esto, guardamos la tabla en formato dBASE (.dbf). ○

○

Menú Archivo > Guardar como ■

Nombre: C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\PE RMANENT\dbf\TEMPERATURA.dbf

■

Tipo: Texto dBASE (dbf)

■

Aceptar, y cerrar la tabla en Calc, ¡pero no cierres el programa, aún vamos a necesitarlo!.

Ahora es necesario hacer un arreglo: abre con Calc la tabla TEMPERATURA.dbf que acabas de guardar, y en la primera celda sustituye ID,N,5,2 por ID,N,5,0. Ahora puedes guardar los cambios.

2. Creación del archivo vectorial con las localidades de las estaciones termométricas: GIS Manager > menú Vector > Generate points > Generate points from database, para abrir el módulo v.in.db. ○

Input table name = TEMPERATURA

○

Name of column containing x coordinate = LONGITUD

○

Name of column containing Y coordinate = LATITUD

○

Must refer to an integer column = ID

○

Name for output vector map = ESTACIONES_TERMOMETRICAS

○

Run

○

Cierra v.in.db

3. Para calcular la ecuación de regresión necesitamos conocer los valores de elevación y distancia al mar de cada estación termométrica (latitud y longitud ya los conocemos), y añadir esos valores a la tabla de datos de temperatura. ○

GIS Manager > menú Vector > Update point attributes from raster > Sample raster map at point locations, para abrir el módulo v.what.rast. ■

Name of input vector points... = ESTACIONES_TERMOMETRICAS

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 20

Modelos de distribución y estudios de flora amenazada ■

Name of an existing raster map... = ELEV

■

Column name = ELEV

■

Run

■

Name of an existing raster map... = DIS_MAR

■

Column name = DIS_MAR

■

Run

■

Cierra v.what.rast.

Guión de prácticas

4. A continuación realizaremos el análisis de regresión múltiple sobre la tabla del archivo vectorial ESTACIONES_TERMOMETRICAS, para componer la ecuación de regresión y calcular el mapa de temperatura sin corregir. ○

Abre con Calc el archivo C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\PERMANENT\dbf\ESTACIONES_TERMOMETRICAS.dbf

○

Pulsa sobre cualquiera de las celdas en blanco debajo de los datos, y en el menú Insertar, busca la entrada Función, para abrir el asistente Funciones.

○

Busca en la lista la función ESTIMACION.LINEAL, selecciónala, y pulsa Siguiente: ■

Datos_Y = selecciona los valores de la columna TMED_ANU.

■

Datos_X = selecciona los valores de las variables LATITUD, LONGITUD, ELEV, y DIS_MAR.

■

tipo_lineal = 1

■

estadística = 1

■

Aceptar. Debe resultar una matriz similar a la de la Figura 2, que detalla el significado de los valores más importantes.

Figura 2: Matriz de resultados de la función ESTIMACIÓN LINEAL de Calc

Según estos resultados, la ecuación de regresión que explica la temperatura media anual en función de las variables latitud, longitud, elevación y distancia al mar sería: TMED_ANU=1387.94-31.85*LATITUD+8.42*LONGITUD-0.05*ELEV+0*DIS_MAR ○

Para calcular el mapa de temperatura sin corrección de residuos teclea y ejecuta en la línea de comandos de GRASS la siguiente expresión:

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 21

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

r.mapcalc TMED_ANU_reg=1387.94-31.85*LATITUD +8.42*LONGITUD-0.05*ELEV 5. Para calcular los valores puntuales de los residuales, creamos una nueva columna en la tabla abierta en Calc, al final, con el nombre TMED_ANr,N,20,6 y en la primera celda escribimos la ecuación siguiente, en la que restamos el valor estimado por la ecuación de regresión (sección entre paréntesis) al valor de temperatura real de cada registro (celdas de la columna TMED_ANU): =F2-(1387,94-31,85*B2+8,42*C2-0,05*D2+0*E2) ○

Ahora arrastramos esta ecuación a todas las celdas inferiores para calcular el residual de cada registro. Después borramos la matriz de resultados de la función ESTIMACION.LINEAL. Puedes guardar la tabla, para que los datos de los residuales queden escritos en el archivo vectorial ESTACIONES_TERMOMETRICAS.

6. Interpolación espacial de los residuos. Utilizaremos un algoritmo distinto al empleado para calcular la precipitación. En este caso es IDW. ○

GIS Manager > Raster > Interpolate surfaces > IDW from vector points, para abrir el módulo v.surf.idw. ■

Name of input vector map = ESTACIONES_TERMOMETRICAS

■

Name for output raster map = TMED_ANU_res

■

Number of interpolation points = 12

■

Attribute table column with values to interpolate = TMED_ANr

■

Run

■

Cierra v.surf.idw.

7. Suma del mapa obtenido por regresión múltiple y el mapa de residuos para obtener el mapa de temperatura final: Teclea y ejecuta en la línea de comandos la siguiente expresión: r.mapcalc TMED_ANU=TMED_ANU_reg+TMED_ANU_res Ahora que ya no los necesitas, puedes borrar los archivos TMED_ANU_res y TMED_ANU_reg ejecutando la siguiente expresión en la línea de comandos: g.remove TMED_ANU_reg,TMED_ANU_res Hemos visto que el cálculo de mapas de temperatura es laborioso, y por motivos de tiempo, no vamos a calcular ninguno más. Sin embargo, sería ideal disponer de, al menos, el mapa de temperaturas medias máximas del mes más cálido (columna TMED_MAX en la base de datos), el de temperaturas medias mínimas del mes más frío (columna TMED_MIN), y el rango anual de temperatura (TMED_MAX – TMED_MIN).

2.7 LA CORRELACIÓN ENTRE VARIABLES En cualquier tipo de análisis estadístico, introducir variables predictoras excesivamente correlacionadas entre sí puede llevar a obtener resultados difíciles de Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 22

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

interpretar. Cuando se prepara un conjunto de variables para alimentar modelos de distribución, también es necesario tener en cuenta esta cuestión. Y resulta normal que muchas de las variables estén correlacionadas entre sí, porque están derivadas de un mismo modelo de elevaciones, porque existen similitudes entre ellas, como en las bandas Landsat. A continuación vamos a analizar la correlación espacial (utilizando el coeficiente de correlación de Pearson) que existe entre las variables que hemos generado. GRASS dispone de una interfaz gráfica para un módulo, r.covar , que calcula la matriz de correlaciones de n mapas, pero no la exporta de modo que la podamos utilizar. Vamos a utilizar la versión de línea de comandos de r.covar, para exportar la matriz de correlaciones a un archivo de texto. La sintaxis de la orden es la siguiente (la sección iniciada con “>>” se utiliza para redirigir el resultado a un fichero en la dirección especificada): r.covar -r map=mapa_1,mapa_2,...,mapa_n >> C:/CURSO_MODELOS/matriz_correlacion.txt Los mapas que introduciremos en el análisis de correlación son únicamente aquellos que tenemos intención de utilizar como variables predictivas en los modelos de distribución. Por tanto, los mapas ALMERIA, DIS_MAR, LATITUD, LONGITUD, MASCARA_ELEV, MASK y OR_360 no será necesario añadirlos en la expresión. 1. Una buena idea para reducir el tiempo de cálculo es reducir temporalmente la resolución a la mitad ejecutando en la línea de comandos: g.region res=0:00:06 2. Teclea y ejecuta (es muy posible que en la ventana output aparezca que la operación se ha detenido con un error. Realmente no es así, y en un par de minutos el cálculo estará terminado): r.covar -r map=DIV_TOPO,ELEV,ITP,LSAT_1,LSAT_2,LSAT_3,LSAT_4,LSAT_5,LSAT_ 6,NDVI,OR_EO,OR_NS,PEND,POS_TOPO,PRE_ANU,PRE_MAX,PRE_MIN,PRE_R AN,RAD_INV,RAD_VER,TMED_ANU >> C:/CURSO_MODELOS/matriz_correlacion.txt 3. Restaura la resolución original; teclea y ejecuta: g.region res=0:00:03 4. Abre con Calc el archivo matriz_correlacion.txt, que está en C:/CURSO_MODELOS. En el diálogo Importación de texto, marca la opción Espacio en Opciones de separación, desmarcando las demás. En la opción A partir de línea, escribe 40. ○

Selecciona las celdas con datos, y sustituye los puntos por comas (menú Editar > Buscar y reemplazar).

○

Crea una nueva primera fila y columna en blanco; en la primera fila, escribe secuencialmente los nombres de las variables tal y como estaban para el cálculo de la matriz de correlación: DIV_TOPO ELEV ITP LSAT_1 LSAT_2 LSAT_3 LSAT_4 LSAT_5 LSAT_6 NDVI OR_EO OR_NS PEND POS_TOPO PRE_ANU PRE_MAX PRE_MIN PRE_RAN RAD_INV

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 23

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

RAD_VER TMED_ANU. ○

Copia la fila con los nombres de las variables, y pégala en la primera columna, utilizando la opción Transponer del Pegado especial. Ahora tienes la matriz de correlaciones completa. Para mejorar la visibilidad, puedes borrar todos los valores por encima o por debajo de la diagonal de “unos”.

○

Es el momento de establecer un umbral máximo de correlación entre dos variables. No existe una norma concreta al respecto, pero lo habitual es establecer este umbral dentro del intervalo [0.5, 0.8]. En este caso lo estableceremos en 0.7 (incluyendo -0.7). Para localizar estos valores en la matriz, lo mejor es utilizar la función Formato condicional, para destacar los valores por encima del umbral. ■

Pulsa la tecla F11 para abrir el cuadro Estilo y formato.

■

Crea un nuevo estilo de celda, usando el botón derecho del ratón sobre el cuadro (zona blanca) y pulsando Nuevo. ●

Pestaña Administrar ○

●

Pestaña Fondo ○

● ■

Nombre = CORRELACION Selecciona un color a tu gusto para destacar los valores por encima del umbral.

Aceptar

Selecciona todas las celdas con valores, y en el menú Formato, abre el diálogo Formato condicional. ●

Condición 1: El valor de la celda mayor o igual que 0,7; Estilo de celda = CORRELACION

●

Condición 2: El valor de la celda menor o igual que -0,7; Estilo de celda = CORRELACION

●

Aceptar.

Si observamos las celdas destacadas, observamos que: ●

Las variables PEND y DIV_TOPO tienen una correlación igual 0.78; seleccionaremos PEND como representante de esta pareja.

●

Las variables TMED_ANU y ELEV tienen una correlación igual -0.97; seleccionaremos TMED_ANU como representante de esta pareja.

●

Las variables OR_NS y RAD_INV tienen una correlación igual -0.78; seleccionaremos RAD_INV como representante de esta pareja.

●

Las variables PRE_ANU y PRE_MAX tienen una correlación igual 0.88; seleccionaremos PRE_MAX como representante de esta pareja.

●

Las variables PRE_RAN y PRE_MAX tienen una correlación igual 0.99; seleccionaremos PRE_MAX como representante de esta pareja.

●

Las capas Landsat tienen todas un índice de correlación mayor de 0.7. Esta cuestión la resolveremos en el siguiente apartado.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 24

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Observando estos resultados, resulta que para mantener una correlación entre capas inferior al umbral establecido, deberíamos descartar para trabajar las capas: DIV_TOPO, ELEV, OR_NS, PRE_ANU y PRE_RAN. En adelante, no trabajaremos más con ellas. Entenderás que fue necesario tomarse el trabajo en generarlas, para llegar a este punto... Una cuestión importante: si para el trabajo que estamos desarrollando nos interesara expresamente conocer la importancia de la precipitación anual (PRE_ANU) en la distribución de la especie, no deberíamos descartar esta variable, sino la correlacionada con ella.

2.8 ANÁLISIS DE COMPONENTES PRINCIPALES Un problema frecuente, derivado del análisis de imágenes multiespectrales, como las Landsat con las que estamos trabajando, es el de la correlación existente entre ellas; representan mucha información redundante, y esto no aporta nada nuevo, e incrementa el tiempo de cálculo de los modelos de distribución. El análisis de componentes principales aplicado a las bandas Landsat busca comprimir toda la información contenida en las bandas originales, en un nuevo conjunto más pequeño, sin perder la información original más significativa. GRASS posee un módulo específico para hacer esta reducción de datos, i.pca: 1. GIS Manager > Imagery > Transform image > Principal Components, para abrir el módulo i.pca.

●

○

Name of input raster maps = LSAT_1,LSAT_2,LSAT_3,LSAT_4,LSAT_5,LSAT_6

○

Name for output raster map = cpLSAT (cp, de componente principal)

Run

El resultado es un conjunto de 6 imágenes con el prefijo cpLSAT, que representan los factores obtenidos mediante el análisis de componentes principales. Tienen una correlación baja entre sí, y de la primera a la última, recogen sucesivamente la varianza del conjunto de datos. Habitualmente, el primer factor ya recoge un porcentaje importante de la misma (80% o más). Nosotros utilizaremos cpLSAT_1 y cpLSAT_2 como sustitutas del conjunto de datos completo. Puedes borrar las restantes: g.remove cpLSAT.3,cpLSAT.4,cpLSAT.5,cpLSAT.6 Ahora eliminaremos los puntos de los nombres de las dos nuevas variables: g.rename cpLSAT.1,cpLSAT_1 g.rename cpLSAT.2,cpLSAT_2 En este momento ya hemos terminado de generar las variables con las que calibraremos los modelos de distribución. El listado completo es el que sigue: ●

ITP

●

NDVI

●

OR_EO

●

PEND

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 25

Modelos de distribución y estudios de flora amenazada ●

POS_TOPO

●

PRE_MAX

●

PRE_MIN

●

RAD_INV

●

RAD_VER

●

TMED_ANU

●

cpLSAT_1

●

cpLSAT_2

Guión de prácticas

Sin embargo, aún tenemos algún preprocesamiento más que hacer...

2.9 CAMBIO DE RESOLUCIÓN DE LAS VARIABLES 2.9.1 Introducción Disminuir la resolución espacial de las variables en este caso no debería ser necesario, porque lo ideal es buscar siempre una máxima resolución de trabajo. Sin embargo, para este curso, en algún momento será necesario acelerar el tiempo de cómputo de los modelos de distribución, por lo que vamos a generar conjuntos de variables a dos resoluciones distintas: 90 metros y 1000 metros. Hasta el momento se ha trabajado con una resolución horizontal de 90 metros (0º 00' 03''), tal y como está establecida en la región de trabajo de GRASS. Esta sería la resolución de las variables para los modelos de distribución si las exportáramos en este momento. GRASS permite cambiar la resolución de trabajo muy rápidamente, ejecutando el comando g.region res=n, siendo n la resolución en formato g/m/s (por ejemplo, 0:00:03, que es la resolución actual de trabajo). Por tanto, si cambiamos la resolución de trabajo con g.region res=0:00:06, a partir de este momento estaríamos trabajando a una resolución en torno a los 180 metros, y las variables se exportarían en esa resolución. Esta operación rápida de cambio de resolución afecta levemente a la calidad de las variables si el cambio en resolución es pequeño, pero puede ser grave si el cambio de resolución es significativo (por ejemplo, de 0:00:03 a 0:00:30), porque el algoritmo de interpolación que utiliza GRASS por defecto, llamado Vecino Más Próximo (Nearest Neighbor) funciona de modo muy simple (para ganar velocidad), sesgando la información: imagina cuatro celdas con los valores de altitud (observa que se trata de una variable contínua) 1, 2, 3 y 4, teniendo en cuenta que la resolución espacial es 0:00:03: 12 34 Al cambiar la resolución a 0:00:06, estas celdas se transformarían en una sola. Cuando GRASS utiliza el interpolador Nearest Neighbor, el valor de la nueva celda sería igual al de la celda superior izquierda (porque es la “vecina más próxima” según GRASS lee el mapa para cambiar la resolución) del conjunto de cuatro celdas, que en este caso es 1. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 26

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Imagina que utilizamos un interpolador que en lugar de seleccionar como representante de las cuatro celdas la primera de ellas, hiciera el promedio de las cuatro; en este caso, el valor de la nueva celda sería 2.5. Aparentemente este nuevo valor es más representativo de la altitud en la zona que el 1 que obtuvimos con el método Nearest Neighbor. Supongamos ahora una variable cualitativa, como podría ser un mapa de usos del suelo, en el que cada valor representa un uso diferente: 12 33 A simple vista, si tuviéramos que representar el valor de estas cuatro celdas en una sola, el más apropiado sería 3 (la moda), porque es el valor dominante. 2.9.2 Cambiando resoluciones GRASS dispone de un módulo específico capaz de hacer estos cambios de resolución basados en estadísticas de celdas: r.resamp.stats (GIS Manager > Raster > Develop map > Resample using aggregate statistics). Este módulo transforma la capa que se le indique A LA RESOLUCIÓN ACTUAL DE LA REGIÓN. Por tanto, si queremos cambiar capas de resolución, antes hay que cambiar la resolución de la región; teclea y ejecuta en la línea de comandos: g.region res=0:00:30 Ahora abre el módulo r.resamp.stats: ●

Name for input raster map = ITP

●

Name for output raster map = 1000m_ITP

●

Aggregation method = Average

●

Run

Ahora deberías repetir el proceso con las demás capas, para pasarlas todas a una resolución de 1000 metros, pero para ahorrar trabajo, en la carpeta GUIONES, está el guión CAMBIA_RESOLUCION.txt, que preparará todas las variables que vamos a utilizar para los modelos de distribución a la resolución de 1000m. Para ejecutar el guión: 1. Menú de inicio de Windows > Todos los programas > GRASS > GRASS Command Line, para abrir la versión de GRASS en línea de comandos, que permitirá la ejecución del guión. 2. Comprueba que la Location y el mapset son correctos (ALMERIA_latlong y PERMANENT), y pulsa secuencialmente las teclas ESC y ENTER. 3. Ordena a la línea de comandos que trabaje en el directorio GUIONES: cd C:\CURSO_MODELOS\GUIONES 4. Ejecuta el guión, escribiendo la siguiente expresión y pulsando ENTER: sh CAMBIA_RESOLUCION.txt

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 27

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

2.10 EXPORTACIÓN DE LAS VARIABLES AMBIENTALES El trabajo de generar variables ambientales para alimentar a los modelos de distribución ya está realmente terminado. Ahora vamos a exportarlas a un formato que el programa de modelización “entiendan”: ESRI ASCII GRID, o ARC/INFO ASCII (extensión .asc). ●

GIS Manager > menú File > Export raster map > ESRI ASCII GRID, para abrir el módulo r.out.arc. ○

Name of an existing raster map layer = nombre_mapa

○

Name of an output ARC-GRID map = C:/CURSO_MODELOS/GEODATOS/VARIABLES/nombre_mapa.asc

○

Number of decimal places = 2

●

Run

●

Deberías repetir la acción sucesivamente cambiando nombre_mapa por el nombre de cada una de las variables, para las tres resoluciones disponibles. Como el proceso es largo, en la carpeta GUIONES está el guión EXPORTAR_VARIABLES.txt, que envía las variables en formato .asc (ESRI ASCII) a C:\CURSO_MODELOS\GEODATOS\VARIABLES dentro de la cuál están las carpetas 90m, y 1000m, que almacenarán respectivamente las variables de cada resolución.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 28

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

SECCIÓN 3: PREPARACIÓN DE LOS REGISTROS DE PRESENCIA (Tiempo estimado: 20 minutos)

3.1 Tipos y calidad de los datos Los modelos de distribución deben alimentarse con registros de presencia de la especie de trabajo. El origen de estos registros puede ser variado; de mayor a menor calidad: 1) marcas de GPS recogidas en muestreos sistemáticos; 2) marcas de GPS recogidas casualmente durante salidas de campo; 3) polígonos o puntos dibujados sobre ortoimágenes de alta resolución en un SIG; 4) marcas a mano sobre un mapa 1:10000 o 1:50000; 5) cuadrículas UTM de 1x1 o de 10x10 Km; 6) topónimos. En el caso que nos ocupa, de variables con alta resolución espacial, únicamente son útiles las tres primeras opciones, que están indicadas para tamaños de celda a partir de los 10 metros. En la recogida de datos con el GPS es importante que la señal de los satélites sea clara, respetar los tiempos de toma de la georreferencia, que el datum esté correctamente configurado y comprobar que el sistema EGNOS (o WAAS) esté activado (permite una resolución espacial en torno a los 3-6 metros). Los puntos y polígonos sobre ortoimágenes de alta resolución son apropiados cuando se conoce perfectamente el terreno, y se pueden localizar con precisión puntos de referencia. La utilidad de esta técnica está fuera de duda cuando se trabaja con especies de porte grande, como arbustos o árboles. Las marcas a mano sobre mapa, trabajando en el terreno, están sujetas a errores de apreciación grandes, con un rango de fiabilidad entre 50 y 100 metros sobre cartografías 1:10000, y entre 100 y 300-400 metros sobre cartografías 1:50000. En estos casos hay que adaptar la resolución espacial de las variables al radio de error estimado para cada georreferencia, y siempre es más fiable trabajar sobre una resolución algo superior al doble de este radio. Las citas de cuadrículas UTM de 1x1, como la mayoría de los registros de herbario (o de GBIF), deben ser utilizadas con variables de al menos 1 kilómetro de resolución horizontal, pero sería importante adaptar las celdas de las variables a la malla UTM 1X1, para evitar errores que pueden llegar a ser de 2 kilómetros en caso de una mala coincidencia entre capas. Lo mismo puede aplicarse con citas en cuadrículas UTM de 10x10, con variables de 10 kilómetros de resolución. Para solucionar el problema de los topónimos (con el que es mejor no encontrarse), se han diseñado los Nomenclátor, que permiten ubicar en un mapa topónimos concretos, o proporcionan directamente las coordenadas más aproximadas. En la página de herramientas geográficas del nodo español de GBIF podéis encontrar algunas de estas herramientas (http://www.gbif.es/herramgeo.php). En la ruta C:\CURSO_MODELOS\GEODATOS\PRESENCIAS puedes encontrar la hoja de cálculo PRESENCIAS.ods. Tiene tres hojas: GPS, con datos tomados en campo mediante GPS; GBIF, con registros de herbario (cuadrículas 1x1) y TALLER, en la que prepararemos los distintos archivos. A partir de los datos de estas hojas trabajaremos los ejemplos del curso. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 29

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

3.2 Números mínimo y máximo de presencias Sobre el tamaño de muestra mínimo para hacer un modelo útil no existe acuerdo, y puede variar según la técnica de modelización empleada. Como norma general, es difícil trabajar con normalidad utilizando menos de 10 registros de presencia. Se pueden conseguir modelos bastante aceptables con tamaños de muestra entre 10 y 30, siendo esta última cifra el mínimo óptimo. Sin embargo esto dependerá de la coherencia ecológica que tengan entre sí esos puntos de presencia. Imagina 10 puntos de presencia de una especie de amplia distribución y gran plasticidad ecológica, con distintos valores de las condiciones ambientales medidos en cada localidad. Difícilmente podrá sacar ningún algoritmo un modelo válido en estas circunstancias. Sin embargo, 5 únicas localidades coherentes entre sí dan datos mínimos pero suficientes para calibrar un modelo de cierta utilidad. Por tanto se puede decir que la calidad del modelo no solo va a depender del tamaño de la muestra, sino de la coherencia de la misma, la naturaleza de la especie, y el modo en el que el algoritmo de modelado trata los datos. No puede hablarse de un número máximo de registros de presencia, porque no es habitual trabajar con más de 1000 puntos. Hasta donde nuestra experiencia nos dice, algunos algoritmos no presentan problemas trabajando con más de 9000 registros de una misma especie.

3.3 Densidad de los datos Es una cuestión a la que se le suele prestar poca atención, pero que en algunos algoritmos de modelado presenta una gran importancia. Una agregación espacial desequilibrada de los datos de presencia puede degradar el resultado, por una ponderación por exceso para las condiciones ambientales que se dan en las áreas en las que hay mayor densidad de registros de presencia (ver Figura 3).

Figura 3: Posibilidades de equilibrar la densidad de muestras de presencia

Esta situación se da fácilmente cuando transformamos un polígono en puntos con una resolución elevada. Si el polígono coexiste con otros registros puntuales, aparece un desequilibrio en la densidad de la muestra (CASO 1, en la Figura 3). En este Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 30

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

caso es recomendable disminuir la densidad de puntos dentro del polígono, bien haciendo un muestreo aleatorio dentro del mismo (CASO 2), o un muestreo sistemático (CASO 3).

3.4 Partición aleatoria de datos Para evaluar los modelos de distribución, una entre las distintas técnicas disponibles para obtener datos independientes es la partición aleatoria de datos. Los registros de presencia se dividen al azar en dos grupos, uno de los cuáles se utilizará para calibrar el modelo (conjunto de calibrado), y otro para evaluarlo (conjunto de evaluación). Cuando se dispone de más de 10 registros de presencia, lo normal es utilizar un 75% para calibrar el modelo y un 25% para evaluarlo (aunque sería recomendable recalibrar el modelo con todos los registros de presencia). Si en número de registros es mucho mayor, se puede llegar a utilizar un 70% de calibrado y un 30% de evaluación, o incluso un 60/40. Un método fácil para seleccionar aleatoriamente los puntos de calibrado y evaluación con Calc consiste en utilizar la función ALEATORIO, como veremos en el siguiente apartado.

3.5 Preparación de datos de presencia El programa con el que vamos a hacer modelos de distribución acepta las coordenadas de las especies en formato de texto separado por tabulaciones, con extensión .txt. Debe ir una especie por archivo de texto, con el siguiente formato: #id 1 2 3 ...

label long lat linaria nigricans linaria nigricans linaria nigricans

abundance ­2.35582 37.07454 ­2.38859 37.11790 ­2.34076 37.10853

1 1 1

El modo más rápido de crear un archivo de texto con los campos separados por tabulaciones es copiar y pegar los datos directamente desde una hoja de cálculo al archivo de texto plano en blanco: 1. En la carpeta C:\CURSO_MODELOS\GEODATOS\PRESENCIAS crea un archivo de texto en blanco (pulsa el botón derecho del ratón sobre el explorador abierto en la carpeta, abre el menú nuevo, y pulsa sobre documento de texto), nombrándolo Linaria_nigricans.txt. Ábrelo con Notepad++ (botón derecho sobre el icono del archivo, y pulsar sobre Edit con Notepad++). 2. Abre la hoja de cálculo PRESENCIAS.ods, y copia todos los registros de la especie Linaria nigricans de la hoja GPS en la segunda fila y columna (celda B2) de la hoja TALLER. 3. En las celdas de la primera fila escribe la cabecera del fichero: #id, label, long, lat, abundance 4. En la columna de la derecha (la siguiente a abundance, columna F) pega y Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 31

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

arrastra hacia abajo la función ALEATORIO, hasta que todos los registros tenga su valor aleatorio. 5. La función aleatorio tiene un comportamiento incómodo: sus valores se actualizan cada vez que modificamos algo en la hoja de cálculo; para fijarlos, cópialos y pégalos en la columna de la derecha usando la opción Pegado Especial, desmarcando Pegar todo y Fórmulas, y marcando Números. Elimina la columna original en la que insertaste la función ALEATORIO. 6. Selecciona todas las columnas con datos, y en el menú Datos, selecciona Ordenar. ○

Pestaña Opciones: ■

○

Pestaña Ordenar por criterios: ■

○

El intervalo contiene etiquetas de columnas = marcado Ordenar según = Columna F

Aceptar

7. En la columna A (cabecera #id) crea números correlativos para identificar cada registro. 8. Rellena la columna abundance con valores 1. 9. Elimina la columna F, en la que están los valores aleatorios. 10. Separaremos el conjunto de datos en dos conjuntos, uno con el 60%, y otro con el 40% (111 y 75 respectivamente). En la última columna, hasta el registro 112, escribe calibrado. Desde el registro 113 hasta el final, en la misma columna escribe evaluación. 11. Cambia el nombre de la hoja de TALLER a LNIGRICANS. 12. Copia las cinco columnas (olvida por el momento esa en la que has escrito calibrado y evaluación; más adelante veremos que significado y utilidad tiene todo esto) desde la cabecera hasta el registro 112 y pégalos en el documento de texto Linaria_nigricans.txt. Ahora puedes guardar el archivo de presencias de que has creado, y la hoja de cálculo.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 32

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

SECCIÓN 4: MODELOS DE DISTRIBUCIÓN CON OPENMODELLER (Tiempo estimado: 20 minutos) El programa OpenModeller es gratuito, y de código abierto. Su principal objetivo es proporcionar a los investigadores una herramienta para generar modelos espaciales de nicho ecológico. Permite generar y comparar múltiples modelos simultáneamente. Incluye varios algoritmos diferentes, y su arquitectura modular facilita la incorporación de nuevos algoritmos con un escaso coste de desarrollo. El interfaz gráfico del programa permite de modo fácil preparar las capas de variables, gestionar los parámetros de los distintos algoritmos, diseñar experimentos con múltiples modelos, y visualizar los resultados, proporcionando un completo informe de cada uno de ellos. Para comenzar vamos a crear varios modelos de distribución de la planta Linaria nigricans, a una resolución de 1000 metros. Durante la ejecución aprovecharemos para explicar las particularidades de cada uno de los modelos que están implementados en OpenModeller. Inicia el programa, que comenzamos.

4.1 Preparación de las variables 1. Menú Edit > Preferences ○

Pestaña Saving: ■

Data Directory = C:/CURSO_MODELOS/RESULTADOS

■

Pulsa OK.

2. Menú Edit > Layer sets: ○

Name = ALMERIA_1000m

○

Pulsar sobre el símbolo “+” ■

Se abre el diálogo Raster Layer Selector ●

Como Base Dir busca y selecciona C:/CURSO_MODELOS/GEODATOS/VARIABLES

●

Despliega el contenido de la carpeta 1000m, selecciona todas las capas y pulsa OK.

○

Asegúrate de que la entrada de la sección Mask layer corresponde a alguna de las capas del conjunto que has seleccionado. Esta es la capa que indicará a OpenModeller el área de trabajo y la resolución.

○

Pulsa Apply y Close.

○

Ahora repite la operación para las variables de 90m; nombra el conjunto de capas como ALMERIA_90m. En este momento tenemos todas las capas preparadas para su utilización con OpenModeller.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 33

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

4.2 Ejecución de un experimento de prueba 1. Menú File > New Experiment, para abrir el diálogo Experiment Designer: ○

Name = L_nigricans_1000m

○

Coordenate System = Lat/Long WGS84

○

Ocurrence Data: ■

Pulsar sobre el símbolo “+”; buscar y seleccionar la carpeta C:/CURSO_MODELOS/GEODATOS/PRESENCIAS. Debe aparecer el fichero Linaria_nigricans.txt seleccionado.

○

Agorithms: selecciona todas las entradas (excepto AquaMaps y las tres últimas, que son variantes de Support Vector Machines).

○

Model Settings

○

■

Creation Layers = ALMERIA_1000m

■

Projection Layers = ALMERIA_1000m

■

Output Format = cualquiera de las capas de C:/CURSO_MODELOS/GEODATOS/VARIABLES/1000m

Ok. Preguntará si quieres ejecutar el experimento ahora; contesta Yes.

Al estar las capas a baja resolución, el experimento se va a ejecutar con cierta rapidez. Cuando la barra de progreso haya llegado al 100%, pulsa OK, y observa el informe de resultados. Vamos a conocer un poco más la lógica interna de cada uno de los algoritmos.

4.3 Lógica interna de los distintos algoritmos 4.3.1 Bioclim Para cada una de las variables se define una envuelta bioclimática cuadrangular para la especie, calculando la media y la desviación estándar de los valores de los registros de presencia sobre la variable, asumiendo una distribución normal. De este modo, para cada variable se tiene en cuenta una envuelta bioclimática óptima para la especie, representada por el intervalo [m-c*s, m+c*s], donde m es la media, s es la desviación estándar, y c un parámetro que define qué porcentaje del rango de desviación estándar queda dentro de la envuelta. Por defecto está establecido en 0.674 que corresponde con un 50% del rango de la desviación estándar, mientras que un valor 3.000 corresponde con un 99.7% del rango de desviación estándar incluido dentro de la envuelta óptima. Este parámetro, por tanto, controla la extensión de la envuelta bioclimática óptima. Adicionalmente se define una envuelta subóptima, a partir de los valores máximo y mínimo de los puntos de presencia sobre la variable.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 34

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

En la Figura 4 queda gráficamente explicado el concepto de envueltas óptima y subóptima. Incrementando el valor del parámetro c, se incrementa el área azul claro (envuelta óptima), aumentando el número de celdas con predicción favorable hacia la presencia de la especie.

Figura 4: Descripción gráfica del funcionamiento de BIOCLIM.

Cada celda del mapa de idoneidad resultante se clasifica según los siguientes criterios: ●

Apropiada para la especie (celdas en color verde): los valores de todas las variables quedan dentro de la envuelta bioclimática óptima para la especie.

●

Marginalmente apropiada para la especie (celdas en color rojo): los valores de al menos una de las variables quedan fuera de la envuelta bioclimática óptima, pero dentro de la subóptima.

●

No apropiada (celdas en color azul): los valores de al menos una de las variables quedan fuera de la envuelta subóptima para la especie.

4.3.2 Climate Space Model Desgraciadamente, la documentación de este algoritmo está incompleta y es ciertamente confusa. Solo puede decirse con certeza que es un modelo de distribución basado en el análisis de componentes principales. El algoritmo estandariza los valores de cada una de las variables ambientales sobre el conjunto de registros de presencias y calcula la matriz de covarianzas. Posteriormente calcula los valores propios y el vector propio de la matriz de covarianzas. Del vector propio se eliminan secuencialmente todos los factores cuyo eigenvalue está por debajo del definido por el parámetro “Number of standard deviations” del modelo (método Broken-Stick). Para generar el mapa de idoneidad del hábitat, el algoritmo en primer lugar estandariza las variables, y las multiplica por sus respectivos vectores propios para generar los componentes. Los parámetros para calibrar el algoritmo son: Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 35

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

●

Number of random eigenvalues: entero, dentro del rango [1, 1000], con un valor típico de 8. Es el número de veces que el algoritmo baraja al azar las filas de la matriz de variables para obtener un valor propio aleatorio, que se usa como referencia para la selección del número de componentes a retener.

●

Number of standard deviations: real, dentro del rango [-10, 10], con un valor típico de 2.0. Cuando los valores propios de la matriz barajada al azar se suman, este número es añadido a la media de todos los valores propios. Cualquier componente cuyo valor propio está por encima de el umbral establecido se retiene.

●

Minimum number of components in model: entero, dentro del rango [1, 20], con un valor típico de 3. Determina el número mínimo de componentes en el modelo.

El mapa de idoneidad resultante muestra un gradiente de colores, indicando los más cálidos zonas idóneas para la especie. 4.3.3 Envelope Score Este algoritmo, similar a Bioclim, usa los valores mínimos y máximos observados de cada variables para definir las envueltas bioclimáticas. La probabilidad de ocurrencia de la especie en cada celda se define según la expresión p = nº variables dentro del criterio min-max / nº total de variables. Si se tienen en cuenta únicamente las celdas con probabilidad 1, el resultado es similar al obtenido con Bioclim, 4.3.4 Environmental Distance Este algoritmo se basa en métricas de similaridad ecológica. Puede utilizar cuatro métricas diferentes: Euclídea, Mahalanobis, Manhattan/Gower y Chebyshev. Los parámetros que presenta son: ●

Nearest 'n' points (puntos vecinos): Con rango [0, oo], y un valor típico de 1. Se refiere al número de puntos cercanos que se utilizan como referencia cuando se calculan las distancias ecológicas. Cuando el parámetro está establecido en 1, las distancias se miden a el punto de presencia más cercano. Cuando el parámetro está establecido en 0, las distancias se miden respecto a la media de todos los puntos.

●

Maximum distance: con rango [0, 1] y valor típico de 0.1, describe la máxima distancia de referencia para definir el espacio ecológico fuera del cual las condiciones se consideran no apropiadas para la especie. El valor 1 corresponde con la máxima distancia posible entre dos puntos en el espacio ecológico, por lo que utilizar este valor supone que todos los puntos del espacio geográfico tendrán asociada una probabilidad de presencia en el resultado. La probabilidad de presencia desciende linealmente con la distancia.

OpenModeller trae preparado el algoritmo separando las cuatro métricas, por lo que en el menú de modelos a ejecutar podemos seleccionarlas por separado. Distancia euclídea (http://en.wikipedia.org/wiki/Euclidean_distance): Se refiere a la distancia “convencional” más corta entre dos puntos de un espacio de n-dimensioBlas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 36

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

nes. Mide la distancia entre los valores de las variables de la celda y los valores de las variables del punto de presencia más cercano, expresándola como un índice de similaridad. Distancia de Chebyshev (http://en.wikipedia.org/wiki/Chebyshev_distance): Es la máxima distancia posible entre dos puntos, en cualquier dimensión de los ejes de coordenadas. Se denomina también “distancia del tablero de ajedrez” (ver Figura 5). La siguiente figura es bastante explicativa al respecto; la distancia del rey a las esquinas es la misma que al lado del tablero al que pertenece la esquina.

Figura 5: Ejemplo gráfico de como se mida la distancia de Chebyshev

Distancia de Mahalanobis (http://es.wikipedia.org/wiki/Distancia_de_Mahalanobis): Se utiliza para determinar la similitud entre dos variables aleatorias multidimensionales. Se diferencia de la distancia Euclídea en que tiene en cuenta la correlación entre las variables aleatorias. Distancia de Manhattan (http://en.wikipedia.org/wiki/Manhattan_distance): Es la suma de las diferencias absolutas de las coordenadas de los dos puntos entre los que se quiere medir la distancia. La Figura 6 resume el concepto.

Figura 6: Ejemplo gráfico explicando diferencias en la medida de distancias de Manhattan y euclídea

4.3.5 GARP GARP (Genetic Algorithm for Rule Set Prediction) es un algoritmo genético que genera modelos que describen las condiciones ecológicas bajo las que la especie puede mantener poblaciones viables. Un algoritmo genético es un método de inteligencia artificial inspirado en el concepto de la evolución guiada por la selección natural. La idea clave es tratar de desarrollar Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 37

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

soluciones a los problemas utilizando el mismo mecanismo que los organismos vivos utilizan para evolucionar. Con este objetivo se crea un conjunto de soluciones potenciales y alternativas a un problema (por ejemplo, distintos genotipos generados por mutaciones), y posteriormente, y de modo iterativo, se modifica y testea este conjunto, hasta que una solución óptima se encuentra (supervivencia del genotipo más apropiado al ambiente). GARP busca iterativamente correlaciones no aleatorias entre la presencia de las especies y las variables ambientales, utilizando distintos tipos de reglas. Cada regla implementa un método diferente de construir modelos de nicho ecológico. Las reglas son atómica, rango regresión logística, envuelta bioclimática, y envuelta bioclimática inversa (atomic, range, logistic regression, bioclimatic envelope y negative bioclimatic envelope). Una regla atómica es una conjunción de valores únicos de las variables ambientales. Un ejemplo de regla atómica en lenguaje natural sería: Si la temperatura media anual es igual a 19.5ºC y la precipitación anual es igual a 380mm, entonces la especie está presente. Una regla de rango es una generalización de la regla Bioclim. Según esta regla, determinadas variables pueden considerarse irrelevantes, y tener cualquier valor. La regla de rango permite negación, para aplicarse fuera de determinados rangos. Por ejemplo: Si la temperatura media anual no es mayor de 23ºC y es menor o igual que 26ºC, entonces la especie está presente. Una regresión logística es una forma de ecuación de regresión en la que el resultado es transformado en una probabilidad. La regla de regresión logística es una adaptación de los modelos de regresión logística para funcionar como reglas. Resulta útil cuando la especie tiene una respuesta ecológica basada en los gradientes ambientales. La regresión logística proporciona la probabilidad p de que una regla sea aplicada. p se calcula según la expresión p=1/(1-e-y), donde y es la suma de la ecuación lineal que representa la regla a+(a*a)+b+(b*b)+...+n+(n*n), donde a...n representa las variables y sus coeficientes. Si p es mayor de 0.75, entonces la regla logística predice presencia. Una regla de envuelta bioclimática es una conjunción de rangos de todas las variables tenidas en consideración. Un ejemplo de regla de envuelta bioclimática sería: Si la temperatura media anual está entre 16ºC y 19ºC, y la precipitación anual está entre 300 y 410mm, entonces la especie está presente. Una regla de envuelta bioclimática negativa es la inversa a la anterior, definiendo rangos de ausencia. En GARP, una población es un conjunto de reglas individuales para predecir la presencia o ausencia de la especie en una celda. Las reglas están compuestas de cromosomas, cada uno de los cuales codifica los coeficientes de las distintas variables. El éxito de cada individuo de la población se mide en cada iteración (generación) según la significación estadística de la regla prediciendo presencia y ausencia, determinado su éxito reproductivo. La siguiente generación está formada por los individuos (reglas) con mayor éxito reproductivo. Este procedimiento se repite haste que un criterio de parada es alcanzado. GARP almacena en un archivo las mejores reglas, que forman el conjunto de soluciones al problema.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 38

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Estas reglas son aplicadas a cada celda del territorio (una sola regla o muchas por celda), seleccionando la de mayor perspectiva de éxito para determinar la presencia o ausencia de la especie. El resultado es un mapa que indica las áreas apropiadas (presencia), no apropiadas (ausencia) y sin datos (sus valores quedaron fuera de las reglas seleccionadas). Al repetir el proceso en múltiples ocasiones, se observa que los resultados nunca coinciden; el sistema no es determinista, y cada modelo es diferente (porque GARP introduce mutación estocástica de generación a generación). Un único modelo por tanto suele ser poco significativo, y el tratamiento más actual es el ensamblado de múltiples resultados (sumándolos) para formar un mapa predictivo. En la versión que utilizamos de OpenModeller, hay cuatro implementaciones distintas de este algoritmo, dos de ellas basadas en el software Desktop Garp, y otras dos, nueva implementaciones del código exclusivas para OpenModeller, que consigue resultados mejor ajustados. GARP (single run) – DesktopGARP implementation: genera una única simulación binaria de distribución, que solo indica presencia y ausencia, aunque en ocasiones podemos encontrar áreas sin predicción, si sus valores quedaron fuera del rango de la regla que define el modelo. Los parámetros de esta implementación son: ● Max generations: valor entero, en el rango [1, oo], con un valor típico de 400. Determina el número de generaciones a lo largo de las cuales evoluciona la regla de clasificación. ●

Convergence limit: valor real, en el rango [0, 1], con un valor típico de 0.01. Es la medida de cuanto cambia un set de reglas de una iteración a otra. En el tiempo el límite de convergencia decrece, y decrece más rápido cuando el set de reglas se vuelve estable. Es necesario configurarlo con cuidado, porque pude ocasionar largos tiempos de cálculo y modelos muy sobreajustados.

●

Population size: valor entero, en el rango [1, 500], con un valor típico de 50. Determina el tamaño del conjunto de reglas.

●

Resamples: valor entero, en el rango [1, 100000], con un valor típico de 2500.

GARP (single run) – OpenModeller implementation: Igual a la anterior, pero utiliza un algoritmo mejorado para obtener resultados más precisos. GARP with best subsets – DesktopGARP implementation: es una ampliación de la anterior, preparada para generar múltiples modelos binarios, evaluarlos y unirlos para formar un único modelo continuo. Comparte con la anterior los cuatro parámetros de la misma, configurados por defecto con los mismos valores, y añade varios más: ●

Training proportion: valor real, en el rango [0, 1], con un valor típico de 0.5. Determina la proporción de registros de presencia que se utilizará para entrenar el modelo, mientras la proporción restante se utiliza para evaluar cada iteración (run).

●

Total runs: valor entero, en el rango [1, 10000], y un valor típico de 10. Determina el número total de modelos binarios que se generan durante el proceso

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 39

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

de modelado. ●

Hard Omission Threshold: valor real, en el rango [0, 100]. Determina el error de omisión máximo aceptable (error que consiste en dejar fuera del área de presencia del modelo un registro de presencia de la especie) en porcentaje del número de registros de presencia. Si se establece en 100, se permiten modelos con cualquier valor de omisión.

●

Models Under Omissión Threshold: valor entero, en el rango [0, 10000], con un valor típico de 20. Determina el número de modelos que deben quedar por debajo del Hard Omission Threshold. Debería ser inferior al número total de “Runs”, pero no está configurado de este modo por defecto (desconozco el motivo, tal vez un error).

●

Commission Threshold: valor real, en el rango [0, 100], y valor típico 0.5. Determina el porcentaje de modelos de distribución que deben retenerse teniendo en cuenta el error de comisión (error de incluir en el área de presencia del modelos puntos con ausencia de la especie. En este caso se sustituyen ausencia por puntos aleatorios).

●

Commission Sample Size: valor entero, en el rango [1, oo], con valor típico 10000. Determina el tamaño de la muestra utilizada para calcular el error de comisión.

GARP with best subsets – OpenModeller implementation: Igual a la anterior, pero utiliza un algoritmo mejorado para obtener resultados más precisos. 4.3.6 Support Vector Machines Las Support Vector Machines son métodos de aprendizaje artificial utilizados para clasificación. Teniendo en cuenta el conjunto de datos de entrada como dos vectores en un espacio de n-dimensiones, una SVM tratará de construir un hiperplano (objeto geométrico mínimo para dividir un espacio en dos) que maximice las diferencias entre los dos vectores. Los objetos se clasifican según su distancia al hiperplano. Cuanto mayor es la distancia entre el objeto y el hiperplano, más robusta es la clasificación.

Figura 7: Representación gráfica de un hiperplano separando dos conjuntos de casos, que resume el funcionamiento del algoritmo SVM

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 40

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

4.4 Algunas consideraciones sobre los resultados Aunque no estés familiarizado con el mundo de los modelos de distribución, probablemente habrás llegado a alguna de estas conclusiones (y tal vez, algunas más...¡compártelas!). ●

Hay muchos modos diferentes de calcular modelos de distribución.

●

Distintos algoritmos producen resultados diferentes, a pesar de utilizar las mismas variables y puntos de presencia para calibrar los modelos.

●

Los modelos resultantes pueden ser binarios (booleanos, de presencia-ausencia), o continuos, expresando la idoneidad del hábitat según un gradiente.

●

Los modelos predicen áreas de presencia más extensas de lo que podemos inferir observando los puntos de presencia de la especie.

Y supongo que también algunas cuestiones: ●

¿Cuáles son más útiles, los modelos binarios, o los continuos?

●

¿Qué representan los gradientes de los modelos continuos?

●

¿Por qué son tan diferentes unos modelos de otros?

●

¿Cuál de los algoritmos proporciona mejores resultados?

●

Y todo esto, aparte de hacer mapitas de colores, ¿tiene alguna utilidad REAL?

En las siguientes secciones trabajaremos distintos ejemplos prácticos a través de los que profundizaremos en estas cuestiones.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 41

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

SECCIÓN 5: EVALUACIÓN DE MODELOS DE DISTRIBUCIÓN (Tiempo estimado: 1 hora y 30 minutos) La evaluación estadística de los modelos de distribución es una de las cuestiones que más controversia causa en el debate relativo a la calidad de los modelos de distribución. Evaluar un modelo de distribución es medir su calidad, comprobar hasta que punto se ajusta a datos reales, para de algún modo determinar hasta que punto es capaz de representar con fiabilidad información desconocida. Actualmente no hay un consenso sobre que técnica de evaluación es la más adecuada, y cada autor utiliza la que prefiere. En esta sección se pretende proporcionar unos conocimientos previos y herramientas suficientes para afrontar una evaluación de modelos de distribución lo más correcta posible.

5.1 Evaluación de modelos binarios (presencia/ausencia) 5.1.1 Cálculo de la sensibilidad, error de omisión y error de comisión Imagina un modelo de distribución binario, como Bioclim, y distribuidos sobre el mismo, una serie de registros de presencia con los que hemos calibrado el modelo. En la Figura 8 puedes observar un ejemplo.

Figura 8: Modelo de distribución binario y presencias de la especie modelada

En este caso el modelo ha clasificado correctamente un 60% de los casos (3 aciertos), por lo que se dice que tiene una sensibilidad igual a 0.6 (presencias acertadas / número total de presencias). Sin embargo, dos registros de presencia han quedado fuera del área de hábitat idóneo según el modelo; han sido omitidos, por lo que este tipo de error se denomina error de omisión (falso negativo). Los errores de omisión pueden deberse a: ●

Incapacidad de algoritmo para discriminar el hábitat idóneo de la especie.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 42

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

●

Las variables clave en la distribución del organismo no están representadas en las que hemos utilizado para evaluar el modelo.

●

La identificación del taxón es errónea.

●

La georreferencia de la localidad es errónea.

●

El registro de presencia representa un individuo fuera de su hábitat idóneo (población sumidero).

El error de omisión es relativamente grave, y es necesario estudiar las causas concretas que lo motivan. Utilizar la sensibilidad parece un modo sencillo de evaluar un modelo, pero a continuación veremos los problemas que origina. Observa la Figura 9: ambos “modelos” tienen un 0% de error de omisión (sensibilidad = 1)...¿cuál de ellos representa mejor la distribución no conocida de la especie?.

Figura 9: Dos ejemplos de modelos "malos"; sobrepredicción (comisión) en a), y sobreajuste (omisión) en b).

El modelo a) probablemente sobreestima el área de presencia potencial de la especie (si imaginas que trabajamos sobre un área de cientos de kilómetros de lado, claro), y está incluyendo zonas en las que no se ha reportado la presencia de la especie, y que podrían no ser apropiadas para la misma. Es lo que se denomina error de comisión, que podría o no ser un error (aún no lo sabemos), y tiene diferentes causas (las dos primeras son errores “suaves”, mientras que la última supone un error real): ●

El área cumple con las condiciones idóneas para la especie, pero no ha sido muestreada. La especie tal vez está allí.

●

El área es apropiada para la especie, pero circunstancias biogeográficas han impedido que la especie la ocupe, o se ha extinguido en esa localidad.

●

El área realmente no es apropiada para la persistencia de la especie.

El modelo b), en cambio, está “sobreajustado” (solo predice presencia en el entorno de los puntos de calibrado), y probablemente está omitiendo lugares en los que es posible la presencia de la especie (error de omisión, otra vez, pero no podemos medirlo), por lo que carece de utilidad para hacer cualquier tipo de predicción.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 43

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Caso práctico: evaluación de modelos binarios midiendo la sensibilidad

Vamos a trabajar sobre esta cuestión con un caso real. Utilizaremos tres modelos binarios (Bioclim y las dos implementaciones de GARP “single run”) de Linaria nigricans, a alta resolución (90 metros) para calcular la sensibilidad de cada uno de ellos, y compararlos entre sí. 1. Prepara un nuevo experimento en OpenModeller, utilizando el conjunto de capas ALMERIA_90m tanto en Creation Layers como en Projetion Layers. Nombra el experimento como L_nigricans_binarios_90m. Marca el algoritmo BIOCLIM, y las dos implementaciones “single run” de GARP. Ten paciencia, porque con la resolución de estas capas OpenModeller tarda un poco más en hacer los cálculos. 2. Ahora observa los resultados en la pestaña Output Map; mira para cada modelo, usando el zoom, cuantos puntos de presencia quedan dentro y fuera del área clasificada como apta. Observarás que BIOCLIM incluye todas las presencias en el área apta (sensibilidad = 1). En el modelo GARP (DesktopGARP) podrás encontrar algunas presencias fuera del área apta (como el modelo no es determinista, la cifra puede variar; yo encontré 3). Probablemente la implementación OpenModeller de GARP deje fuera unas cuantas más (13). Esta observación visual de los resultados solo tiene un carácter demostrativo, porque lo cierto es que OpenModeller mide la sensibilidad de cada modelo (aunque la llama accuracy, y la expresa en porcentaje). En la pestaña Report de cada modelo se encuentra el apartado Output statistics, que además de la sensibilidad, indica el error de omisión (1 – sensibilidad), el porcentaje de celdas clasificadas como aptas, y el total de celdas clasificadas. Los resultados de los tres modelos deben ser parecidos a los que presenta la Tabla 1: modelo

sensibilidad omisión celdas aptas (%)

BIOCLIM

1

0

27.32

GARP desktop

0.98

0.02

46.23

GARP openmodeller

0.88

0.12

19.82

Tabla 1: Resultados de la evaluación de modelos binarios presentados por OpenModeller

¡ATENCIÓN!: Mientras trabajamos esta cuestión, aprovecharemos el tiempo para que OpenModeller genere varios modelos continuos de alta resolución para el apartado Evaluación de Modelos Continuos. Prepara un nuevo experimento con el nombre L_nigricans_continuos_90m, y marca para ejecución los siguientes algoritmos: Climate Space Model, Envelope Score, Environmental distance (las cuatro variantes), los GARP que no son “single run”, y SVM (solo el primero de ellos). Pon el experimento en ejecución; tardará “un rato largo”. 5.1.2 Partición aleatoria de datos En el ejemplo anterior hemos utilizado para evaluar el modelo los mismos registros de presencia con los que lo calibramos; esta técnica de evaluación, que se llama resustitución, siempre proporciona resultados muy optimistas que nos pueden llevar Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 44

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

a error (como ejemplo, el valor de sensibilidad del modelo b en la Figura 9, o los valores de BIOCLIM y GARP desktop en el caso de estudio). La resustitución no se considera una técnica válida para evaluar modelos, y la alternativa habitual es la partición aleatoria de datos (data-splitting). El conjunto de presencias se divide, al azar, en dos grupos, uno de calibrado y otro de evaluación, con una proporción variable entre 0.5 y 0.8, aproximadamente, y según el número de presencias disponibles (cuanto mayor es el número de presencias, más grande puede ser el conjunto de evaluación). En la Figura 10 siguiente se representa un ejemplo de dos modelos calibrados con los registros de presencia (verde) que se evaluarán con los registros de evaluación (azul). Los registros de evaluación no se han utilizado para calibrar el modelo.

Figura 10: Modelos a) y b) mostrando presencias de calibrado (verde) y presencias de evaluación (azul)

Ahora los valores de sensibilidad son 1 para el modelo a, y 0 para el modelo b. En este caso, la técnica de partición de datos ha resultado útil para determinar que el modelo sobreajustado no tiene ningún poder predictivo, pero aún nos sigue diciendo que el modelo a es “muy bueno” para predecir la presencia de la especie, algo bastante sospechoso. Caso práctico: evaluación de modelos binarios utilizando partición aleatoria de datos

A continuación evaluaremos los modelos BIOCLIM y GARP (ambas implementaciones), usando el método de partición aleatoria de datos. 1. Creación de un nuevo mapset para almacenar los modelos: Inicia GRASS, selecciona la localidad ALMERIA_latlong, y crea un nuevo conjunto de mapas llamado MODELOS, utilizando la opción Create new mapset in selected location. Inicia la sesión en el nuevo mapset. 2. Importación de los modelos: GIS Manager > File > Import raster map > Multiple formats using GDAL, para abrir el módulo r.in.gdal: ○

Pestaña Options: ■

○

Override projection = marcado

Pestaña Required: ■

Raster file to be imported = C:/CURSO_MODELOS/RESULTADOS/mo-

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 45

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

delOutputs/L_nigricans_binarios_90m/Bioclim/Linaria_nigricans_projection.tif ■ ●

Name for output raster map = BIOCLI

Run ○

Pestaña Required: ■

Raster file to be imported = C:/CURSO_MODELOS/RESULTADOS/modelOutputs/L_nigricans_binarios_90m/GARP_(single_run)__DesktopGARP_implementation/Linaria_nigricans_projection.tif

■

Name for output raster map = GARPDK

○

Run

○

Pestaña Required: ■

Raster file to be imported = C:/CURSO_MODELOS/RESULTADOS/modelOutputs/L_nigricans_binarios_90m/GARP_(single_run)__new_openModeller_implementation/Linaria_nigricans_projection.tif

■

Name for output raster map = GARPOM

○

Run

○

Cierra r.in.gdal

3. Adaptación de la región de trabajo a los modelos: teclea y ejecuta en la línea de comandos: g.region rast=BIOCLI 4. El modelo Bioclim importado a GRASS presenta tres valores: 0 para zonas no aptas, 50 para zonas marginalmente aptas, y 100 para zonas completamente aptas. No distinguiremos entre los valores 50 y 100. Los modelos binarios GARP importados desde OpenModeller tienen valores 0 en las áreas clasificadas como aptas, y valores 100 en las áreas clasificadas como no aptas. Ahora recodificaremos los valores de presencia a 1: GIS Manager > Raster > Change categorie values and labels > Recode interactively, para abrir la ventana Interactive rules entry. ○

Map to recode = BIOCLI

○

Recoded map = LNBIOCLI

○

En la ventana escribe: 100:100:1:1 50:50:1:1 0:0:0:0

○

Apply

○

Aceptar en la ventana de error (ignórala, todo ha salido bien)

○

Map to recode = GARPDK

○

Recoded map = LNGARPDK

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 46

Modelos de distribución y estudios de flora amenazada ○

En la ventana escribe: 100:100:1:1 0:0:0:0

○

Apply

○

Aceptar

○

Map to recode = GARPOM

○

Recoded map = LNGARPOM

○

En la ventana mantén los valores anteriores

○

Apply

○

Ok

Guión de prácticas

5. Borra los mapas que ya no nos sirven, tecleando y ejecutando en la línea de comandos: g.remove BIOCLI,GARPDK,GARPOM 6. Ahora importaremos los puntos de evaluación a GRASS: ○

En la carpeta C:\CURSO_MODELOS\GEODATOS\PRESENCIAS está la hoja de cálculo PRESENCIAS.ods; ábrela. En la hoja LNIGRICANS debes tener diferenciados los registros de calibrado y evaluación. ■

Crea una nueva hoja (puedes nombrarla EVALUACION_Ln), pega los registros de evaluación en ella, elimina las dos últimas columnas, y a las restantes ponles de nombre de campo ID, ESPECIE, LONGITUD y LATITUD respectivamente.

■

Crea una nueva columna a la derecha, nómbrala CASO, y rellena las celdas con el valor de texto “presencia”.

■

Crea tres columnas nuevas con los nombres LNBIOCLI, LNGARPDK y LNGARPOM, y rellena todas las celdas con valores 0,0.

■

Guarda el resultado tal cuál está (Archivo > Guardar)

■

Guarda el resultado con el nombre EVALUACION_Ln, en formato dBASE (.dbf), en la carpeta C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\MO DELOS\dbf.

○

Cierra el fichero, vuelve a abrirlo (menú Archivo, Documentos recientes), y en el nombre del primer campo, sustituye el 2 por un 0. Ahora guárdalo y cierra Calc.

○

Para transformar los puntos de evaluación en un archivo vectorial de GRASS: GIS Manager > Vector > Generate points > Generate points from database, para abrir el módulo v.in.db. ■

Input table name = EVALUACION_Ln

■

Driver name = dbf

■

Name of column containing x coordinate = LONGITUD

■

Name of column containing y coordinate = LATITUD

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 47

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

■

Must refer to an integer column = ID

■

Name for output vector map =PUNTOS_EVALUACION_LINARIA

■

Run

■

Cierra v.in.db.

7. Ahora utilizaremos este archivo vectorial para “preguntar” a los modelos cuáles de estos puntos quedan dentro y fuera de las áreas clasificadas como aptas. ○

○

○

GIS Manager > Vector > Update point attributes from raster > Sample raster map at point locations, para abrir el módulo v.what.rast. ■

Name of input vector ... = PUNTOS_EVALUACION_LINARIA

■

Name of existing raster ... = LNBIOCLI

■

Column name = LNBIOCLI

Run ■

Name of existing raster ... = LNGARPDK

■

Column name = LNGARPDK

Run ■

Name of existing raster ... = LNGARPOM

■

Column name = LNGARPOM

○

Run

○

Cierra v.what.rast.

8. Abre con Calc el fichero PUNTOS_EVALUACION_LINARIA.dbf, que está en C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\ MODELOS\dbf. En él tienes los datos necesarios para calcular la sensibilidad de cada modelo (dividiendo la sumatoria de cada columna por 75, que es el número de registros de evaluación). ¡Después del cálculo recuerda cerrar el fichero sin guardar los cambios!. Los resultados deberían ser parecidos a los mostrados en la Tabla 2 (se comparan con los obtenidos al evaluar los modelos solo con los registros de calibrado): sensibilidad sensibilidad celdas aptas (%) (evaluación) (calibrado)

modelo

aciertos

BIOCLIM

63

0.84

1

23.66

GARP desktop

70

0.93

0.99

38.94

GARP openmodeller

62

0.83

0.86

13.40

Tabla 2: Resultados de la evaluación de modelos binarios utilizando registros de evaluación. Observa que BIOCLIM es para el que más han variado los valores de sensibilidad al usar puntos de evaluación.

Evaluando los modelos utilizando la técnica de resustitución, llegábamos a la conclusión de que BIOCLIM era “mejor modelo”, seguido de cerca por GARP desktop, y de lejos por GARP openmodeller. La utilización de un conjunto de puntos de evaluación Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 48

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

independiente muestra que estos resultados eran erróneos, y en realidad, en cuanto a sensibilidad, los tres modelos tienen un desempeño muy similar. Sin embargo, observa los porcentajes de celdas aptas; GARP openmodeller presenta la proporción más baja de celdas clasificadas como aptas, lo que en este caso lo convierte en el modelo con mejor proporción celdas aptas / presencias detectadas. En este caso, sería nuestra elección de “mejor modelo” entre los tres disponibles. 5.1.3 registros de ausencia y matriz de confusión En los ejemplos anteriores, “sabemos” que los modelos están identificando como aptas áreas que posiblemente no lo son. ¿de qué modo podemos evaluar esta circunstancia?. Una posibilidad es disponer de registros de ausencia de la especie. Con la inclusión de las ausencias, ahora aparece un nuevo tipo de acierto: dejar fuera del modelo un registro de ausencia; y ya es posible medir el error de comisión: incluir registros de ausencia dentro del área de presencia potencial del modelo(error de comisión). En este momento podemos resumir los posibles tipos de aciertos y errores de un modelo binario en lo que se denomina matriz de confusión (ver Tabla 3): Datos reales (registros de presencia y ausencia)

Datos simulados (modelo de distribución)

presencia

ausencia

presencia

A

B

ausencia

C

D

Tabla 3: Matriz de confusión, que resume los errores (en gris) y los aciertos posibles al evaluar un modelo de distribución usando presencias y ausencias

Nota que los términos que indican algún tipo de error están marcados en gris En la matriz de confusión los términos son: ●

A: presencias correctamente clasificadas por el modelo.

●

B: ausencias erróneamente clasificadas por el modelo como presencias; falso positivo, o error de comisión.

●

C: presencias erróneamente clasificadas por el modelo como ausencias; falso negativo o error de omisión.

●

D: ausencias correctamente clasificadas por el modelo.

●

N: suma de todos lo casos

NOTA: En la matriz de confusión, además de valores de conteo, se pueden utilizar proporciones. ¿Qué utilidad tiene esta matriz de confusión?: A partir de la matriz de confusión pueden calcularse diversas medidas de precisión de los modelos. Algunas, como la sensibilidad (A/(A+C)) o la especificidad (capacidad de discriminar ausencias; D/ (B+D)), solo hacen uso de dos de los términos de la matriz, mientras que hay algunas que utilizan los cuatro términos. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 49

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Una de las utilizadas con mayor frecuencia en la evaluación de modelos de distribución es Kappa, que se calcula según la siguiente expresión: K = [(a+d)-(((a+c)*(a+b)+(b+d)*(c+d))/N)] / [N-(((a+c)*(a+b)+(b+d)*(c+d))/N)] Kappa representa la probabilidad de que el modelo se ajuste a la muestra de evaluación más allá de lo esperado por azar. Sus valores están en el rango [-1, 1] y pueden interpretarse según el siguiente baremo: ●

K < 0.4

- escaso

●

0.4 < K < 0.75

- bueno

●

K > 0.75

- excelente

En esta página (http://www.glue.umd.edu/~dchoy/thesis/Kappa/#) se ilustra el concepto de Kappa de modo interactivo. Ahora observa la Figura 11y la Tabla 4; te ayudarán a aclarar ideas:

Figura 11: Tres modelos generados con distintos algoritmos mostrando registros de calibrado, evaluación y ausencias

Teniendo en cuenta el ejemplo de la figura, estas serían las matrices de confusión de cada uno de los modelos: Datos reales (registros de presencia y ausencia) modelo a presencia Datos simulados (modelo de distribución) ausencia

modelo b

modelo c

pres.

aus.

pres.

aus.

pres.

aus.

5

5

0

0

3

1

0

0

5

5

2

4

Tabla 4: Matrices de confusión de los modelos de la Figura 11; pres. = presencia; aus. = ausencia.

En la Tabla 5 se observan los valores de sensibilidad, especificidad y kappa de los tres modelos de la figura anterior. Aunque Kappa ha sido utilizado frecuentemente como medida de precisión para modelos de distribución, y algunos autores han determinado que su desempeño es satisfactorio, se sabe que es sensible al tamaño de muestra (peores resultados con tamaños de muestra pequeños), y proporciona resultados difíciles de interpretar cuando los tamaños de muestra no están balanceados.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 50

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

modelo sensibilidad especificidad kappa a

1

0

0

b

0

1

0

c

0.6

0.8

0.4

Tabla 5: Valores de evaluación de los modelos a), b) y c).

Caso práctico: evaluación de modelos binarios utilizando partición aleatoria de datos y registros de ausencia

A continuación vamos a evaluar los tres modelos que creamos en el primer caso práctico de evaluación, pero añadiendo cosas nuevas: registros de ausencia, y las medidas de precisión especificidad y kappa. 1. En la carpeta C:\CURSO_MODELOS\GEODATOS\PRESENCIAS tienes la hoja de cálculo AUSENCIAS.ods; ábrela con Calc. 2. En la carpeta C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\ MODELOS\dbf tienes el fichero EVALUACION_Ln.dbf; ábrelo también. 3. Copia todas las filas de datos de AUSENCIAS.ods (no copies los nombres de columna), y pégalas tras el último registro de la tabla EVALUACION_Ln., de modo que las columnas casen entre sí. Continúa los valores de ID hasta el último registro, y los ceros de las tres últimas columnas. Ahora tienes en una misma tabla los registros de presencia y ausencia de Linaria nigricans. Cierra AUSENCIAS, y guarda los cambios en EVALUACION_Ln. 4. Ahora crearemos el archivo vectorial con los puntos de presencia y ausencia: ○

GIS Manager > Vector > Generate points > Generate points from database, para abrir el módulo v.in.db. ■

Input table name = EVALUACION_Ln

■

Name of column containing x coordinate = LONGITUD

■

Name of column containing y coordinate = LATITUD

■

Must refer to an integer column = ID

■

Name for output vector map =PUNTOS_EVALUACION_LINARIA

■

Allow overwrite = activado

■

Run

■

Cierra v.in.db

5. Ahora utilizaremos este archivo vectorial para “preguntar” a los modelos cuáles de estos puntos quedan dentro y fuera de las áreas clasificadas como aptas, tal y como hicimos en el caso práctico anterior. ○

GIS Manager > Vector > Update point attributes from raster > Sample raster map at point locations, para abrir el módulo v.what.rast. ■

Name of input vector ... = PUNTOS_EVALUACION_LINARIA

■

Name of existing raster ... = LNBIOCLI

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 51

Modelos de distribución y estudios de flora amenazada ■

Column name = LNBIOCLI

■

Run

■

Name of existing raster ... = LNGARPDK

■

Column name = LNGARPDK

■

Run

■

Name of existing raster ... = LNGARPOM

■

Column name = LNGARPOM

■

Run

■

Cierra v.what.rast.

Guión de prácticas

6. En la carpeta C:\CURSO_MODELOS\EVALUACION tienes la hoja de cálculo KAPPA. Está preparada para introducir los valores de la matriz de confusión de cada uno de los modelos que evaluaremos. Ahora abre el archivo PUNTOS_EVALUACION_LINARIA que está en C:\CURSO_MODELOS\ GEODATOS\GRASSDB\ALMERIA_latlong\MODELOS\dbf. En la columna correspondiente a cada modelo están los valores que necesitas para calcular Kappa. Sigue las siguientes reglas para rellenar las celdas correspondientes en la hoja KAPPA: ○

Un valor 1 en un registro de presencia indica una presencia acertada; el número de presencias acertadas es A en la matriz de confusión.

○

Un valor 0 en un registro de presencia indica una presencia clasificada erróneamente como ausencia (error de omisión, o falso negativo); el número de falsos negativos es C en la matriz de confusión.

○

Un valor 1 en un registro de ausencia indica una ausencia clasificada erróneamente como presencia (error de comisión, o falso positivo); el número de falsos positivos es B en la matriz de confusión.

○

Un valor 0 en un registro de ausencia indica una ausencia correctamente clasificada; el número total es D en la matriz de confusión.

El resultado debe ser similar al que muestra la Figura 12:

Figura 12: Resultados de evaluación de los modelos BIOCLIM, GARPDK y GARPOM utilizando ausencias.

Los modelos LNBIOCLI y LNGARPOM tienen la misma sensibilidad (en este caso), es decir, clasifican presencias con la misma precisión. Sin embargo, LNGARPOM presenta una especificidad superior, es decir, clasifica las ausencias mejor que LNBIOCLI. Por tanto, su índice Kappa es mucho mejor. Tal y como sospechábamos desde el principio, LNGARPOM es el modelo mejor ajustado de los tres.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 52

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Sin embargo hay una cuestión a tener en cuenta, y es importante: el valor de Kappa (y el de cualquier otra medida de precisión) depende de la localización espacial de los registros de presencia y ausencia. Las presencias tienen un carácter inequívoco (salvo error de determinación); la especie ha sido observada, y no hay duda al respecto. Sin embargo, la localización de las ausencias depende del criterio del investigador: pueden determinarse ausencias en áreas adyacentes a las poblaciones conocidas, en áreas alejadas, en la puerta de casa o en mitad del Mediterráneo...y cada uno de estos criterios va a condicionar el valor de kappa. Además un registro de ausencia puede ser realmente una falsa ausencia si se da alguno de los siguientes casos: ●

La planta estaba en el lugar, pero no fue detectada, porque era poco aparente, o aún no había germinado;

●

El hábitat es apropiado, pero una barrera o limitación dispersiva ha impedido que la planta llegue allí;

●

Se ha extinguido recientemente por un suceso ambiental estocástico, etc...

Los modelos de distribución con los que estamos trabajando no utilizan registros de ausencia durante la fase de calibrado, por lo que “diseñar” muestras de ausencia para evaluarlos supone un trabajo extra que no redunda en la precisión de los modelos, y este “diseño” está sujeto a los errores que anteriormente hemos comentado. 5.1.4 registros aleatorios como sustitutos de las ausencias Llegados a este punto, podemos obviar los problemas intrínsecos de las ausencias y utilizarlas teniendo en cuenta las posibles consecuencias, o podemos olvidarnos de ellas, cambiar de concepto y utilizar puntos aleatorios como sustitutos de las ausencias. Los puntos aleatorios se disponen al azar sobre el territorio, y permiten evaluar hasta que punto un modelo puede distinguir entre un patrón de distribución real (representado por los puntos de presencia del conjunto de evaluación) de un patrón de distribución aleatorio (representado por los puntos aleatorios). Los puntos aleatorios tienen la ventaja de que no existen problemas conceptuales ni posibles errores relacionados con ellos; pueden estar en cualquier parte del área de trabajo, y todos son igual de válidos. El reemplazo de registros de ausencia por puntos aleatorios exige redefinir algunos conceptos, como la matriz de confusión descrita en la Tabla 3 (ver Tabla 6). Notarás que el término B, que anteriormente representaba un error (error de comisión), ahora ya no lo representa (aunque tampoco es un acierto). Esto supondrá un cambio en el modo de hacer los cálculos, e interpretar los resultados. Trataremos de aclarar esto a continuación. Datos reales (registros de presencia y puntos aleatorios)

Datos simulados (modelo de distribución)

presencia

aleatorio

presencia

A

B

ausencia

C

D

Tabla 6: Matriz de confusión modificada para trabajar con puntos aleatorios en lugar de ausencias. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 53

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

En la Figura 13 ilustramos un problema que puede derivarse de evaluar un modelo de distribución con dos grupos distintos de puntos aleatorios.

Figura 13: Modelo con presencias de evaluación y dos grupos distintos de puntos aleatorios.

Los valores de sensibilidad, especificidad y kappa del modelo, según lo evaluemos con un grupo de aleatorios u otro quedan reflejados en la Tabla 7: grupo aleatorios sensibilidad especificidad kappa A

0.6

0

-0.4

B

0.6

1

0.6

Tabla 7: Evaluación del modelo de la Figura 13 según dos grupos distintos de puntos aleatorios.

Sí, el mismo modelo ofrece resultados muy diferentes, de hecho, “más” diferentes que los que obtendríamos utilizando ausencias, porque independientemente del criterio utilizado para diseñarlas, estas nunca caerían dentro del área de distribución conocida de la especie. La solución idónea es generar muchos puntos aleatorios (cientos, o miles, según el tamaño de área de trabajo), y computar los valores de evaluación cientos o miles de veces con subgrupos de puntos aleatorios (con el mismo tamaño que la muestra de presencias de evaluación), para ofrecer como resultado una media y una dispersión de los valores de evaluación. Veamos como ejemplo la Figura 14; para el modelo de la imagen, si tomamos los puntos aleatorios en grupos del mismo tamaño que el conjunto de puntos de presencia, los casos posibles son los establecidos en la Figura 15.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 54

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Figura 14: Modelo con presencias de evaluación y puntos aleatorios.

Figura 15: Valores posibles de evaluación que pueden obtenerse en el caso de la Figura 14.

Puedes ver que el valor máximo posible de Kappa es 0.6 (debido al error de omisión del modelo), y el menor, -0.4, queda muy por debajo de lo recomendable. Aunque estos son los casos posibles, la media de Kappa (que estará comprendida entre -0.4 y 0.6) y su desviación estándar estarán determinadas por la relación entre el área de trabajo y el área apta según el modelo, la distribución de los puntos aleatorios, la selección de grupos de aleatorios, y el número de iteraciones de cálculo. Si utilizamos esta técnica para comparar modelos, es necesario que en todos se efectúe el cálculo bajo las mismas condiciones. A continuación realizaremos un ejercicio práctico probando esta aproximación. Caso práctico: evaluación de modelos binarios utilizando partición aleatoria de datos, puntos aleatorios y técnicas de remuestreo

Ya disponemos de registros de presencia para evaluar los modelos (archivo vectorial PUNTOS_EVALUACION_LINARIA), pero contiene presencias y ausencias. Ahora queremos olvidarnos de las ausencias: GIS Manager > Vector > Query with attributes, para abrir v.extract. ●

Name of input vector map = PUNTOS_EVALUACION_LINARIA

●

Name for output vector map = PRESENCIAS

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 55

Modelos de distribución y estudios de flora amenazada ●

WHERE conditions of SQL...= CASO = 'presencia'

●

Run

●

Cierra v.extract.

Guión de prácticas

Ahora necesitamos un fichero vectorial con puntos, muchos puntos aleatorios: 1. Si tienes GRASS abierto en el mapset MODELOS, necesitamos una máscara; para que GRASS busque la máscara (MASK) y la copie desde el mapset PERMANENT al mapset MODELOS, teclea y ejecuta la siguiente expresión: g.copy MASK,MASK 2. Ahora hay que crear puntos aleatorios, pero como GRASS los genera para todo el área de trabajo (incluso fuera de la máscara), tendremos que seguir varios pasos: creamos puntos vectoriales para todo el área de trabajo; convertimos el archivo de puntos vectorial en raster, y después convertimos el archivo raster, compuesto por puntos, en vectorial de nuevo. Es en este último paso donde GRASS si aplica la máscara, y tenemos como resultado puntos aleatorios para toda la provincia de Almería: ○

○

○

Vector > Generate points > Generate random points, para abrir el módulo v.random. ■

Name for output vector map = ALEATORIOS_1

■

Number of points to be created = 22000

■

Run

■

Cierra v.random.

File > Map type conversion > Vector to raster, para abrir el módulo v.to.rast. ■

Name of input vector map = ALEATORIOS_1

■

Name for output raster map = ALEATORIOS_1

■

Source of raster values = val

■

Type = point

■

Run

■

Cierra v.to.rast.

File > Map type conversion > Raster to vector, para abrir el módulo r.to.vect. ■

Name of input raster map = ALEATORIOS_1

■

Name for output raster map = ALEATORIOS

■

Feature type = point

■

Run

■

Cierra r.to.vect.

3. Borra los archivos que ya no necesitas ejecutando secuencialmente: g.remove vect=ALEATORIOS_1 Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 56

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

g.remove ALEATORIOS_1 Ahora dispones de un archivo vectorial con unos 10000 puntos aleatorios cubriendo la provincia de Almería. Los registros de presencia ya los tienes en el archivo vectorial PRESENCIAS. Estas presencias ya llevan asociado su valor para cada uno de los tres modelos, por lo que no es necesario repetir la consulta. Los puntos aleatorios aún no los tienen, pero antes de hacer la consulta, hay que añadirles las columnas correspondientes con los nombres de cada uno de los modelos. En la carpeta C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\MODELOS\ dbf está el archivo ALEATORIOS.dbf; ábrelo con Calc. Elimina las columnas value,N,11,0 y label,C,17, y añade tres columnas llamándolas LNBIOCLI,N,20,2, LNGARPDK,N,20,2 y LNGARPOM,N,20,2. Rellena todas las celdas de estas nuevas columnas con valores 0 (cero), y guarda el fichero. Ahora utilizaremos este archivo vectorial consultar los valores de los puntos aleatorios sobre los modelos. ○

GIS Manager > Vector > Update point attributes from raster > Sample raster map at point locations, para abrir el módulo v.what.rast. ■

Name of input vector ... = ALEATORIOS

■

Name of existing raster ... = LNBIOCLI

■

Column name = LNBIOCLI

■

Run; ten un poco de paciencia, ¡son muchos puntos!;

■

Name of existing raster ... = LNGARPDK

■

Column name = LNGARPDK

■

Run

■

Name of existing raster ... = LNGARPOM

■

Column name = LNGARPOM

■

Run

■

Cierra v.what.rast.

Para preparar los datos para el remuestreo de aleatorios y cálculo de kappa: 1. Abre con Calc los ficheros ALEATORIOS.dbf y PRESENCIAS.dbf, que están en C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\ MODELOS\dbf. 2. En la carpeta C:\CURSO_MODELOS\EVALUACION\KAPPA crea dos archivos de texto vacíos; uno llámalo PRESENCIAS.txt y otro, ALEATORIOS.txt. (¡NO BORRES LOS DEMÁS ARCHIVOS QUE HAY EN LA CARPETA!). 3. Abre PRESENCIAS.txt con Notepad++; pega en él los valores de las columnas LNBIOCLI,N,20,6, LNGARPDK,N,19,2 y LNGARPOM,N,19,2 de la hoja de cálculo PRESENCIAS, y guárdalo. Debe quedar así: 0 1 0

0 1 1

0 1 1

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 57

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

... 4. Abre ALEATORIOS.txt; pega en él los valores (NO LA CABECERA) de las columnas LNBIOCLI,N,20,6, LNGARPDK,N,19,2 y LNGARPOM,N,19,2 de la hoja de cálculo ALEATORIOS, y guárdalo. El aspecto del contenido del archivo debe ser igual al anterior (con distintos valores). Puedes cerrar las hojas de cálculo, por ahora no las necesitaremos. 5. En el mismo directorio (C:\CURSO_MODELOS\EVALUACION\KAPPA) crea un archivo de texto llamado MODELOS.txt, con los nombres de los modelos escritos en una sola columna, tal que así: BIOCLIM GARPDK GARPOM Para calcular Kappa remuestreando los puntos aleatorios, se ha preparado un programa para Octave. Para ejecutar programas escritos en este lenguaje, inicia el programa Octave: Todos los programas > GNU Octave > Octave. Se abrirá una línea de comandos, a la que debemos ordenar moverse a la carpeta en la que están los archivos PRESENCIAS.txt, ALEATORIOS.txt y MODELOS.txt, además del programa de evaluación (KAPPA.m), y dos archivos más que este necesita. Teclea en la línea de comandos de Octave la siguiente expresión y después pulsa ENTER: cd C:\CURSO_MODELOS\EVALUACION\KAPPA Ahora, para ejecutar el programa de evaluación, simplemente teclea el nombre del programa, KAPPA y pulsa ENTER. El programa te pedirá los siguientes parámetros: ●

El nombre de la especie que está representada en los modelos de distribución (Linaria nigricans).

●

Número de iteraciones: es el número de veces que el programa remuestreará los puntos aleatorios para calcular kappa. El máximo recomendado es 10000, que según el ordenador con el que trabajes, puede llevar de 1 a 3 minutos de cálculo. En este caso de estudio 1000 iteraciones son suficientes para distinguir entre modelos sin problemas. Habrá que incrementarlo en caso de modelos con un desempeño similar.

Pulsa ENTER. Tras hacer los primeros cálculos aparecerán en la ventana del programa algunos resultados. AVISO: Si la pantalla de resultados muestra un cartel azul en su base, pulsa la tecla q para desbloquearla. El resultado más intuitivo es el gráfico de cajas, que indica la distribución de los valores Kappa de cada modelo. Una vez observados los índices de evaluación de cada modelo, el programa da la oportunidad de comparar dos de los modelos (los dos de valores más altos, por ejemplo) para comprobar si existen diferencias significativas Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 58

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

entre ellos. Debes introducir sucesivamente, según el programa lo pida, los números correspondiente a los modelos que quieres evaluar, en orden de izquierda a derecha, tal y como están representados en el gráfico. El programa almacena los resultados en C:\CURSO_MODELOS\EVALUACION\ KAPPA, y se componen de: ●

●

Un archivo de texto llamado RESULTADOS_KAPPA.txt, que guarda un resumen de los resultados de la ejecución. Los valores correspondientes a cada modelo están en el mismo orden en el que se introdujeron las columnas en los archivos PRESENCIAS.txt y ALEATORIOS.txt. Los parámetros que muestra este archivo son: ○

Nombre de la especie

○

Nombres de los modelos evaluados

○

Número de presencias

○

Número de puntos aleatorios

○

Iteraciones computadas

○

Kappa medio: este es el valor más importante

○

Desviación estándar de Kappa

○

Kappa máximo: no es muy relevante, porque puede darse en un solo caso

○

Kappa mínimo: tiene el mismo grado de significación que el anterior

○

Sensibilidad: es un valor importante, pro no define cual es el mejor modelo

○

Especificidad: importante; cuanto mayor es, menor es el error de comisión del modelo

○

Desviación estándar de la especificidad

○

Tabla ANOVA comparando todos los modelos entre sí

○

Primer modelo: el primero de los modelos que seleccionaste para una comparación

○

Segundo modelo: el segundo de los modelos que seleccionaste para una comparación

○

Tabla ANOVA comparando los dos modelos anteriores

Una imagen con el gráfico de cajas comparando los valores Kappa de los distintos modelos. El gráfico correspondiente a cada modelo indica, de arriba a abajo: outliers, máximo, 3er cuartil, límite de confianza superior para la mediana, mediana, límite de confianza inferior para la mediana, 1er cuartil, mínimo, outliers.

NOTA SOBRE EL PROGRAMA “KAPPA”: Para evaluar diferencias significativas entre los valores kappa de los distintos modelos, el programa aplica un primer test ANOVA de medidas repetidas, usando todos los modelos, con el algoritmo como factor. Si el estadístico F es menor que el p-valor, podemos concluir que no hay diferencias significativas entre los valores Kappa de los distintos modelos. Como mejor modelo seleccionaríamos aquel con Kappa más alto. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 59

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Una vez examinados los resultados, cierra Octave tecleando exit en la línea de comandos.

5.2 Evaluación de modelos continuos La evaluación de modelos continuos tiene dos posibilidades: 1. Aplicación de un umbral de corte para transformar el modelo continuo en binario (de los valores continuos del modelo, se selecciona uno, según un criterio concreto, por encima del cual se definen las áreas de presencia, y por debajo, las áreas de ausencia), y utilizar alguna de las métricas de evaluación que hemos tratado en el apartado anterior. 2. Evaluación del modelo con sus valores continuos, aplicando nuevas métricas de evaluación, independientes del umbral de corte. La primera cuestión plantea un problema que se tratará posteriormente: el criterio para aplicar un umbral de corte, que encierra algunas complejidades. Sin embargo, una vez aplicado un punto de corte, quedaría acudir a la medida de error Kappa expuesta anteriormente. La segunda posibilidad vamos a tratarla a continuación: 5.2.1 La curva ROC En primer lugar ilustraremos la evaluación de modelos continuos con un ejemplo gráfico, que puedes ver en la Figura 15:

Figura 15: Modelo continuo mostrando registros de presencia (evaluación) y registros de ausencia.

En el modelo se presentan 5 niveles de idoneidad del hábitat (5 representa la idoneidad máxima, y 1 la idoneidad mínima), 8 presencias (del conjunto de datos de Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 60

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

evaluación) y 8 ausencias. En la Figura 16 puedes observar las distintas medidas de evaluación tomadas según los distintos umbrales de corte posibles para transformar este modelo continuo en binario:

Figura 16: Valores de evaluación del modelo de la Figura 15, teniendo en cuenta todos los posibles umbrales de corte.

Observa que en lugar de la especificidad, tal y como hemos calculado hasta ahora, calculamos 1-especificidad, que es la fracción de falsos positivos. La sensibilidad (fracción de aciertos) va disminuyendo progresivamente, según el umbral de corte se hace más alto; cada vez más presencias quedan fuera del área idónea según el modelo. Progresivamente disminuye también el número de ausencias dentro de área idónea (de ahí la disminución de 1-especificidad). Ambas circunstancias provocan el comportamiento de kappa: crece hasta que el umbral de corte es 3, y vuelve a decrecer hasta que el umbral se hace igual a 5. Si en un gráfico representamos los valores de sensibilidad en el eje y, y los de 1-especificidad en el eje x, obtenemos un gráfico como el de la Figura 17:

Figura 17: Representación gráfica de los valores de sensibilidad y 1-especificidad del modelo de la Figura 15.

Si en este gráfico unimos los puntos mediante líneas, y trazamos la diagonal correspondiente obtenemos una curva que delimita un área bajo ella, tal y como puedes observar en la Figura 18: Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 61

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Figura 18: Representación gráfica de la curva ROC correspondiente al modelo de la Figura 15. El área sombreada en rojo es la que indica el valor AUC del modelo.

El resultado es una curva ROC (de “Receiver Operating Characteristic”), cuyo área inferior (AUC = Area Under Curve) indica la probabilidad de que el modelo, ante una presencia y una ausencia escogidas al azar, clasifique con un valor de idoneidad mayor a la presencia que a la ausencia. En este caso, suponiendo que el AUC de esta curva fuera 0.74, para el modelo evaluado, de cada 100 veces que sean seleccionadas al azar una presencia y una ausencia, en 74 ocasiones dará un valor mayor de idoneidad al registro de presencia. El valor de AUC oscila entre 0.5 (valor de la diagonal dibujada en negro) y 1 (vértices de la curva cercanos a la esquina superior izquierda). Los valores cercanos a 0.5 indican que el modelo no puede distinguir presencias de ausencias más allá de lo esperado por azar (“mal modelo”), mientras que valores cercanos a 1 indican una discriminación perfecta entre casos de presencia y ausencia. De modo general se consideran malos los valores por debajo de 0.70 (este varía según los autores). Cuando en lugar de ausencias se utilizan puntos aleatorios, la interpretación de AUC cambia ligeramente, tomándose como la probabilidad de que el modelo discrimine correctamente entre una presencia y un punto aleatorio tomados al azar. Como en el caso de cálculo de Kappa con puntos aleatorios, es necesario generar una cantidad importante de ellos, y realizar técnicas de remuestreo para calcular un valor promedio de AUC y una desviación estándar, que caracterizarán la precisión del modelo evaluado. En este contexto, el máximo valor de AUC será siempre menor que 1, porque los puntos aleatorios pueden caer dentro de áreas con alta idoneidad de hábitat. Si la distribución de la especie cubre una proporción a de las celdas del área de trabajo, entonces el valor máximo de AUC será 1-a/2. Pero el valor a es desconocido (no conocemos la distribución real de la especie), por lo que no podemos saber como de Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 62

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

óptimo es un valor de AUC. Sin embargo, sí que podemos utilizar los valores AUC obtenidos mediante este método para distinguir entre modelos mejor o peor ajustados a los registros de evaluación. Aunque AUC se ha utilizado muy frecuentemente como medida de precisión de modelos de distribución, algunos autores le encuentran ciertas desventajas a este método: ●

La curva ROC resume la precisión del modelo teniendo en cuenta regiones del espacio ROC en los que el usuario raramente trabaja, como los extremos de la curva, en los que las tasas de error son elevadas.

●

Pondera por igual los errores de comisión y omisión.

●

No proporciona información acerca de la distribución espacial de los errores.

●

Las áreas de trabajo amplias resultan en valores de AUC más altos. La precisión de distintos modelos de la misma especie no puede hacerse si las áreas de trabajo difieren en extensión.

●

No pueden compararse modelos de distintas especies.

●

Los registros de ausencia tienen un alto grado de incertidumbre, y tienen gran importancia en el cálculo de AUC.

La crítica más contundente procede de Lobo et al. (2007). En este trabajo se concluye que AUC proporciona información sobre lo generalista o especialista que es una especie teniendo en cuenta el rango de variables predictoras utilizadas, pero no proporciona información sobre el ajuste del modelo. Sin embargo los mismos autores especifican que siempre que se trate de una misma especie y el mismo área de trabajo, AUC permite distinguir modelos por su grado de ajuste con los datos de evaluación. Por el momento no se ha propuesto ninguna alternativa sólida a este método para evaluar modelos de distribución. Caso práctico: Cálculo de la curva ROC para los modelos de Linaria nigricans

A continuación vamos a realizar un proceso que es tan tedioso como valioso: evaluaremos mediante la curva ROC nueve modelos de distribución de alta resolución de la especie Linaria nigricans, para establecer cual de ellos se ajusta mejor a la distribución espacial de la especie. En primer lugar comprobaremos si OpenModeller ha terminado de ejecutar todos los modelos que le encargamos al comienzo de la sección dedicada a la evaluación de modelos. Para evaluar los modelos que hemos generado es necesario seguir estos pasos: 1. Añadir las columnas con los nombres de los nuevos modelos a los archivos vectoriales PRESENCIAS y ALEATORIOS del mapset MODELOS. ○

En C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\ MODELOS\dbf están los archivos ALEATORIOS.dbf y PRESENCIAS.dbf. Ábrelos con Calc, y añádeles a ambos los siguientes nombres de columna: ■

CSM,N,20,2

(Climate Space Model)

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 63

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

■

ES,N,20,2

(Envelope Scores)

■

EDEUC,N,20,2

(Environmental Distance Euclidean)

■

EDCHE,N,20,2

(Environmental Distance Chebyshev)

■

EDMAH,N,20,2

(Environmental Distance Mahalanobis)

■

EDMAN,N,20,2

(Environmental Distance Manhattan)

■

GARPDK,N,20,2

(GARP Desktop Garp Implementation)

■

GARPOM,N,20,2

(GARP OpenModeller Implementation)

■

SVM,N,20,2

(Support Vector Machines)

○

Rellena todas las celdas de las nuevas columnas con el valor 0 (cero).

○

Guarda y cierra ambos ficheros.

2. Importación y consulta sobre los modelos: Como el procedimiento de importación de modelos y consulta de los valores de las presencias y puntos aleatorios ya ha sido trabajado, se ha preparado un guión de GRASS para realizar la tarea automáticamente. Abre GRASS Command Line, cerciorándote de que lo abres en la localidad ALMERIA_latlong, mapset MODELOS. ○

Una vez abierto, muévete al directorio con los guiones; teclea y ejecuta: cd C:\CURSO_MODELOS\GUIONES

○

Para ejecutar el script, teclea y ejecuta (paciencia, tardará un poco; ¡son muchos puntos aleatorios!): sh MODELOS_CONTINUOS.txt

3. Preparación de los archivos de evaluación: ○

Crea dos ficheros de texto en blanco llamados ALEATORIOS.txt y PRESENCIAS.txt en la carpeta CURSO_MODELOS\EVALUACION\AUC.

○

En C:\CURSO_MODELOS\GEODATOS\GRASSDB\ALMERIA_latlong\MODELOS\dbf están los archivos ALEATORIOS.dbf y PRESENCIAS.dbf. Ábrelos con Calc

○

Copia los datos de las columnas correspondientes a los nuevos modelos de los .dbf a los .txt correspondientes (recuerda no incluir las cabeceras).

○

Guarda los cambios en los archivos de texto y ciérralos. También puedes cerrar los abiertos en Calc, sin guardar cambios.

○

En la carpeta CURSO_MODELOS\EVALUACION\AUC crea un archivo de texto en blanco llamado MODELOS.txt. Escribe los nombres de los modelos en una sola columna, en el mismo orden en el que estaban en las tablas de ALEATORIOS y PRESENCIAS, tal que así: CSM ES EDEUC EDCHE EDMAH

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 64

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

EDMAN GARPDK GARPOM SVM 4. Cálculo de AUC: Se ha preparado un programa similar al utilizado para calcular kappa, por lo que para evaluar los modelos debes ejecutarlo: ○

Inicia octave, y muévelo a la carpeta C:\CURSO_MODELOS\EVALUACION\AUC: cd C:\CURSO_MODELOS\EVALUACION\AUC

○

iniciar al programa de cálculo: AUC

○

El programa pide el nombre de la especie y el número de iteraciones que se utilizarán para computar AUC. Este programa es más lento que el usado para calcular Kappa (los cálculos son más complejos), por lo que 1000 iteraciones serán suficientes (en caso de no obtener diferencias significativas entre la pareja de modelos con mayor AUC, deberíamos volver a evaluar incrementando el número de iteraciones).

○

Lo más probable es que la ventana de resultados se bloquee; pulsa q para desbloquearla.

○

Una vez observado el gráfico de cajas, decide la pareja de modelos que mayores valores de AUC presentan (o puede ayudarte del archivo RESULTADOS_AUC.txt, mirando los valores de AUC medio más altos), e introduce sus números (según su orden) para comprobar si existen diferencias significativas entre ellos, tal y como hiciste en la evaluación de modelos binarios. En este caso los modelos 3 y 6 deberían ser los más cercanos.

5. Interpretación de los resultados: El programa escribe un archivo de resultados (RESULTADOS_AUC.txt) que proporciona la siguiente información: ○

Nombre de la especie

○

Lista de los modelos evaluados

○

Número de presencias

○

Número de puntos aleatorios

○

Número de iteraciones

○

Media de las presencias: valor medio de los registros de presencia sobre el modelo

○

Desviación estándar de las presencias: desviación de la media anterior

○

Punto de corte según media y desviación: es una recomendación de umbral de corte para transformar el modelo continuo en binario. Se calcula restando al valor medio de las presencias sobre el modelo su desviación estándar.

○

AUC máxima de cada modelo

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 65

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

○

AUC media: es el valor más importante, que caracteriza la precisión de cada modelo

○

Desviación estándar de AUC

○

AUC mínima

○

Tabla ANOVA comparando entre sí todos los modelos

○

Primer modelo de la comparación

○

Segundo modelo de la comparación

○

Tabla ANOVA comparando los dos modelos anteriores

Has aprendido lo suficiente (y no es poco) para evaluar modelos de distribución según marcan los cánones. Si has llegado hasta aquí sin odiarme...¡enhorabuena!.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 66

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

SECCIÓN 6: TRANSFORMACIÓN DE MODELOS CONTINUOS EN BINARIOS (Tiempo estimado: 10 minutos) Aunque los modelos continuos proporcionan información ecológica de gran valor, para diversas aplicaciones prácticas son más apropiados los modelos binarios. Un ejemplo claro de esta necesidad es el diseño de reservas de flora, que exige unos límites geográficos perfectamente delimitados. La transformación de un modelo continuo (gradiente de idoneidad de hábitat o probabilidad de presencia) en un modelo binario (presencia/ausencia) exige la adopción de un criterio, a partir del cual se selecciona un valor umbral, por encima y debajo del cual se consideran, respectivamente regiones aptas (presencia potencial) y no aptas (ausencia potencial) para la especie. Una posibilidad comúnmente utilizada en trabajos de modelización es la adopción de criterios arbitrarios (0.5 es uno de los más utilizados), pero distintos estudios comparativos demuestran que es una mala opción. Existen métodos diferentes, de los cuáles los más utilizados son: ●

Maximización de Kappa

●

Idoneidad media de los registros de presencia

●

Maximización de la suma sensibilidad-especificidad

●

Métodos basados en la curva ROC

La complejidad de cálculo de algunos métodos no garantiza los mejores resultados; un método sencillo que proporciona resultados, que al menos son tan certeros como los obtenidos mediante métodos más complejos, es el segundo de la lista: establecer el umbral de corte en el valor medio de los registros de presencia (hablando siempre de registros de evaluación, no utilizados para calibrar los modelos). Sin embargo la adopción de un criterio concreto depende de la utilidad que se va a dar al modelo. Suponiendo que trabajamos en el diseño de una reserva para una especie amenazada en un área remota sobre la que no podemos hacer prospecciones, es deseable un umbral “generoso”, que garantice una alta sensibilidad (inclusión de muchos de los registros de presencia), para no dejar fuera de la reserva poblaciones potenciales. Trabajando con especies invasoras, se da el mismo caso, siendo necesario pronosticar nuevas áreas de expansión, para lo que utilizaremos un criterio relajado. Un caso contrario puede ser el diseño de una prospección de campo para encontrar poblaciones de una especie amenazada. Como interesa optimizar el tiempo de muestreo, lo indicado es utilizar un criterio restrictivo, para reducir el área total de búsqueda. Similar es la utilización de un modelo de distribución para estimar un tamaño poblacional; si no utilizamos un criterio restrictivo, estaremos sobreestimando el tamaño poblacional. En todo caso, independientemente del criterio, siempre es importante describir el Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 67

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

método que se ha utilizado y justificar su elección. En el fichero de resultados de la evaluación de modelos continuos mediante AUC están los valores de idoneidad promedio de las presencias (apartado MEDIA_PRESENCIAS), y los valores medios menos la desviación estándar (apartado PUNTO_DE_CORTE_MEDIA_menos_DESVIACION). Ambos pueden utilizarse como umbrales de corte, según los trabajos a realizar. En la Figura 19 se muestra un resumen de como funcionan estos dos criterios en los modelos que hemos generado:

Figura 19: Comparación de dos criterios de selección de umbral de corte para transformar modelos continuos en binarios. En gris se destacan los valores del modelo que mejor valor AUC obtuvo durante la evaluación.

Presta atención a los valores de omisión de cada criterio. El criterio media es un 15% más restrictivo que el criterio media - desviación, pero la idoneidad media de los registros omitidos es relativamente alta. Este criterio deja fuera de cualquier modelo al que se aplique un número considerable de presencias reales, que además tienen valores altos de idoneidad. Pero la omisión baja lleva consigo un incremento considerable del área de presencia potencial tal y como puedes ver en la Figura 20.

Figura 20: Transformación del modelo continuo EDMAN en binario según los dos criterios comparados en la Figura 19. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 68

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Como puedes observar, para reducir la omisión del modelo binario un 15% se incrementa el área potencial en 71244 celdas (¡un 861%!). Es crucial ser cauto en la elección del umbral de corte, sobre todo cuando se desconoce la biología propia de la especie. Caso práctico: transformación del mejor modelo continuo de Linaria nigricans en binario

Una vez decidido un valor apropiado para el umbral de corte, con GRASS es rápido y sencillo crear modelos binarios a partir del modelo continuo: ●

Raster > Change category values and labels > Recode interactively: ○

Map to recode = EDMAN

○

Recoded map = EDMAN_bin (el valor de corte según la media de las presencias)

○

En la ventana escribe lo siguiente, sustituyendo n por el valor MEDIA_PRESENCIAS correspondiente al modelo EDMAN en el archivo RESULTADOS_AUC.txt:

n:100:1:1 0:n:0:0 ○

Apply

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 69

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

SECCIÓN 7: APLICACIONES PRÁCTICAS DE LOS MODELOS DE DISTRIBUCIÓN (Tiempo estimado: 1 hora y 10 minutos) En este apartado se ilustrarán aplicaciones prácticas de los modelos de distribución, trabajando sobre varios ejemplos. Por cuestiones de tiempo, en la mayoría de los casos se simplificará el proceso de modelado, omitiendo fases como la evaluación, o utilizando una resolución menor (1000m).

7.1 CARTOGRAFÍA DE POBLACIONES (Tiempo estimado: 20 minutos) Los objetivos de esta aplicación son medir la superficie potencialmente ocupada por la especie Linaria nigricans, y obtener dos cartografías vectoriales que puedan visualizarse sobre imágenes aéreas y/o satelitales: una a escala de detalle (1:10000), y otra a escala de reconocimiento (1:100000). Los modelos de distribución de especies se diseñaron originalmente para obtener cartografías de especies de distribución desconocida. Sin embargo, al trabajar con cartografías generadas mediante estos métodos, es importante conocer y comunicar las limitaciones de los resultados. No todas las especies son susceptibles de ser cartografiadas mediante estos métodos, y esto es especialmente cierto para aquellas que tienen patrones de dispersión estocásticos, o tienen relaciones poco definidas con las variables ambientales utilizadas como predictores. La especie sobre la que se va a trabajar cumple condiciones para ser bien cartografiada, porque está asociada a una componente litológica bien reflejada en las bandas Landsat, y su mecanismo de dispersión la ha llevado allí donde existía hábitat apropiado para soportar nuevas poblaciones. El material de partida es el mejor modelo de alta resolución obtenido para Linaria nigricans en el caso práctico de la sección 5.2 (Evaluación de modelos continuos). Este modelo ya fue transformado en binario (EDMAN_76) utilizando el criterio basado en el valor de idoneidad promedio de las presencias de evaluación en el caso práctico de la sección 6. 7.1.1 Medición del área de ocupación Se hace directamente sobre el modelo binario. ●

Raster > Reports and statistics > Sum area by map and category, para abrir el módulo r.report. ○

Raster map to report on = EDMAN_bin

○

k = marcado (para obtener el área en kilómetros cuadrados)

○

h = marcado (para obtener el área en hectáreas)

○

c = marcado (para obtener el número de celdas ocupadas)

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 70

Modelos de distribución y estudios de flora amenazada ○

Run

○

Cierra r.report.

Guión de prácticas

Los resultados se mostraran en la ventana de salida, en formato texto. Los co rrespondientes al área ocupada tienen el valor 1 en el campo description. 7.1.2 Cartografía de detalle La cartografía de detalle debe tener toda la información proporcionada por el modelo. Para obtener la cartografía de detalle: 1. Para preparar la capa EDMAN_bin para vectorizarla, es necesario en primer lugar transformar los valores 0 (cero) en Null: Raster > Develop map > Null values, para abrir r.null. ○

Raster map for wich edit null file = EDMAN_bin

○

List of cell values to be set to null = 0

○

Run

2. Para vectorizar la capa: File > Map type conversions > Raster to vector, para abrir r.to.vect. ○

Name of input raster map = EDMAN_bin

○

Name of output vector map = Linaria_nigricans_DETALLE

○

Feature type = area

○

Run

7.1.3 Cartografía de reconocimiento La cartografía a escala de reconocimiento debe mostrar un menor grado de detalle. Si visualizas el modelo binario EDMAN_bin verás que tiene muchas celdas aisladas con presencia potencial de la especie. Aunque estas celdas podrían ser candidatas para una prospección, en el caso de una cartografía de reconocimiento aportan demasiada información, y es preferible eliminarlas filtrando el modelo. 1. Es necesario volver a modificar los valores Null del modelo EDMAN_bin, transformándolos en 0 (cero): Raster > Develop map > Null values, para abrir r.null. ○

Raster map for wich edit null file = EDMAN_bin

○

The value to replace the null value = 0

2. Filtrado del modelo binario utilizando un filtro de moda (es el apropiado para capas binarias y cualitativas): Raster > Neighborhood analysis > Moving window, para abrir el módulo r.neighbors. ○

Use circular neighborhood = marcado

○

Name of input raster map = EDMAN_bin

○

Name of output raster map = EDMAN_bin_moda_5

○

Neighborhood operation = mode

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 71

Modelos de distribución y estudios de flora amenazada ○

Neighborhood size = 5

○

Run

○

Cierra r.neighbors.

Guión de prácticas

3. Para vectorizar EDMAN_bin_moda_5, es necesario transformar los valores 0 (cero) en Null de nuevo: Raster > Develop map > Null values, para abrir r.null. ○

Raster map for wich edit null file = EDMAN_bin_moda_5

○

List of cell values to be set to null = 0

○

Run

○

Cierra r.null.

4. Para vectorizar la capa: File > Map type conversions > Raster to vector, para abrir r.to.vect. ○

Smooth corners of area features = marcado

○

Name of input raster map = EDMAN_76_moda_5

○

Name of output vector map = Linaria_nigricans_RECONOCIMIENTO

○

Feature type = area

○

Run

○

Cierra r.to.vect.

5. Para preparar una visualización sobre una imagen Landsat en color verdadero, es necesario hacer la composición de color “fundiendo” varias bandas. En GIS Manager pulsa el icono a la derecha de abrir nueva capa raster (se llama Add RGB or HIS layer). Selecciónalo en la ventana de mapas activos, y en la configuración de mapas activos establece la siguiente configuración (recuerda que las bandas Landsat están en el mapset PERMANENT): ○

red = LSAT_3

○

green = LSAT_2

○

blue = LSAT_1

6. Añade a la vista la capa Linaria_nigricans_RECONOCIMIENTO, y juega con los colores de línea, relleno (fill areas) y transparencia, para conseguir una buena visualización de la cartografía. 7.1.4 Comparación con una cartografía real En la carpeta C:\CURSO_MODELOS\GEODATOS\PRESENCIAS\Linaria_nigricans dispones del archivo vectorial LINARIA_NIGRICANS_CARTOGRAFIA.shp. Es una cartografía de detalle de las poblaciones de Linaria nigricans. Impórtalo a GRASS: File > Import > Import vector map > Multiple formats using OGR, para abrir v.in.ogr. ●

Pestaña Options: ○

Override dataset projection = marcado

○

OGR datasource name = C:/CURSO_MODELOS/GEODATOS/PRESENCIAS/Linaria_nigricans/LINARIA_NIGRICANS_CARTOGRAFIA.shp

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 72

Modelos de distribución y estudios de flora amenazada ○

Name of output vector map = Linaria_nigricans_REAL

○

Type = boundary

○

Run

○

Cierra v.in.ogr.

Guión de prácticas

Visualiza y compara la cartografía real con la de reconocimiento y detalle que has generado a partir del modelo de distribución. La Figura 21 ofrece una comparación entre las grandes poblaciones conocidas de la planta.

Figura 21: Comparación entre las cartografías obtenidas usando un modelo de distribución como base, y una cartografía real de gran precisión con la distribución conocida de Linaria nigricans.

7.1.5 Discusión Ya conocemos la relevancia que tiene el umbral de corte en el área de presencia final de un modelo de distribución binario. La resolución del modelo también tiene un papel importante en el área final, y cuanto mayor sea la resolución, mejor será la medida final de área. Estas medidas de área siempre deben tomarse como área de presencia potencial. Para tomarlas como realistas habría que realizar verificaciones en campo. Una posibilidad para mejorar la estimación de área es utilizar para medirla solo aquellos polígonos en los que hay presencia registrada de la planta. Con la resolución actual de trabajo (unos 90 metros), el resultado de la cartografía de detalle aparecerá muy cuadriculado. Para obtener cartografías de detalle finas, es necesario recurrir a resoluciones mayores en las variables predictoras. Al comparar los resultados con la cartografía real, debes tener en cuenta que para generar la real, varios equipos de personas han trabajado en campo, durante varios años consecutivos y utilizando GPS e imágenes aéreas de alta resolución como base topográfica. El resultado que se obtiene con la cartografía de reconocimiento de Linaria nigricans es relativamente bueno, porque esta especie tiene características apropiadas para modelar su distribución, pero ojo, esto no siempre va a ser así.

7.2 BÚSQUEDA DE NUEVAS POBLACIONES DE Linaria nigricans (Tiempo estimado: 35 minutos) El objetivo de esta aplicación es obtener un mapa útil para programar prospecciones de campo encaminadas a localizar poblaciones desconocidas de Linaria nigricans. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 73

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Esta aplicación podría considerarse relacionada con la anterior, pero en la cartografía el objetivo es representar la distribución más restringida posible, mientras que para buscar nuevas poblaciones hay que flexibilizar los criterios para tener en cuenta áreas en las que, por ahora, no hay registros de presencia de la especie. El primer criterio para considerar una nueva población es la distancia a otras poblaciones conocidas. Puede seleccionarse una distancia arbitrariamente, o utilizar un criterio razonado, pensando en las capacidades dispersivas de la especie, o en la distancia entre las poblaciones actualmente conocidas. En este caso tomaremos una distancia mínima de 10 kilómetros, porque Linaria nigricans aparentemente dispersa ocasionalmente a larga distancia, y presenta un área de distribución amplia. Para definir esta distancia utilizaremos como referencia la cartografía real de la planta (Linaria_nigricans_REAL) utilizada en el ejercicio anterior. El segundo criterio es de idoneidad del hábitat, y para definirlo vamos a utilizar un ensamblado de modelos. Se transformarán todos los modelos de Linaria nigricans en binarios, utilizando el umbral de corte definido por el criterio más relajado (media menos desviación estándar). Se sumarán todos los modelos, y se tomarán como localidades potenciales aquellas en las que la mayoría absoluta de los modelos predice presencia. 7.2.1 Procedimiento de búsqueda de nuevas poblaciones Los pasos serán: calcular las distancias a las poblaciones conocidas, reclasificar el modelo de distancias para utilizar solo las mayores a 10 kilómetros, ensamblar los modelos de distribución de la planta, y cruzar el ensamblado con la capa de distancias. 1. Cálculo de las distancias a las poblaciones conocidas de Linaria nigricans. ○

○

Transformación el archivo vectorial Linaria_nigricans_REAL: GIS Manager > File > Map type conversions > Vector to raster, para abrir el módulo v.to.rast. ■

Name of input vector map = Linaria_nigricans_REAL

■

Name for output raster map = Linaria_nigricans_REAL

■

Source of raster values = val

■

Type = area

■

Run

■

Cierra v.to.rast.

Mapa de distancias a las poblaciones: GIS Manager > Raster > Terrain analysis > Cost surface, para abrir el diálogo r.cost. ■

Use the Knight's move = marcado

■

Name of raster map containing grid cell cost information = COSTE (en el mapset PERMANENT)

■

Name for output raster map = DISTANCIA_LINARIA

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 74

Modelos de distribución y estudios de flora amenazada

○

■

Starting points raster map = Linaria_nigricans_REAL

■

Run

■

Cierra r.cost.

Guión de prácticas

Mapa de distancias mayores de 10 kilómetros a cualquier población conocida de Linaria nigricans: Las unidades del mapa de distancias están en decenas de metros. Por tanto, cualquier celda con un valor superior a 1000 estará al menos a 10 kilómetros de una población conocida de la planta. Solo tenemos que reclasificar el mapa DISTANCIA_LINARIA para que tome valor 1 en todas las áreas con valor mayor a 1000, y 0 en todas las áreas con valor menor: GIS Manager > Raster > Change category values and labels > Recode interactively, para abrir la ventana Interactive rules entry. ■

Map to recode = DISTANCIA_LINARIA

■

Recoded map = DISTANCIA_LINARIA_10KM

■

En la ventana escribir = 1000:10000:1 0:999:0

■

Apply

■

Cierra Interactive rules entry.

(El número 10000 es la máxima distancia, redondeada, que aparece en el mapa; el valor aparece durante la creación del mapa, en la ventana Output, como Peak cost value) 2. Ensamblado de modelos. En primer lugar hay que transformar en binarios todos los modelos, utilizando como umbral de corte para cada modelo el que aparece en el apartado PUNTO_DE_CORTE_MEDIA_menos_DESVIACION en el archivo RESULTADOS_AUC.txt (en CURSO_MODELOS/EVALUACION/ AUC). nn representa el valor de umbral de corte para cada modelo:

○

■

CSM

n1

■

ES

n2

■

EDEUC

n3

■

EDCHE

n4

■

EDMAH

n5

■

EDMAN

n6

■

GARPDK

n7

■

GARPOM

n8

■

SVM

n9

Recodificación de los modelos: GIS Manager > Raster > Change category values and labels > Recode interactively, para abrir la ventana Interactive rules entry.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 75

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

■

Map to recode = CSM

■

Recoded map = CSM_M-SD (M-SD, de “mean minus standar deviation”)

■

En la ventana escribir = n1:100:1 0:n1-1:0

■

Apply

■

Aceptar

■

Map to recode = ES

■

Recoded map = ES_M-SD

■

En la ventana escribir = n2:100:1 0:n2-1:0

■

○

Apply, y así sucesivamente hasta haber recodificado los 7 modelos restantes.

Ensamblado de modelos. Los modelos se agregarán mediante suma: Raster > Overlay maps > Map series para abrir el módulo r.series: ■

Name of input raster map(s) = CSM_M-SD,ES_M-SD,EDEUC_MSD,EDMAH_M-SD,EDCHE_M-SD,EDMAN_M-SD,GARPDK_MSD,GARPOM_M-SD,SVM_M-SD

■

Name of output raster map = ENSAMBLADO_M-SD

■

Aggregate operation = sum

■

Run

■

Cierra r.series.

3. Eliminación de áreas que quedan fuera del criterio de mayoría absoluta (aquellas de ENSAMBLADO_M-SD con valores entre 0 y 4): Raster > Develop map > Null values, para abrir r.null. ○

Raster map for wich edit null file = ENSAMBLADO_M-SD

○

List of cell values to be set to null = 0,1,2,3,4

○

Run

4. Aplicación del criterio de mínima distancia: multiplicamos el modelo ENSAMBLADO_M-SD por la capa que contiene el criterio de distancia DISTANCIA_LINARIA_10KM: En la línea de comandos teclea y ejecuta la siguiente expresión (no olvides las comillas en el mapa ENSAMBLADO, o la calculadora confundirá el guión de M-SD con un signo de sustracción): r.mapcalc NUEVAS_POBLACIONES_Ln=”ENSAMBLADO_M-SD”* DISTANCIA_LINARIA_10KM ○

El mapa resultante tiene valores 0 que no nos interesan. Para quitarlos:

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 76

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Raster > Develop map > Null values, para abrir r.null. ■

Raster map for wich edit null file = NUEVAS_POBLACIONES_Ln

■

List of cell values to be set to null = 0

■

Run

■

Cierra r.null.

5. Preparación para visualización. Le pondremos a la capa NUEVAS_POBLACIONES_Ln colores que nos faciliten identificar las áreas con mayor consenso (mayor número de modelos predicen presencia): Raster > Manage map colors > Color tables, para abrir el módulo r.colors. ○

Pestaña Required: ■

○

Pestaña Colors: ■

○

Name of input raster map = NUEVAS_POBLACIONES_Ln Type of color table = gyr

Run

El mapa resultante puede superponerse con transparencia sobre una imagen satélite, o mejor, sobre imágenes aéreas para preparar mapas guía sobre los que diseñar prospecciones.

Figura 22: Ejemplo de un mapa diseñado para la búsqueda de nuevas poblaciones de la planta amenazada Linaria nigricans.

7.2.2 Discusión Este tipo de aplicaciones han sido probadas con éxito por varios investigadores, entre los que nos incluimos. Lo habitual es disponer de un número limitado de Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 77

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

presencias al comienzo, y según se van desarrollando las labores de prospección, las nuevas localidades de presencia que se van obteniendo, deberían utilizarse para recalibrar los modelos y hacer nuevos mapas. Este tipo de realimentación es muy costosa si utilizamos la técnica de ensamblado de modelos, por lo que en este caso, sería más operativo trabajar con un único modelo como guía (el mejor obtenido en una prueba previa, por ejemplo). El ensamblado de modelos, que se basa en buscar coincidencias de criterio entre distintas técnicas, se propone últimamente como una técnica predictiva robusta, y programas de modelado como OpenModeller, que permiten correr varios modelos de una sola vez, facilitan esta labor. Sin embargo, es necesario ser cauto al aplicar esta técnica; como sabemos, varios de los modelos de Linaria nigricans obtuvieron valores de AUC relativamente bajos comparados con el resto (CSM, ES, EDCHE, GARPDK y GARPOM). Y de los cuatro mejores modelos, tres de ellos están basados en un mismo algoritmo (Environmental Distance), por lo que un ensamblado puede ocultar información sesgada.

7.3 ENSAMBLADO DE BIODIVERSIDAD PARA ASISTIR EN EL DISEÑO DE RESERVAS DE FLORA (Tiempo estimado: 40 minutos) El objetivo de esta práctica es obtener mapas de biodiversidad como material de partida para el diseño de redes de reservas de flora. La tarea real debería consistir en modelar la distribución de todas las especies posibles de flora de Almería, utilizando varios modelos de distribución diferentes, y seleccionando la técnica más precisa para cada especie. Con el umbral de corte de cada modelo, se transformarían todos en binarios, y se sumarían, como en la práctica anterior, para obtener el número potencial de especies de flora en cada celda del área de trabajo. Después se establecería un umbral de corte en un número de especies mínimo, y se delimitarían las reservas según características como el tamaño de parche, cercanía a centros urbanos, antropización del medio, etc. Después se comprobaría la información en campo para finalmente abordar el diseño final de la red. Esta tarea es de grandes dimensiones, por lo que vamos a trabajar sobre una simplificación que ilustra los conceptos más importantes. Eliminaremos el modelado con múltiples algoritmos, utilizando solo Environmental Distance Manhattan, porque es rápido en su ejecución, y ya sabemos que proporciona buenos resultados. Trabajaremos con las variables de 1000 metros, para utilizar datos de presencia de GBIF, porque no disponemos de suficientes puntos de presencia de alta resolución. Omitiremos la fase de evaluación de modelos, porque la evaluación individual de un número tan alto de modelos exige mucho tiempo. Finalmente aplicaremos el mismo umbral de corte a todos los modelos, estableciéndolo en 50. ¡Observa que estas prácticas son poco o nada recomendables en un trabajo de producción real!. En este ejercicio aprovecharemos las características de GRASS para realizar el trabajo automáticamente. Para ejecutar la mayoría de las acciones automáticamente crearemos un guión. Es difícil construir un guión tan grande sin que surja algún fallo, así que presta atención a cada paso.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 78

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

7.3.1 Preparación de los modelos de distribución 1. Prepara un nuevo experimento en OpenModeller: ○

Nombre = BIODIVERSIDAD_POTENCIAL.

○

Occurrence data = selecciona ESPECIES_ALMERIA_GBIF.txt, que está en C:\CURSO_MODELOS\ GEODATOS\PRESENCIAS. Contiene los datos de presencia de 225 especies.

○

Algorithm = Environmental Distance Manhattan

○

Model Settings:

○

■

Creation layers = ALMERIA_1000m

■

Projection layers = ALMERIA_1000m

Una vez terminados los modelos, ve a la carpeta C:\CURSO_MODELOS\ RESULTADOS\modelOutputs\BIODIVERSIDAD_POTENCIAL\Environmental_Distance_-_Manhattan, y copia (copiar y pegar, no arrastrar) todos los archivos a la carpeta (tienes que crearla) C:\MODELOS. Esto es para acortar la ruta de acceso en los guiones y facilitar su lectura; ¡OpenModeller crea rutas de acceso muy largas!.

2. Los modelos tardarán en ejecutarse unos minutos, por lo que mientras, en GRASS, prepara un nuevo Mapset llamado BIODIVERSIDAD en la localidad ALMERIA_latlong. Entra en el nuevo mapset, y vamos a prepararlo: ○

Config > Region > Change region settings: ■

Pestaña Existing: ●

○

Set region to match this raster map = 1000m_ITP@PERMANENT

■

Run

■

Cierra r.region.

Raster > MASK, para abrir r.mask: ■

Raster map to use as MASK = 1000m_ITP@PERMANENT

■

Run

■

Cierra r.mask.

3. Preparación del guión para GRASS: 225 modelos son muchos para importarlos a mano, por lo que vamos a crear un guión para importar las imágenes, procesarlas y preparar los resultados: ○

En el menú de Windows, pulsa en Ejecutar, teclea cmd y pulsa sobre Aceptar para abrir la línea de comandos de Windows XP.

○

Teclea y ejecuta la siguiente expresión, para mover la terminal a la carpeta en la que están almacenados los resultados de los modelos: cd C:\MODELOS

Para obtener una lista de los archivos .tiff que necesitamos importar a GRASS y guardarla en CURSO_MODELOS\GUIONES, teclea y ejecuta: dir *.tif >> C:\CURSO_MODELOS\GUIONES\lista_modelos.txt ○

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 79

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

○

Cierra la terminal de Windows.

○

Ve a la carpeta CURSO_MODELOS\GUIONES, y abre con Calc el fichero lista_modelos.txt teniendo en cuenta las siguientes opciones:

○

■

A partir de línea 6

■

Ancho fijo: separa los nombres de los archivos .tif del resto de elementos que van en el fichero.

■

Aceptar.

Ahora haz las operaciones necesarias para obtener un resultado como el que se muestra en la Figura 23. (NOTA: Elimina las dos últimas líneas de la hoja de cálculo, tienen datos que no interesan. Observa que en la columna C, output= lleva un espacio al principio, !no lo olvides!)

Figura 23: Aspecto de la hoja de cálculo para generar la primera parte del guión de GRASS "BIODIVERSIDAD_POTENCIAL.txt" ○

Guarda la hoja como formato .ods con el nombre GUION_BIODIVERSIDAD_POTENCIAL en C:\CURSO_MODELOS\GUIONES

○

Vuelve a guardarla en la misma ubicación como formato .csv, con el nombre BIODIVERSIDAD_POTENCIAL.txt. ■

En Separador de campos, borra lo que venga.

■

En Separador de texto, también borra lo que aparezca.

■

Aceptar.

Si todo ha ido bien, en la carpeta GUIONES ya tienes la primera parte del guión, que importará a la base de datos de GRASS todos los modelos. Si abres el guión con Notepad++, cada línea debería haber quedado así, con espacios delante de -o, input, output r.in.gdal ­o input=C:/MODELOS/Adiantum_capillus_projection.tif  output=Adiantum_capillus ○

La siguiente operación es transformar los modelos en binarios, aplicándoles un umbral de corte con valor 50. Añadiremos las operaciones necesarias al guión anterior, para realizarlas todas en una sola ejecución. ■

Crea en la carpeta CURSO_MODELOS\GUIONES un archivo de texto con el siguiente contenido, y guárdalo como UMBRAL_CORTE.txt: 50:100:1:1 0:49:0:0

■

En la hoja de cálculo GUION_BIODIVERSIDAD_POTENCIAL.ods in-

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 80

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

serta una hoja nueva llamándola CORTE para crear la siguiente parte del guión. Debe quedar como indica la Figura 24 (ten en cuenta los espacios delante de output y rules):

Figura 24: Aspecto de la hoja de cálculo para crear la parte del guión que transforma los modelos continuos en binarios.

○

■

Exporta en formato .csv con el nombre CORTE.txt.

■

Con Notepad++ copia el contenido de CORTE.txt debajo del contenido de BIODIVERSIDAD_POTENCIAL.txt; ahora nuestro guión importará las imágenes y les aplicará el umbral de corte del archivo UMBRAL_CORTE.txt.

■

Guarda los cambios en BIODIVERSIDAD_POTENCIAL.txt, pero no cierres el archivo, porque vamos a añadir dos operaciones más.

Ahora prepararemos la suma de los modelos binarios. Crea una nueva hoja de cálculo llamada CONCATENAR, y pega en ella las columnas de la hoja CORTE con todos los nombres de los modelos, y la que tiene el sufijo _bin (de “binario”). Uniremos el nombre del modelo al sufijo con la función CONCATENAR, tal y como aparece en la Figura 25. Arrastra la función hasta haber unido todos los nombres de modelo con su sufijo.

Figura 25: Función CONCATENAR, para unir nombres de modelos y sufijos. ■

Crea una nueva hoja de cálculo llamada SUMA, en la que pegarás la nueva columna que has creado con la función CONCATENAR, pero utilizando la opción de Pegado Especial, con las siguientes opciones: ●

Selección: ○

Pegar todo = desmarcado

○

Texto = marcado

○

“Todo lo demás” = desmarcado

●

Transponer = marcado

●

Aceptar

■

Si todo ha ido bien, tienes los nombres de los modelos binarios en una sola fila. Guarda la hoja de cálculo como GUION_BIODIVERSIDAD_POTENCIAL.ods. de nuevo, porque tenemos que cambiar los separadores de texto en la exportación.

■

Ahora exporta la hoja SUMA como .csv, con el nombre SUMA.txt.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 81

Modelos de distribución y estudios de flora amenazada ●

Separador de campo = , (coma)

●

Separador de texto = ninguno

●

Aceptar

Guión de prácticas

■

Vuelve a guardar la hoja de cálculo como GUION_BIODIVERSIDAD_POTENCIAL.ods.

■

Abre el archivo SUMA.txt con Notepad++, y al principio y al final de la expresión de suma escribe lo que aparece en negrita en la siguiente expresión. El resultado BIODIVERSIDAD_POTENCIAL_ESPECIES es la suma de todos los modelos binarios:

r.series input=Adiant...mays_bin  output=BIODIVERSIDAD_POTENCIAL_ESPECIES method=sum

○

■

Ahora copia la línea completa, y pégala en BIODIVERSIDAD_POTENCIAL.txt.

■

Vuelve al archivo SUMA.txt, y elimina todos los sufijos “_bin”; de este modo tendremos la misma expresión anterior, pero en este caso para sumar los modelos continuos, y obtener un modelo de biodiversidad basado en la idoneidad. En la expresión de suma, cambia el resultado de BIODIVERSIDAD_POTENCIAL_ESPECIES a BIODIVERSIDAD_POTENCIAL_IDONEIDAD.

■

Ahora copia la línea completa, y pégala en BIODIVERSIDAD_POTENCIAL.txt, debajo de la anterior expresión de suma.

Ejecución del guión: es hora de probar si lo hemos hecho bien. Inicia GRASS Command Line. Comprueba que está en el mapset BIODIVERSIDAD, y pulsa ESC y ENTER consecutivamente para entrar. ■

Para ejecutar el guión, mueve la terminal hasta la carpeta de guiones y ordena la ejecución: cd C:\CURSO_MODELOS\GUIONES sh BIODIVERSIDAD_POTENCIAL.txt

Ahora que se está ejecutando el guión, piensa en el tiempo que hubiera llevado ejecutar manualmente estas acciones. Las técnicas para escribir guiones rápidamente al principio son complejas para quién no tiene excesivo manejo con el PC, pero a la larga permiten ahorrar mucho tiempo. Es posible que la terminal informe de algún error en el guión. Es importante descubrir de que se trata, para corregirlo y volver a ejecutarlo. Los resultados de este guión son BIODIVERSIDAD_POTENCIAL_ESPECIES, que indica el número potencial de especies por kilómetro cuadrado, y BIODIVERSIDAD_POTENCIAL_IDONEIDAD, que es la sumatoria de idoneidad de todos los modelos. Ambos son estimadores de biodiversidad basados en modelos de distribución, pero presentan poca correlación entre sí (0.50), por lo que representan criterios muy diferentes. Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 82

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

Una comparación de estos modelos con la biodiversidad aparente puede ayudarnos a definir criterios. El mapa de biodiversidad aparente indica el número real de especies encontradas en cada celda del area de trabajo, obtenido a partir de los registros de presencia utilizados para hacer los modelos. Para hacer este mapa de biodiversidad aparente, es necesario importar los archivos vectoriales de las especies individuales (que son creados automáticamente por OpenModeller), transformarlos en raster con valor 1 en las celdas ocupadas por registros de presencia, y sumarlos, tal y como hicimos en el guión anterior. Como ya has aprendido a crear un guión para tratamiento masivo de mapas, no es necesario que prepares este, porque ya está en la carpeta GUIONES; se llama BIODIVERSIDAD_APARENTE.txt, y solo tienes que ejecutarlo. Teclea y ejecuta: cd C:\CURSO_MODELOS\GUIONES sh BIODIVERSIDAD_APARENTE.txt Una vez concluya la ejecución, cierra la línea de comandos, y abre GRASS con su interfaz gráfica para visualizar los resultados de ambos guiones (Figura 26).

Figura 26: Comparación de mapas de biodiversidad: a) biodiversidad aparente, o número real de especies por km2; b) biodiversidad potencial, en número potencial de especies por km2: b) sumatorio de idoneidad, para todas las especies modeladas.

Ahora que tienes los resultados, te propongo como ejercicio que selecciones criterios de protección y diseñes una red de reservas lo más razonada posible. 7.3.2 Discusión Tratar de generar una red de reservas a partir de los datos reales de presencia de las especies da una imagen sesgada de la biodiversidad de Almería, y reduce enormemente el número de enclaves con posibilidades de ser incluidos en una red de reservas. Por ejemplo, tomando como criterio un número mínimo de 20 especies para que una celda entre a formar parte de la red, solo 11 km 2 serían protegidos, y muy probablemente, gran cantidad de localidades con alta diversidad que no hayan sido muestreadas quedarían fuera de la red. Además las reservas seleccionadas serían puntuales, sin otra estructura espacial que la definida por un radio concreto alrededor de la celda de presencia. Utilizar la biodiversidad potencial como criterio, según el mismo umbral de 20 especies, supondría la protección de 316 km2, y permitiría trazar reservas con una estructura espacial determinada, más acordes con la realidad del paisaje. Sin embargo el criterio de número de especies excluye localidades con pocas especies que pueden ser importantes por su grado de amenaza o interés biogeográfico Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 83

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

(endémicas, o de distribución restringida, como ejemplos). El grado de amenaza y exclusividad de especies concretas debería ser un criterio más, y podría implementarse tratando por separado estas especies para definir una red de reserva prioritaria, que después se uniría al segundo grupo de especies más genéricas, para formar una red completa. Aún así, la calidad del modelo está supeditada a la calidad del muestreo. el ejemplo claro lo supone la ausencia de áreas a proteger en las montañas (Gádor, Sierra Nevada, etc.), porque pocas especies del conjunto de datos de GBIF han sido recogidas en ellas. Este mismo trabajo, ejecutado a alta resolución, y llevando a cabo todas las tareas de gestión de calidad de los resultados, proporcionaría resultados espacialmente muy precisos, útiles para delinear una red de reservas que ofrecería un grado de protección razonable para áreas no muestreadas. Sin embargo es necesario insistir en que, incluso utilizando la técnica más depurada, un sesgo importante en los datos, y una elección errónea de criterios puede llevar a resultados erróneos.

7.4 EVALUACIÓN DEL IMPACTO POTENCIAL DEL CAMBIO CLIMÁTICO EN LA DISTRIBUCIÓN DE UNA ESPECIE AMENAZADA (Tiempo estimado: 30 minutos) El objetivo de esta práctica es conseguir una toma de contacto con la técnica de proyección de modelos de distribución sobre escenarios de cambio climático, porque una profundización en esta temática exige un curso por sí solo. Prepararemos un modelo de alta resolución (90m) de distribución potencial actual y otro de distribución potencial futura de Linaria nigricans, para compararlos entre sí, y evaluar el impacto potencial del cambio del clima en las poblaciones de la planta. En los últimos años es frecuente encontrar trabajos que predicen la distribución futura de especies animales y vegetales ante distintos escenarios de cambio climático, basados en modelos de distribución de especies y simulaciones del comportamiento del clima para lo que resta de siglo 21. Antes de comenzar con esta práctica debes entender que es un campo teórico en el que la discusión está a la orden del día. Hay una enorme cascada de incertidumbres que afecta a todo el proceso de modelización: ●

Los escenarios de emisiones (Special Report on Emission Scenarios) del IPCC: no se conoce cuál de las líneas evolutivas (A1, A2, B1, B2) es las más probable.

●

Los modelos de predicción climática (Modelos de circulación global océanoatmósfera): se basan en las emisiones de los distintos escenarios SRES, tienen resoluciones muy groseras (celdas de cientos de kilómetros) y es posible que no tengan en cuenta todas las variables que afectan al clima a largo plazo.

●

Los métodos de regionalización de los modelos de predicción climática (transformación de baja resolución a alta resolución espacial): están en fase de desarrollo, y hasta cierto punto se desconoce su capacidad para ajustarse a predicciones a largo plazo.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 84

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

●

Los métodos de modelado de distribución: los distintos algoritmos, con distintas lógicas internas, asumen comportamientos distintos de la misma especie ante cambios en las condiciones climáticas.

●

El nicho ecológico de las especies: no se conoce con exactitud hasta que punto el nicho de una especie permanece estable ante cambios ambientales, y los modelos de distribución aún no asumen esta posibilidad (consideran un nicho ecológico estático). Si una especie es capaz de adaptarse, las predicciones de un modelo de distribución serán erróneas.

Esto es solo un sumario rápido de esta cascada de incertidumbres, y cualquier trabajo como el que desarrollaremos a continuación exige un estudio minucioso de cada una de ellas, sobre todo de las que entran en nuestro dominio, que son las relativas a los aspectos biológicos, relacionados con la ecología de la especie estudiada. 7.4.1 Escenarios de Cambio Climático Para generar mapas de temperatura y precipitación de distintos escenarios de cambio climático se han utilizado los datos regionalizados proporcionados por la AEMET (http://www.aemet.es/es/elclima/cambio_climat/escenarios). Particularmente disponemos para esta práctica de los siguientes mapas: precipitación acumulada anual, temperatura media máxima y temperatura media mínima, correspondientes al escenario CGCM-A2, periodo 2055-2070. Cuando se trabaja con modelos de cambio climático, las variables obtenidas por teledetección (Landsat en este caso) no tienen sentido, porque tienen una evolución temporal que no podemos predecir. Por tanto no utilizaremos para los modelos las variables cpLSAT_1, cpLSAT_2 y NDVI. Las variables para esta práctica están en C:\GEODATOS\VARIABLES, en dos carpetas: PRESENTE_90m y FUTURO_90m. Introdúcelas como conjuntos de variables en OpenModeller, con los mismos nombres. 7.4.2 Proyección de modelos de distribución con OpenModeller Prepara un nuevo experimento de OpenModeller, para conseguir el modelo de distribución de Linaria nigricans correspondiente a l: ●

Nombre = L_nigricans_PRESENTE

●

Occurrence Data = pulsa el signo “+” verde, y añade el archivo Linaria_nigricans.txt (se borrarán los registros de GBIF).

●

Algoritmo = Environmental Distance - Manhattan

●

Creation layers = PRESENTE_90m

●

Projection layers = PRESENTE_90M

●

Ok, para ejecutar el experimento.

Cuando el modelo se ejecute, prepara un nuevo proyecto, para obtener los modelos con la distribución potencial futura de la especie. ●

Nombre = L_nigricans_FUTURO

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 85

Modelos de distribución y estudios de flora amenazada ●

Occurrence Data = Linaria nigricans

●

Algoritmos= Environmental Distance - Manhattan

●

Creation layers = PRESENTE_90m

●

Projection layers = FUTURO_90M

●

Ok, para ejecutar el experimento.

Guión de prácticas

7.4.3 Análisis de los modelos 1. En la localidad ALMERIA_latlong prepara un nuevo mapset llamado PROYECCION_LINARIA. Adapta la región de trabajo a los mapas de 90m: g.region rast=MASK 2. Importación de los modelos: usa el módulo r.in.gdal para importar los modelos. Recuerda que están en CURSO_MODELOS/RESULTADOS/modelOutputs/. Nómbralos Ln_P y Ln_F respectivamente. 7.4.3.1 Distribución potencial actual y futura La idea es conseguir un mapa de distribución potencial actual y distribución potencial futura para calcular el área ocupada en cada momento del tiempo. Es necesario establecer el umbral de corte del modelo Ln_P para transformarlo en binario, y aplicar el mismo criterio al modelo Ln_F. Posteriormente se miden las áreas y se establecen las comparaciones deseadas.

Para poder conocer el umbral de corte, necesitas acceso al archivo de puntos de presencia de Linaria nigricans, que está en el mapset MODELOS. En este momento no tienes acceso al mismo; para abrir este acceso: Config > GRASS Working environment > Mapset acces; marca MODELOS para poder acceder a él, y pulsa OK. Copia el archivo vectorial PRESENCIAS al mapset PROYECCION_LINARIA ejecutando: g.copy vect=PRESENCIAS@MODELOS,PRESENCIAS El resto del procedimiento ya lo has realizado anteriormente: añade una nueva columna al dbf PRESENCIAS del mapset PROYECCION_LINARIA, y muestrea el modelo Ln_P con los puntos de presencia usando v.what.rast. Después abre el .dbf PRESENCIAS, y calcula el promedio que han obtenido sobre el modelo los puntos de evaluación. Utiliza este valor para transformar los modelos Ln_P y Ln_F en binarios. Después calcula el área de ambos, y piensa en como puedes conocer las áreas de persistencia de Linaria nigricans utilizando los mapas resultantes. 7.4.3.2 Diferencial de idoneidad: frentes de avance y retroceso Ante un cambio del clima, se espera un desplazamiento de las condiciones óptimas para la especie. Considerando una población extensa de una especie cualquiera, determinados individuos, por su localización geográfica más marginal, van a experimentar un incremento de estrés fisiológico por el cambio de condiciones. En estas localizaciones se apreciarán antes los efectos del cambio en las poblaciones de la especie; es lo que se denomina frente de retroceso. En cambio, otros individuos experimentarán una mejora en las condiciones ambientales, y estarán en una situaBlas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 86

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

ción privilegiada para que sus descendientes colonicen localidades en las que aparecen condiciones apropiadas; es lo que se denomina frente de avance. A partir de modelos de distribución es sencillo determinar frentes de avance y retroceso, siempre que la variación de condiciones entre los escenarios presente y futuro esté dentro de un margen razonable (por ejemplo, cuando hay un lapso de tiempo relativamente pequeño entre ambos momentos). La única operación a realizar es restar al modelo Ln_F los valores de idoneidad del modelo Ln_P. r.mapcalc DIFERENCIAL=Ln_F-Ln_P Colorea el mapa resultante con r.colors, seleccionando la tabla differences; el mapa resultante es fácil de interpretar: Las áreas de presencia de Linaria nigricans con valor positivo (rojo), en las que se va a incrementar la idoneidad para la especie, se consideran frentes de avance, y las zonas con valores negativos (azul), en las que van a degradarse las condiciones, se consideran frentes de retroceso. Si superpones la cartografía Linaria_nigricans_REAL (mapset MODELOS) que hay en la carpeta CURSO_MODELOS/GEODATOS/PRESENCIAS/Linaria_nigricans, podrás hacer una evaluación visual de cuáles son las poblaciones que mejores posibilidades de persistencia tienen si solo tenemos en cuenta los riesgos derivados del cambio climático. 7.4.4 Discusión Los resultados de estas proyecciones están sometidos a muchas incertidumbres, pero hasta el momento son las herramientas más objetivas de las que disponemos para evaluar los cambios que pueden darse en la distribución vegetal como consecuencia de los cambios que se están produciendo en el clima. Una única proyección sobre un solo escenario climático sirve para establecer el comportamiento potencial de la especie ante uno de los muchos escenarios posibles que pueden darse en el futuro, pero no aporta información sobre la incertidumbre asociada a la predicción. Un método para disminuir este margen de incertidumbre es generar modelos con distintos algoritmos sobre distintos escenarios climáticos.

7.5 ENSAMBLADO DE MODELOS PARA PROYECCIONES DE DISTRIBUCIÓN EN ESCENARIOS DE CAMBIO CLIMÁTICO (Tiempo estimado: 30 minutos) El objetivo de esta práctica es generar modelos de distribución actual y futura de una especie (cualquiera del listado de GBIF) sobre variables de 1000 metros de resolución, utilizando distintos escenarios de cambio climático y ensamblado de modelos. Para abreviar eliminaremos dos pasos fundamentales: partición de datos y evaluación de los modelos. No me cansaré de avisar: en un trabajo de producción real, sería una muy mala idea omitir ambos. 7.5.1 Escenarios climáticos Se utilizan los mismos escenarios presente y futuro del ejercicio anterior, pero en su versión de baja resolución (1000m). Las capas están en las carpetas Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 87

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

PRESENTE_1000m y FUTURO_1000m. Añádelas en OpenModeller como conjuntos de capas con los mismos nombres. 7.5.2 Proyección de modelos de distribución Prepara un nuevo experimento de OpenModeller, para conseguir los modelos de distribución correspondientes a las condiciones actuales: ●

Nombre = ENSAMBLADO_PRESENTE

●

Occurrence Data = selecciona solo una especie del listado de GBIF: tienes que añadir el archivo de texto con los datos.

●

Algoritmos = marca los siguientes: ○

Climate Space Model

○

Environmental Distance

○

Environmental Distance Chebyshev

○

Environmental Distance Mahalanobis

○

Environmental Distance Manhattan

○

GARP with best subsets OpenModeller

○

SVM

●

Creation layers = PRESENTE_1000m

●

Projection layers = PRESENTE_1000m

●

Ejecuta el experimento.

Ahora debes preparar un nuevo proyecto, para obtener los modelos con la distribución potencial futura de la especie. ●

Nombre = ENSAMBLADO_FUTURO

●

Occurrence Data = selecciona la misma especie que en el caso anterior

●

Algoritmos = marca los mismos que en el caso anterior

●

Creation layers = PRESENTE_1000m

●

Projection layers = FUTURO_1000_

7.5.3 Análisis de los modelos entre escenarios 1. Preparación de un nuevo mapset de GRASS: En la localidad ALMERIA_latlong prepara un nuevo mapset llamado PROYECCION_ENSAMBLADO. Adapta la región de trabajo uno de los mapas de 1000m: g.region rast=1000m_ITP@PERMANENT 2. Importación de los modelos: usa el módulo r.in.gdal para importar los modelos. Recuerda que están en CURSO_MODELOS/RESULTADOS/modelOutputs. Nómbralos con las mismas iniciales que en los ejercicios anteriores (CSM, EDEUC, EDCHE, EDMAH, EDMAN, GARPOM, SVM), poniéndoles el sufijo “_P” a los del preBlas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 88

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

sente, y “_F” a los del futuro. 7.4.3.2 Ensamblado de modelos

Para evaluar las diferencias en cuanto a distribución potencial de la planta entre los distintos escenarios, prepararemos un ensamblado con los modelos PRESENTE, y otro con los modelos FUTURO. Se utiliza un ensamblado de modelos porque es un modo de limitar la incertidumbre debida al algoritmo de modelado. 1. Transformación de los modelos en binarios. En este caso no disponemos de presencias de evaluación para determinar los umbrales de corte a partir de los que transformar los modelos en binarios, y con los algoritmos derivados de Environmental Distance, hacerlo a partir de las presencias de calibrado no sirve (todas tienen valor 100). Ante esta circunstancia, y para abreviar, utilizaremos, de nuevo un umbral de corte arbitrario establecido en 50. Transforma en binarios todos los modelos (usando el módulo Recode interactively) con este umbral, poniéndoles el sufijo “_bin” a los nombres originales. 2. Suma los modelos “_F_bin” entre sí, y los modelos “_P_bin” entre sí. 3. Utilizaremos de nuevo un criterio de mayoría absoluta: si al menos cuatro modelos consideran la presencia de la especie, el lugar se define como apto. Reclasifica los resultados de las sumas anteriores, de modo que los lugares aptos tengan valor 1, y los no aptos valor 0. Para reclasificarlos, usa el módulo Interactively rules entry, poniendo en la ventana de entrada de reglas: 0:3:0:0 4:6:1:1 Con estos dos mapas que acabas de crear puedes medir las áreas de ocupación actuales y futuras, para calcular potenciales pérdidas o ganancias, y puedes componer un mapa que muestre visualmente la información, cambiando colores y transparencias de las capas. 7.4.4 Discusión Ya se han comentado las debilidades de estas técnicas, y todas sus aplicaciones están limitadas por las mismas circunstancias. Aún así, es un campo de investigación muy activo, y se está utilizando este método para el diseño de redes de reserva y corredores ecológicos que tengan en cuenta variaciones en la distribución de las especies debidas al cambio climático. En cualquier caso, utilizando datos climáticos regionalizados como los que proporciona la AEMET, variables de alta resolución, un tamaño apropiado de tamaño de muestra, y métodos robustos de evaluación, selección y análisis de modelos, los resultados muestran una significación mucho mayor que la que hayas podido obtener con el rápido ejemplo que has trabajado.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 89

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

7.6 EJERCICIOS PROPUESTOS 7.6.1 Cartografía de poblaciones En la hoja de cálculo PRESENCIAS, en la pestaña GPS tienes registros de GPS de distintas plantas de la provincia de Almería. Genera una cartografía de detalle (90m) De una de ellas, utilizando todos los pasos necesarios para asegurar la calidad del proceso de modelado: partición aleatoria, modelado con múltiples algoritmos, evaluación y selección de modelos y selección de un umbral de corte adecuado. Mide el área de ocupación potencial, y prepara una visualización de la cartografía sobre una composición de color natural obtenida de las bandas Landsat. 7.6.2 Búsqueda de nuevas poblaciones A partir de los resultados de la actividad anterior, prepara una cartografía apropiada para diseñar prospecciones de campo con el objetivo de encontrar nuevas poblaciones de la especie con la que trabajas. 7.6.7 Explora libremente Piensa en lo que has aprendido, y piensa en modos de extraerle jugo. Si tienes dudas, usa los datos para experimentar, y conocer el funcionamiento y utilidad de los modelos. Inventa nuevos métodos de análisis de los modelos y aplicaciones prácticas en tus trabajos de conservación de plantas.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 90

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

BIBLIOGRAFÍA Allouche O, Tsoar A, Kadmon R (2006). Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). Journal of Applied Ecology, 43, 1223-1232. Araújo MB, New M (2006). Ensemble forecasting of species distribution. Trends in Ecology and Evolution, 22, 1 Araújo MB, Guisan A (2006). Five (or so) challenges for species distribution modelling. Journal of Biogeography, 33, 1677-1688. Anderson RP, Lew D, Peterson AT (2003). Evaluating predictive models of species' distributions: criteria for selecting optimal models. Ecological Modelling, 162, 211232. Barry S, Elith J (2006). Error and uncertainty in habitat models. Journal of Applied Ecology, 43, 413-423. Cohen J (1960). A coefficient of agreement for nominal scales. Educational and Phychological Measurement, 41, 687-699. Elith J, Burgman MA, Regan HM (2002). Mapping epistemic uncertainties and vague concepts in predictions of species distribution. Ecological Modelling, 157, 313-329. Fielding AH, Bell JF (1997). A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation, 24, 38-49. Guisan A, Thuiller W (2005). Predicting species distribution: offering more than simple habitat models. Ecology Letters, 8, 993-1009. Guisan A, Zimmermann en (2000). Predictive habitat distribution models in ecology. Ecological Modelling, 135, 147-186. Hirzel AH, Le Lay G, Helfer V, Randin C, Guisan A (2006). Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling, 199, 142-152. Jiménez-Valverde A, Lobo JM, Hortal J (2008). Not as good as they sem: the importance of concepts in species distribution modelling. Diversity and Distributions, DOI: 10.1111/j.1472-4642.2008.00496.x Jiménez-Valverde A, Lobo JM (2007). Threshold criteria for conversion of probability of species presence to either-or presence-absence. Acta oecologica, 31, 361-369. Liu C, Berry PM, Dawson TP, Pearson RG (2005). Selecting thresholds of occurrence in the prediction of species distributions. Ecography, 28, 385-393. Lobo JM, Jiménez-Valverde A, Real R (2007). AUC: a misleading measure of the preformance of predictive distribution models. Global Ecology and Biogeography, Online Article, DOI: 10.1111/j.1466-8238.2007.00358.x Manel S, Williams HC, Ormerod SJ (2001). Evaluating presence-absence models in ecology: the need to account for prevalence. Journal of Applied Ecology, 38, 921-931. Phillips SJ, Anderson RP, Schapire RE (2006). Maximum entropy modeling of speBlas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 91

Modelos de distribución y estudios de flora amenazada

Guión de prácticas

cies geographic distributions. Ecological Modelling, 190, 231-259. Raes N, Steege H (2007). A null-model for significance testing of presence-only species distribution models. Ecography, 30, 727-736. Rodríguez JP, Brotons L, Bustamante J, Seoane J (2007). The application of predictive modelling of species distribution to biodiversity conservation. Diversity and Distributions, 13, 243-251. Soberón J (2007). Grinnellian and Eltonian niches and geographic distributions of species. Ecology Letters, 10, 1115-1123. Thuiller W, Albert C, Araújo MB, Berry PM, Cabeza M, Guisan A, Hickler T, Midgley GF, Paterson J, Schurr FM, Sykes MT, Zimmermann EN (2008). Predicting global change impacts on plant species' distributions: Future challenges. Perspectives in Plant Ecology, Evolution and Systematics, 9, 137-152.

Blas Benito de Pando; Unidad de Conservación Vegetal, Universidad de Granada

página 92