UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO POSGRADO EN CIENCIA E INGENIERÍA DE LA COMPUTACIÓN “OBTENCIÓN EN PARALELO DE ISOYETAS CON OROGRAFÍA” T E S I S QUE PARA OBTENER EL GRADO DE MAESTRA EN CIENCIAS (COMPUTACIÓN) P R E S E N T A: ROSA ETNA CERVANTES CAMACHO DIRECTOR DE TESIS: DR. EDGAR GARDUÑO ÁNGELES México, D. F. 2011. UNAM – Dirección General de Bibliotecas Tesis Digitales Restricciones de uso DERECHOS RESERVADOS © PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL Todo el material contenido en esta tesis esta protegido por la Ley Federal del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México). El uso de imágenes, fragmentos de videos, y demás material que sea objeto de protección de los derechos de autor, será exclusivamente para fines educativos e informativos y deberá citar la fuente donde la obtuvo mencionando el autor o autores. Cualquier uso distinto como el lucro, reproducción, edición o modificación, será perseguido y sancionado por el respectivo titular de los Derechos de Autor. i Este trabajo está dedicado a Eduardo, Rosa, Carlos, Augusto, Carla e Ilse con todo mi agradecimiento por literalmente sostenerme la vida y el alma. Ustedes me ayudaron a convertir un lejano anhelo en realidad. iii Resumen Este trabajo se realizó como tesis de maestría en Ciencia e Ingeniería de la Computación en la UNAM. Los métodos de interpolación de precipitación que se proponen en esta tesis son Vecinos Voronoi y Vecinos Voronoi con Orografía, 4 o 16 estaciones cercanas por el inverso de la distancia y por el inverso de la distancia al cuadrado, 4 o 16 estaciones cercanas con Orografía. Las interpolaciones aquí propuestas estiman el valor de precipitación entre estaciones de medición o pluviómetros sobre el territorio dado y para un periodo específico. No son métodos de predicción de clima. Las mediciones están dispuestas de manera arbitraria, dispersa y no regular sobre el territorio. Para la implementación y experimentación realizada en esta tesis se tomaron los datos de precipitación diaria de Aguas Superficiales de la Comisión Nacional del Agua (CONAGUA), usados bajo un permiso especial. Se dispuso de toda la historia completa de los datos. Para los cálculos se consideraron precipitaciones anuales sobre décadas específicas. La consideración de los datos del conjunto de 1,460 estaciones estuvo sujeta a que las series pasaran un control de calidad que toma en cuenta las recomendaciones en la materia de la Organización Mundial Meteorológica (OMM). La orografía utilizada fue obtenida del modelo digital de elevación con resolución de 30 m, obtenido a partir de radar satelital, publicado por la Agencia Espacial Estadounidense (NASA) en junio 2009 y que es del dominio público. Se diseñaron los métodos de interpolación propuestos, se programaron éstos y los métodos de interpolación más utilizados: Kriging y Spline de curvatura continua sujeto a tensión, que son los métodos utilizados por Matlab©, Surfer© y GMT®, con la finalidad de tener uniformidad para hacer la evaluación de desempeño de los métodos de interpolación propuestos. Se hicieron experimentos de interpolación con datos anuales sobre diferentes periodos y se seleccionaron los de la década de 1971-1980 para hacer las valoraciones. La calidad de las interpolaciones se midió contra el estándar en cálculo de precipitación areal, el método de Thiessen, pero también se hizo una medición cuantitativa con la técnica de verificación cruzada. Los nueve métodos propuestos, Thiessen con orografía, Vecinos Voronoi y con orografía, 4/16 estaciones cercanas ponderando por el inverso de la distancia o por el inverso de la distancia al cuadrado y 4/16 estaciones cercanas con orografía producen interpolaciones de precipitación con las siguientes características: 1. Local. El valor en un punto de la superficie está determinado por unas cuantas estaciones, el número de ellas la determina el método. Vecinos Voronoi toma en cuenta la estación más cercana y las que quedan a la vista. El método de las 4 estaciones cercanas toma en cuenta el valor de las cuatro estaciones con menor distancia euclidiana sobre (X, Y). El de 16 las dieciséis correspondientes. 2. Adaptable. A la disponibilidad de mediciones en la zona. En las áreas donde hay mucha medición, los datos utilizados para la interpolación son más cercanos al punto; en donde no hay medición cercana son utilizados datos lejanos, los disponibles. 3. No lineal. Por la manera de cálculo los métodos producen superficies no lineales. iv 4. Discontinuo. Vecinos Voronoi. La superficie producida por la interpolación a través del método Vecinos Voronoi es discontinua, no diferenciable en todas las aristas del diagrama de Voronoi donde las estaciones contiguas hayan medido valores diferentes. 5. Continuo. 4 /16 estaciones cercanas. Es continua y diferenciable en todos los puntos la superficie producida por la interpolación a través de los métodos 4 / 16 estaciones cercanas. 6. Complejidad lineal en tiempo de ejecución y espacio de almacenamiento. 7. La precisión de la interpolación de los nueve métodos propuestos se describe a continuación: La mínima diferencia porcentual fue prácticamente cero para cualquier método. Ello se debe a que todos honran los puntos de medición. Thiessen con orografía presentó una diferencia máxima porcentual de 600%, del mismo orden que el método de cuatro estaciones cercanas por el inverso de la distancia al cuadrado. Eso es exactamente el mismo orden de magnitud que las diferencias máximas porcentuales en el valor interpolado por los métodos Voronoi, Voronoi con orografía, cuatro estaciones cercanas por el inverso de la distancia, dieciséis estaciones cercanas por el inverso de la distancia y por el inverso de la distancia al cuadrado, cuatro y dieciséis estaciones cercanas con orografía que están alrededor de 400%. A pesar de ello el promedio de las diferencias porcentuales del valor interpolado por los nueve métodos propuestos es muy bajo, del 32% en promedio y la varianza es tan sólo del 29%. Esto demuestra la mayor calidad y precisión de los métodos propuestos en comparación a los índices obtenidos por los métodos utilizados actualmente, tal y como se detalla en el capítulo de resultados. 8. Producen mapas adecuados para la toma de decisiones ya que pueden hacerse sobre ellos lecturas utilizables, debido a la resolución espacial con la que pueden ser generados (una milésima de grado) y debido a que los mapas generados no producen grandes desviaciones atribuibles a aberraciones matemáticas, tales como las que pueden ser producidas con las interpolaciones Kriging y Spline de curvatura continua, que bajo determinados parámetros éstos métodos globales entregan cosas como precipitaciones negativas y valores mucho más altos que la precipitación máxima histórica registrada. 9. Los mapas con orografía tienen una correlación mayor con la vegetación del territorio vista en las imágenes de satélite y comparada de manera cualitativa. El método que nos parece más adecuado para la interpolación de precipitación es Vecinos Voronoi con Orografía ya que además de las características arriba descritas es geométricamente correcto al aportar toda la información cercana disponible y adaptándose a las características espaciales de la medición. Agradezco al Dr. Edgar Garduño por su trabajo, colaboración y guía. Agradezco al Posgrado de Ciencia e Ingeniería de la Computación (PCIC) de la Universidad Nacional Autónoma de México (UNAM) por haberme permitido participar en su programa de posgrado. Agradezco al Departamento de Ciencias de la Computación (DCC) del Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas (IIMAS). Este trabajo tuvo el apoyo de CONACYT a través del convenio número 290541. El detalle y las descripciones completas de lo que está descrito en este resumen está en el contenido de la tesis. v Agradecimientos Agradezco al Consejo Nacional de Ciencia y Tecnología (CONACYT) por el apoyo brindado a través del convenio número 290541. Agradezco al Posgrado en Ciencia e Ingeniería de la Computación (PCIC) de la Universidad Nacional Autónoma de México (UNAM) por haberme permitido participar en su programa. El alto nivel de sus profesores y administrativos lo hacen ser el mejor en México y tener reconocimiento internacional. Agradezco al Departamento de Ciencias de la Computación del Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas (IIMAS) por permitirme utilizar sus instalaciones y por todo el apoyo con el aporte de recursos materiales para la realización de esta tesis. Agradezco a mis sinodales Dr. Luis Miguel de la Cruz, Dr. Jorge Lira, Dr. Jorge Urrutia y Dr. Héctor Benítez por su colaboración en este trabajo. Todos sus apuntamientos me permitieron entender, aprender y mejorar. Muchas gracias. Agradezco a mi estimado Dr. Edgar Garduño Ángeles por su trabajo y buena guía. Aprendí mucho de su criterio agudo y su integridad perenne. Me enseñó mucho con su trabajo disciplinado y constante; las reuniones semanales de una hora a través de los más de dos años de tutela acumulan una gran cantidad de trabajo reflejada en esta tesis. En todas las sesiones hubo productivas discusiones y agudos cuestionamientos que le agradezco profundamente. Le agradezco mucho su paciencia conmigo y sobre todo haber enfrentado el reto de dirigir esta tesis de características tan sui géneris. Agradezco a la Dra. María Elena Martínez y al Dr. Jorge Márquez por el trabajo de revisión, análisis y aportaciones que ayudaron mucho a enfocar y delimitar algunos aspectos de este trabajo. Agradezco a mis compañeros y profesores, amigos todos, por su colaboración. En especial a Jorge García, Adolfo García y Héctor Benítez por las largas conversaciones, sabios y útiles consejos y por la valiosa compañía. Agradezco a mis colaboradores y amigos de la Organización Mundial Meteorológica (OMM), la Comisión Nacional del Agua (CONAGUA) y el Servicio Meteorológico Nacional (SMN) por su apoyo y valoración de los resultados de esta tesis. Agradezco a Augusto, Carla e Ilse por la inspiración y el ejemplo que me dan sus vidas, pero también por el apoyo que fue indispensable para la realización de este trabajo. Agradezco a Rosa y Carlos por compartirme sus sueños y anhelos, su vida y sus días. Son los padres más maravillosos de la historia de los mundos. Su amor me enseñó a amar la vida y me ha sostenido a través de las dificultades. Gracias por compartirme su maravillosa vida. Agradezco a mi amadísimo Eduardo por la entereza con la que me protege, porque invariablemente enfrenta las dificultades que nos ha presentado la vida con tranquilidad y una sonrisa, por todo su amor y finos cuidados, que transforman este mundo hostil y difícil en un paraíso donde da gusto vivir, me permite soñar y colabora activamente para transformar mis deseos en realidad cotidiana. Agradezco a toda la gente que paga impuestos y a la que cumple con su trabajo en paz. vii Causa primera de mi agradecimiento la origina profundamente Dios por haberme regalado vida y salud para realizar este trabajo. Deseo que sea para su servicio. ix Contenido 1. Introducción ................................................................................................................ 1 1.1. Contexto y motivación ........................................................................................... 1 1.2. Antecedentes ....................................................................................................... 2 1.3. De cómo es la lluvia .............................................................................................. 3 1.3.1. Escala de medición. ....................................................................................... 4 1.3.2. Suficiencia de medición. ................................................................................. 4 1.3.3. Área de cobertura de la medición. ................................................................... 4 1.3.4. Días con lluvia. .............................................................................................. 7 1.3.5. Lluvias de un día. .......................................................................................... 7 1.3.6. Distribución espacial de las estaciones. .......................................................... 10 1.3.7. Altitud de las estaciones. .............................................................................. 10 1.3.8. Generalidades. ............................................................................................ 11 1.4. Hipótesis ............................................................................................................ 11 1.5. Objetivos ........................................................................................................... 12 1.6. Metodología ....................................................................................................... 12 1.7. Las partes que conforman a este documento ........................................................ 13 2. Conceptos preliminares ............................................................................................... 15 2.1. Definiciones ....................................................................................................... 15 2.2. Mapas actuales ................................................................................................... 17 2.3. De la orografía y la precipitación .......................................................................... 18 2.4. Interpolación por Kriging y Spline CCST ................................................................ 20 2.4.1. Kriging........................................................................................................ 21 2.4.2. Spline de curvatura continua sujeto a tensión ................................................ 27 3. Algoritmos ................................................................................................................. 33 3.1. Definición del problema ....................................................................................... 33 3.1.1. Definición de precipitación diaria ................................................................... 37 3.1.2. Datos dispersos ........................................................................................... 37 3.1.3. Cálculo de malla .......................................................................................... 38 3.1.4. Puntos de cálculo ........................................................................................ 38 3.1.5. Distancia geodésica ..................................................................................... 38 3.1.6. Determinación de interioridad de un punto a un polígono cerrado.................... 39 3.1.7. Puntos interiores ......................................................................................... 39 3.1.8. Criterio de cercanía ...................................................................................... 40 3.1.9. Cálculo de la distancia .................................................................................. 41 3.1.10. Determinación de la diferencia de altitud ....................................................... 41 3.2. Vecinos Voronoi .................................................................................................. 41 3.3. Vecinos Voronoi con Orografía ............................................................................. 43 3.4. 4 / 16 Estaciones Cercanas por el inverso de la distancia y por el inverso de la distancia al cuadrado. .................................................................................................... 44 3.5. 4 / 16 Estaciones Cercanas con Orografía ............................................................. 45 3.6. Thiessen con Orografía ....................................................................................... 46 3.7. Cálculos areales .................................................................................................. 47 3.8. Trazado de isoyetas ............................................................................................ 48 3.9. Cálculo de la constante de orografía ..................................................................... 49 x 4. Implementación ......................................................................................................... 51 4.1. Datos utilizados .................................................................................................. 51 4.1.1. Precipitación diaria ...................................................................................... 51 4.1.2. Cartografía digital ........................................................................................ 52 4.1.3. Modelo digital de elevación........................................................................... 52 4.2. Cálculo por malla ................................................................................................ 54 4.3. Programas ......................................................................................................... 56 4.3.1. Programa de generación de mallas de preprocesamiento de orografía.............. 56 4.3.2. Rutina de consulta a la base de datos climatológica. ....................................... 56 4.3.3. Rutinas de los 16 métodos de cálculo de la superficie de precipitación con salidas parciales en imagen, archivos de datos ........................................................................ 56 4.3.4. Procesamiento paralelo ................................................................................ 59 4.3.5. Lista de rutinas y otros programas. ............................................................... 59 4.4. Radio máximo de competencia ............................................................................ 60 4.5. Procesamiento en paralelo ................................................................................... 60 5. Resultados ................................................................................................................. 63 5.1. Resultados de los experimentos realizados ............................................................ 63 5.1.1. Experimentos realizados ............................................................................... 63 5.2. De los métodos de validación............................................................................... 84 5.2.1. Comparación contra datos artificiales ............................................................ 84 5.2.2. Comparación por el método de validación cruzada de datos ............................ 85 5.3. Validación comparando contra con datos artificiales ............................................... 86 5.4. Validación usando referencia cruzada ................................................................... 92 5.5. Validación comparando contra con isoyetas manuales históricas ............................. 94 5.6. Validación cualitativa contra imágenes de satélite .................................................. 96 6. Conclusiones .............................................................................................................. 97 Anexo - Siglas, acrónimos, nombres y marcas ..................................................................... 99 Referencias .....................................................................................................................101 Índice .............................................................................................................................103 xi Lista de Figuras Fig. 1 Área total nacional que sería cubierta si a cada estación se le asigna un radio de cobertura determinado. Son pocas estaciones para un vasto territorio y el radio de cobertura es el área donde tiene validez cada dato. ............................................................................................. 5 Fig. 2 Área de cada una de las 37 Regiones Hidrológicas comparada con la parte cubierta por las estaciones si se consideran diferentes radios de cobertura. La cobertura total de los territorios se alcanza con radios mayores o iguales a 200 km. ................................................................ 6 Fig. 3 Precipitación diaria de la estación Tierra Blanca Veracruz de junio a julio de 1971. En este periodo ha caído la mitad (3,377 mm) de la precipitación anual acumulada (7,216 mm). .......... 8 Fig. 4 En rojo la precipitación diaria de la estación de la Presa Taxhimay en el estado de Hidalgo, en azul la acumulada. Datos de junio y julio de 1971. Es una zona donde las precipitaciones los días del año que sí hay lluvia se parecen en magnitud mucho unas a otras, todas debajo de los 50 mm. ........................................................................................................................... 9 Fig. 5 Gráfica de la intensidad de lluvia diaria respecto al porcentaje de la lluvia anual total. Porcentaje de precipitación anual respecto a la anual por estación de los datos de la década 1971 a 1980. En rojo se muestra la lluvia acumulada por frecuencia de intensidad de lluvia y en azul el valor acumulado. ...................................................................................................... 9 Fig. 6 Gráfica de la intensidad de lluvia media diaria y de la máxima lluvia diaria registrada en ese intervalo comparadas contra el porcentaje de lluvia anual. ................................................... 10 Fig. 7 Distribución espacial de las estaciones de medición pluviométrica. Las zonas urbanas y económicamente muy activas son las que tienen mayor densidad de estaciones. ................... 11 Fig. 8 Mapa oficial del Servicio Meteorológico Nacional, CONAGUA (tomado de su página Web), el cual muestra la precipitación media anual. Está indicado que la lámina media anual es de 773.5 mm. Puede observarse en el Golfo de California que el método matemático de interpolación asume continuidad del terreno y en la Península de Yucatán menor precipitación en la costa. ............ 18 Fig. 9 Se muestra el simulador de la termodinámica de la atmósfera que muestra el efecto del enfriamiento adiabático de la atmósfera. El programa pertenece a la Universidad de Iowa, EUA. . ......................................................................................................................... 19 Fig. 10 Variograma el cual describe la rugosidad espacial de un conjunto de datos. ........................... 22 Fig. 11 Superficie de precipitación obtenida mediante el método Kriging Gaussiano. El color indica cuál es la precipitación media anual con los datos en la década 1971 a 1980. El color verde significa que llueve igual en la península de Baja California que en la Cuenca Lerma e igual que en Guerrero. ..................................................................................................................... 23 Fig. 12 Mismos datos que la figura anterior pero vista en perspectiva. Así se observa que el radio de competencia de cada estación es demasiado pequeño. Los parámetros con los que se interpola el método de Kriging son críticos para la calidad del resultado. ............................................. 24 Fig. 13 Interpolación Kriging de tipo exponencial. Este tipo de interpolación por el método de Kriging es el que da resultados con más sentido. Puede verse perfectamente en el mapa las zonas secas y húmedas. ..................................................................................................................... 25 xii Fig. 14 Acercamiento y vista en perspectiva de la Interpolación de Kriging tipo Exponencial de la figura anterior en donde se aprecia la continuidad de la superficie y cómo cada estación parece un pico o depresión en ella. ........................................................................................................... 25 Fig. 15 Superficie de precipitación interpolada por el método de Kriging Esférico. Pueden apreciarse detalles como que en medio de la selva Lacandona no hay medición y en todos los mapas se cae la superficie interpolada en esa zona. ........................................................................... 26 Fig. 16 Recorte de la vista en perspectiva de la superficie de precipitación interpolada por Kriging Esférico. Se aprecia claramente cómo las estaciones van dando forma a la interpolación y cómo aberraciones matemáticas producen zonas de precipitación negativa. ................................... 27 Fig. 17 Spline de Curvatura Continua Sujeto a Tensión con valor T=0, llamado Laplaciano. Se aprecia la lámina uniformemente deprimida en el noroeste y medianamente alta en el sureste. Las estaciones parecieran tener poco peso individual y más como conjunto hablando espacialmente. ..................................................................................................................... 30 Fig. 18 Método de interpolación de precipitación Spline de Curvatura Continua Sujeto a Tensión con valor T=1, llamado Biarmónico. Presenta una zona de precipitaciones negativas entre Oaxaca y Guerrero, en una zona sin medición. .................................................................................. 30 Fig. 19 El mejor desempeño del método de Spline de Curvatura Continua Sujeto a Tensión se da con el parámetro recomendado de T = 0.75. ................................................................................ 31 Fig. 20 Para el punto en la malla representado por el cruce de las líneas rojas en blanco se muestra el polígono de Voronoi al que pertenece y en gris se muestra la primera vecindad o también llamada 1-vecindad de Voronoi. ......................................................................................... 42 Fig. 21 Triangulación de la malla para el trazado de isoyetas. .......................................................... 48 Fig. 22 Un año de las series de datos diarios de las estaciones Cerro Catedral, Nevado de Toluca y San Juan de las Huertas. La constante de orografía se determinó 0.4 por cada 500 m de diferencia de altura. ..................................................................................................................... 50 Fig. 23 Ubicación de las estaciones pluviométricas con registro diario de la CONAGUA. ..................... 52 Fig. 24 Índice de archivos DEM del Modelo Digital de Elevación ASTER Global de la NASA. Mosaicos de 1º × 1º. Escala de color arbitraria. ..................................................................................... 53 Fig. 25 Malla preprocesada a 0.02 del modelo digital de elevación ASTER Global de la NASA visualizado mediante el programa desarrollado. Escala de color arbitraria. .............................................. 53 Fig. 26 Cálculos resultantes en una malla de 20 km después de ser perfilado y la imagen es el estado previo a la rotulación. Son aproximadamente un millón de puntos calculados aunque la malla contiene en realidad millón y medio 1600 × 950. ................................................................. 54 Fig. 27 Pruebas de cobertura con diferentes radios de competencia de las estaciones........................ 60 Fig. 28 Interpolación por Vecinos Voronoi. ..................................................................................... 65 Fig. 29 Interpolación por Vecinos Voronoi con Orografía. ................................................................ 66 Fig. 30 Interpolación por 4 Estaciones Cercanas por el inverso de la distancia. .................................. 67 Fig. 31 Interpolación por 16 Estaciones Cercanas por el inverso de la distancia. ................................ 67 Fig. 32 Interpolación por 4 Estaciones Cercanas por el inverso de la distancia al cuadrado. ................ 68 Fig. 33 Interpolación por 16 Estaciones Cercanas por el inverso de la distancia al cuadrado. .............. 68 Fig. 34 Interpolación por 4 Estaciones Cercanas con Orografía. ....................................................... 69 Fig. 35 Interpolación por 16 Estaciones Cercanas con Orografía. ...................................................... 69 xiii Fig. 36 Interpolación por Thiessen con Orografía. ........................................................................... 70 Fig. 37 Mapa de isoyetas de Vecinos Voronoi con Orografía a cada 500 mm. .................................... 71 Fig. 38 Mapa de interpolación por Vecinos Voronoi con Orografía con isoyetas a cada 500 mm. ......... 71 Fig. 39 Trazo de las isoyetas por el método de marching triangles. .................................................. 72 Fig. 40 Acercamiento que muestra la ubicación de las estaciones en los mapas. ................................ 72 Fig. 41 Interpolación por Vecinos Voronoi con Orografía mostrada en Global Mapper®. ...................... 73 Fig. 42 Comparativo de la precipitación areal nacional calculada con mallas a 1 km, 2 km, 5 km y 10 km. No hay diferencia significativa entre ellas. ..................................................................... 74 Fig. 43 Comparativo de las áreas totales conformadas por los elementos de las celdas de la malla a resoluciones de 1 km, 2 km, 5 km y 10 km comparadas con las áreas de los polígonos originales. ..................................................................................................................... 75 Fig. 44 Comparativo del cálculo de la precipitación areal por Regiones Hidrológicas calculada por mallas a diferentes resoluciones: 1 km, 2 km, 5 km y 10 km. .......................................................... 75 Fig. 45 Comparativo de las áreas de las Regiones Hidrológicas calculadas por cada elemento de la malla a diferentes resoluciones: 1 km, 2 km, 5 km y 10 km y comparada con las áreas de los polígonos originales. ......................................................................................................... 76 Fig. 46 Mapas de las épocas de secas y de lluvias de cuatro décadas diferentes. Las diferencias se aprecian principalmente en las épocas de lluvias.................................................................. 77 Fig. 47 Gráfica de tiempos separada por procesos para cada uno de los métodos de interpolación cuando se procesa en serie a un sólo CPU. .......................................................................... 78 Fig. 48 Gráfica que muestra la aceleración producida al procesar con 4 CPUs comparada con 1 CPU para todos los métodos de cálculo propuestos. ........................................................................... 78 Fig. 49 Gráfica de aceleración del procesamiento del cálculo de los algoritmos que incluye un Kriging. 79 Fig. 50 Gráfica de tiempos de proceso diferenciada por etapas. Se omitieron los Kriging para tener una visión clara de comparación. .............................................................................................. 79 Fig. 51 Similar a la figura anterior se muestra la Gráfica de Tiempos de proceso diferenciado por etapas que incluye los métodos Kriging. ........................................................................................ 80 Fig. 52 Precipitación areal anual nacional con datos en la década 1971 a 1980 calculada por los diferentes métodos. .......................................................................................................... 80 Fig. 53 Precipitación areal anual por Regiones Hidrológicas comparada entre los diferentes métodos de cálculo. ..................................................................................................................... 81 Fig. 54 Precipitación areal anual por Organismos de Cuenca comparada entre los diferentes métodos. 81 Fig. 55 Precipitación areal anual por Estados parte 1 comparada entre los diferentes métodos. .......... 82 Fig. 56 Precipitación areal anual por Estados parte 2 comparada entre los diferentes métodos. .......... 82 Fig. 57 Precipitación areal anual nacional calculada con datos de diferentes décadas y calculada por el método de Vecinos Voronoi. .............................................................................................. 83 Fig. 58 Comparación de la precipitación areal anual tomando datos de diferentes décadas y calculado por el método de Vecinos Voronoi con Orografía. ................................................................ 83 Fig. 59 Precipitación areal nacional anual por años calculada por Vecinos Voronoi con Orografía. ....... 84 Fig. 60 Error del método Vecinos Voronoi contra Thiessen. .............................................................. 87 Fig. 61 Error del método Vecinos Voronoi con Orografía contra Thiessen. ......................................... 87 xiv Fig. 62 Error del método 4 Estaciones Cercanas por el inverso de la distancia comparado contra Thiessen. ..................................................................................................................... 88 Fig. 63 Error del método 16 Estaciones Cercanas por el inverso de la distancia comparado contra Thiessen. ..................................................................................................................... 88 Fig. 64 Error del método 4 Estaciones Cercanas por el inverso de la distancia al cuadrado comparado contra Thiessen. ............................................................................................................... 89 Fig. 65 Error del método 16 Estaciones Cercanas por el inverso de la distancia al cuadrado comparado contra Thiessen. ............................................................................................................... 89 Fig. 66 Error del método 4 Estaciones Cercanas con Orografía comparado contra Thiessen. ............... 90 Fig. 67 Error del método 16 Estaciones Cercanas con Orografía comparado contra Thiessen. ............. 90 Fig. 68 Error del método Thiessen con Orografía comparado contra Thiessen. .................................. 91 Fig. 69 Error del método Kriging Gaussiano comparado contra Thiessen. .......................................... 91 Fig. 70 Error del método Spline de Curvatura Continua Sujeto a Tensión con T=0.75 comparado contra Thiessen. ..................................................................................................................... 92 Fig. 71 Tres experimentos con el 20% de las estaciones omitidas y 3 con el 30% para medir los valores interpolados contra las mediciones reales para el método Vecinos Voronoi con Orografía. ....... 93 Fig. 72 Comparación de los datos de precipitación y área por Regiones Hidrológicas calculada por el método de Vecinos Voronoi contra los datos manuales CONAGUA 1930 – 1990. 100% quiere decir que son iguales. ....................................................................................................... 94 Fig. 73 Interpolación por Vecinos Voronoi con Orografía sobrepuesta a las imágenes de satélite de Google Earth©. ................................................................................................................. 96 Fig. 74 Acercamiento que muestra la enorme correlación entre la intensidad de precipitación y la vegetación subyacente en las imágenes satelitales............................................................... 96 xv Lista de Tablas Tabla 1 Precipitación diaria comparada con la anual de una estación de desierto. Puede observarse que no es raro que en un día llueva el 50% de la precipitación total anual; y más que la precipitación anual acumulada en varios años. .......................................................................................... 7 Tabla 2 Muestra del cálculo de las áreas con una malla de 0.02 y con una de 0.05. La malla de 0.02 da una diferencia de dos centésimas de punto porcentual mientras que la de 0.05 da ocho centésimas. ..................................................................................................................... 55 Tabla 3 Comparación de las áreas del territorio mexicano calculado con nuestras mallas de 0.02 y 0.05 contra los datos oficiales. La diferencia es muy pequeña. ..................................................... 56 Tabla 4 Parámetros de comando de línea. ....................................................................................... 57 Tabla 5 Parámetro de tipo de cálculo para el programa Etisol. .......................................................... 58 Tabla 6 Archivos que se general al correr Etisol. .............................................................................. 58 Tabla 7 Tabla de métodos de interpolación propuestos y contra los que se compara. ......................... 64 Tabla 8 Tabla de evaluación del desempeño de los métodos de interpolación usando referencia cruzada. ..................................................................................................................... 92 Tabla 9 Tabla de datos donde se comparan la precipitación areal por Región Hidrológica y las áreas de las mismas del método Vecinos Voronoi contra los datos manuales CONAGUA 1930 – 1990. 100% quiere decir que son iguales. .................................................................................... 95 xvii 1 1 1. Introducción En este trabajo se presentan nueve maneras de hacer los mapas de precipitación en México. Siendo la precipitación un fenómeno complejo que sólo puede medirse en algunos pocos puntos del territorio se hace una estimación de cuánto llovió en las zonas donde no se mide a través de un mapa de isoyetas. En este capítulo se describe qué son las isoyetas y los motivos que generaron la investigación que se presenta en esta tesis, se muestran las razones por las cuales se consideró necesario hacer un mejor método de interpolación de precipitación y se muestran los problemas que presentan las interpolaciones que se utilizan actualmente. Se describen características de la precipitación. Luego, se muestra el análisis estadístico realizado sobre los datos de precipitación utilizados, que es la base de datos nacional de precipitación de México, con la finalidad de comenzar a conocer la naturaleza de la información que se quiere representar. Se establece la hipótesis de este trabajo y cuáles son los objetivos que se persiguieron y la manera que se siguió para llegar a alcanzarlos. Se ubica al problema en el área de conocimiento correspondiente y se enuncian las materias o áreas de conocimiento que se emplearon para resolver este problema. Se explica para qué se realizó este trabajo y finalmente se explica el contenido de las siguientes partes del documento. 1.1. Contexto y motivación La red de estaciones pluviométricas de Aguas Superficiales de la CONAGUA es de 1,460 estaciones aproximadamente; cada año se incrementan algunas, otras salen de operación y otras se descomponen por lo que el número no es fijo a través del tiempo. Estas estaciones están distribuidas de manera arbitraria en el territorio nacional, es decir, no siguen ningún patrón de distribución regular ni están sujetos a una malla geométricamente uniforme, son lo que se conoce como un conjunto de datos dispersos. Conocer cuánto y en dónde llueve es importante para la economía, la planeación urbana y de producción, y otras múltiples actividades humanas. Importantes decisiones a nivel gubernamental, empresarial y personal, de acción y prospección, se toman diariamente con base en los mapas de precipitación areal. Para las cuencas grandes la distribución espacial de la lluvia se vuelve el factor determinante para descripción útil para estudios hidrológicos [1]. Bajo 2 las condiciones actuales es imposible hacer medición de precipitación sobre todo el territorio nacional, es por ello que con los puntos de medición disponibles se hace una interpolación de una superficie de precipitación areal. Las curvas de nivel sobre un mapa que representan una superficie de agua se llaman isoyetas. La estimación de la precipitación areal para un territorio definido se acostumbra dar como un número escalar el cual representa la altura de una lámina de agua uniforme con milímetros como unidades. El método por antonomasia para el cálculo de la precipitación areal es el de Thiessen, el cual asigna cada valor medido al área definida por su polígono de Voronoi. Con la llegada de las computadoras esta manera ha cambiado y según nuestra propia investigación los métodos matemáticos que utilizan los programas comerciales son los métodos Kriging y Spline de Curvatura Continua Sujeto a Tensión. La importancia de Kriging como la herramienta para la generación de superficies que describen un fenómeno de este tipo está establecido en [2]. Los profesionales del área utilizan estos métodos en los programas comerciales, tales como Surfer© y Matlab©. Estos dos métodos no representan la precipitación de manera adecuada porque sus estimaciones matemáticas no toman en cuenta la geofísica del territorio, como la presencia de un golfo que afecta mucho el comportamiento de la precipitación en la realidad. México es un país que tiene una gran diversidad de climas en su territorio y esto se debe a diferentes factores tales como la ubicación geográfica y la forma del terreno. Hacia el este se tiene la Sierra Madre Oriental como una cadena montañosa alta que conforma una barrera antes de llegar al Golfo de México. De manera similar, hacia el oeste se tiene la Sierra Madre Occidental que también es una cadena montañosa alta que protege el territorio del mar Pacífico. En el Noroeste se tiene la Península de Baja California, que tiene una cadena montañosa a todo su largo y que no es tan alta excepto en algunos puntos. En el Sureste se tiene la Península de Yucatán que es excepcionalmente llana. Uno de los fenómenos conocidos que se presentan en zonas de mar enfrentadas con cadenas montañosas altas, tal y como se presenta en nuestro territorio nacional, es la precipitación al enfrentarse las masas de humedad provenientes del mar con cadenas montañosas altas, provocando lo que en geofísica se conoce como enfriamiento adiabático de la atmósfera. 1.2. Antecedentes Las isoyetas son líneas de contorno que se dibujan sobre un mapa bidimensional de un terreno para representar el nivel de precipitación por área, matemáticamente estimado con base en un conjunto de mediciones puntuales. Las mediciones pueden ser de periodos de tiempo como diarias, mensuales o anuales. La estimación de los valores se obtiene mediante algoritmos de geometría computacional. El método utilizado en meteorología llamado Thiessen [3] calcula una triangulación de Delaunay con los puntos de medición y se establecen sus polígonos de Voronoi; se pondera el área de dichos polígonos sobre el área total asignándoles el valor medido más cercano, y de esta manera se obtiene la precipitación por área total. Este método es adecuado para calcular la precipitación areal cuando se analizan planicies, sin embargo, el efecto de las cuencas hidrológicas naturalmente especificadas por la orografía del terreno provoca imprecisiones grandes para terrenos sinuosos. México tiene gran cantidad de terrenos sinuosos, donde pueden encontrarse ecosistemas diferentes separados por distancias 3 tan cortas como 200 Km. o menos, pero donde se concentra gran cantidad de actividad humana, por lo que se requiere información de mejor calidad. Las isoyetas son una particularidad de una gráfica de contornos, cuando ésta muestra un terreno y sobre ella se representa la precipitación areal. Las líneas de contorno pueden pensarse como cortes a una función escalar con dominio bidimensional, con planos paralelos al dominio. Los algoritmos para su generación se pueden dividir en 2 familias, aquellos que utilizan una malla para determinar los valores interpolados y algoritmos de subdivisión recursiva. Todos los métodos de interpolación presentan divergencias matemáticas al extrapolar, situación que se presenta en el mapa de isoyetas dado que no se puede medir la precipitación en toda la superficie del país. Esto nos antepone dos retos computacionales: el científico que es el desarrollo de un nuevo algoritmo de generación de isoyetas que tome en cuenta la orografía del terreno, y el tecnológico, que consiste en generar implementaciones adecuadas para un rápido procesamiento, haciendo uso cuando así convenga del procesamiento en paralelo. De esta manera se podrán aprovechar las nuevas tecnologías de cómputo disponibles en las computadoras personales. En específico, las actuales capacidades de procesamiento que ofrecen los procesadores gráficos (GPU) de las nuevas tarjetas de video y su disponibilidad masiva, serán el instrumento de aplicación de los resultados de estos trabajos, esperando aportar disponibilidad y precisión al estimar la precipitación por área con base en los puntos medidos. 1.3. De cómo es la lluvia Con la finalidad de conocer los datos que se quieren representar se describen las características relevantes de cómo es la lluvia, cómo se mide y se hace un análisis estadístico básico de los datos. La Organización Mundial Meteorológica (OMM) define lluvia como el fenómeno meteorológico de precipitación de agua líquida en gotas mayores a medio milímetro. La lluvia es un fenómeno local, por ello si se mide cuánto llueve en un solo día y con esa sola medición se quiere hablar de cuánto llovió en todo un territorio el resultado es impreciso. Por ejemplo, en la Ciudad de México es posible observar que en menos de 20 km las condiciones de precipitación son diferentes. Puede estar lloviendo en Ciudad Universitaria pero en Lindavista no a la misma hora. Aún con fenómenos grandes como huracanes, la precipitación no es uniforme en todo el territorio. Sin embargo cuando se habla de precipitación acumulada en periodos largos, del orden de meses o años, entonces se puede hablar de cuánto ha llovido en el territorio con un verdadero sentido de descripción que tiene correlación con la realidad. Es así como un mapa de precipitación areal puede ayudar a describir el fenómeno de lluvia, información que es relevante para muchas actividades humanas. Los mapas de precipitación areal anual generados a partir de las mediciones que se hacen a través de un conjunto de años, están normalizados y muestran el comportamiento medio del fenómeno en cada zona. La OMM establece que para la definición de los valores normales se deben tener criterios de calidad y densidad de datos. Para que cada estación pueda ser considerada adecuada para 4 incluirse en el cálculo se busca tener una serie de datos de treinta años, que comience por ejemplo el primero de enero de 1931 y termine el 31 de diciembre de 1960. La serie de datos de una estación pluviométrica no le deberán faltar más de 3 registros diarios consecutivos y no deben faltar más de 5 datos en un solo mes, según lo estipula el documento NORMAL_WCDP_No10_WMO_No341 de la OMM para el cálculo de normales. Esto es una relación global de 5/6 de datos (83.3% de datos). Usando datos que cumplan estos lineamientos pueden establecerse los llamados valores normales de precipitación. El valor normal de precipitación para el territorio mexicano es de 780 mm anuales. Este es un valor promedio para todo el territorio por lo que en las zonas desérticas llueve menos y en las selváticas más. Para tratar de entender cómo llueve sobre el territorio nacional, se hará una exploración de los datos diarios de precipitación nacional de 1971 a 1980. 1.3.1. Escala de medición. La lluvia se mide en milímetros de altura de agua acumulada sobre una superficie horizontal en un periodo de tiempo. En los pluviómetros manuales la resolución es de 0.1 mm y se acostumbra a hacer una lectura diaria a las 8:00 h, en estaciones automatizadas la escala mínima puede ser 0.25 mm y se pueden tomar lecturas cada 10 minutos. El mínimo valor registrado es cero y el máximo es casi un metro en un periodo de 24 horas. 1.3.2. Suficiencia de medición. Se analizaron estadísticamente las 1,460 estaciones pluviométricas disponibles cuantificando los días del año que se midió la lluvia. Cuando sí hay medición y no llueve el registro dice 0 mm; cuando no se mide no hay registro y no se sabe si llovió o no. Cuando llueve muy poco, brizna apenas perceptible, se registra como inapreciable. Poniendo como criterio de suficiencia de datos el que una estación deberá tener al menos el 85% de los días del año con registro de lluvia, se encontró que solamente 1,307 estaciones tienen 3 o más años completos. Estas 1,307 estaciones son las candidatas a tomarse en cuenta en los cálculos y análisis que se describen en este trabajo. 1.3.3. Área de cobertura de la medición. Es materialmente imposible medir la precipitación sobre toda la superficie del territorio con las tecnologías actuales. Aunque algún día lo fuera, todos los datos históricos permanecerán como un conjunto de mediciones discretas dispersas en una distribución espacialmente arbitraria y con variaciones de periodo a periodo, ya que no todas las estaciones funcionan todo el tiempo y a veces se abren más y otras se cierran. En realidad son pocas las estaciones que se tienen para la vastedad del territorio. El radio de cobertura establece el área de competencia de cada dato, es decir, que cada medición es válida para tomarse en cuenta por los cálculos a la distancia dictada por el radio de cobertura. Esta es una limitación que nos impide para los métodos propuestos en esta tesis, estimar cuánto llovió en Mérida basándonos en las mediciones que se hicieron en Ciudad Victoria, cosa que los métodos globales sí hacen. Es por ello que a cada estación puede asignársele un radio de competencia de la medición. En la gráfica de la Fig. 1 se muestra el área por radio de cobertura de cálculo de todas las estaciones 5 a nivel nacional. La cartografía digital que se maneja está descrita a través de unas polilíneas que definen unos polígonos del territorio nacional. El área contenida en dichos polígonos es lo que se considera el valor real del área territorial. Entonces en la Fig. 1 primero se muestra el área total de los polígonos que conforman el territorio nacional, un poco menos de dos millones de kilómetros cuadrados. Luego, si se forman círculos a partir de cada estación y se les asigna un radio determinado y se suman las áreas de todos los círculos, la cobertura de tal adición es la que se muestra en los diferentes casos. - 500,000 1,000,000 1,500,000 2,000,000 2,500,000 Area polígonos Area 500km Area 200km Area 100km Area 50km Area 15km Radio cobertura en km A re a en k m 2 Fig. 1 Área total nacional que sería cubierta si a cada estación se le asigna un radio de cobertura determinado. Son pocas estaciones para un vasto territorio y el radio de cobertura es el área donde tiene validez cada dato. En los seguros agrícolas que protegen las inversiones económicas hechas en los cultivos agrícolas de desastres naturales como sequías o inundaciones, las compañías de seguros exigen que haya una estación meteorológica certificada a una distancia máxima de 15 km del terreno sembrado. Si a cada estación de medición de precipitación se le asigna un radio de competencia de dichos 15 km no se alcanza a cubrir ni la cuarta parte del territorio, tal y como se muestra en la extrema derecha de la gráfica de la Fig. 1. Tres cuartas partes del territorio quedan cubiertas con radio de cobertura de 50 km y con radios de cobertura mayores o iguales que 200 km queda cubierta la totalidad del territorio nacional, de manera que no quedan huecos sin dato en el mapa del país. Esto es válido en las Regiones Hidrológicas, como lo muestra la Fig. 2. Las Regiones Hidrológicas son las cuencas hidrológicas del terreno, zonas donde el agua escurre hacia el interior y están mayormente delimitadas por las partes altas de las montañas y los cuerpos de agua. Las 37 Regiones Hidrológicas que usamos en los cálculos son las definidas por la CONAGUA. 6 0.0% 20.0% 40.0% 60.0% 80.0% 100.0% 120.0% RH01 RH02 RH03 RH04 RH05 RH06 RH07 RH08 RH09 RH10 RH11 RH12 RH13 RH14 RH15 RH16 RH17 RH18 RH19 RH20 RH21 RH22 RH23 RH24 RH25 RH26 RH27 RH28 RH29 RH30 RH31 RH32 RH33 RH34 RH35 RH36 RH37 Región Hidrológica % d e ár ea t o m ad a en c u en ta e n c ál cu lo Rel.area polig Rel.area 500 Rel.area 200 Rel.area 100 Rel.area 50 Rel.area 15 Fig. 2 Área de cada una de las 37 Regiones Hidrológicas comparada con la parte cubierta por las estaciones si se consideran diferentes radios de cobertura. La cobertura total de los territorios se alcanza con radios mayores o iguales a 200 km. 7 1.3.4. Días con lluvia. La estadística de las 1,460 estaciones pluviométricas disponibles sobre los días con precipitación mayor o igual a 1 mm, denominados aquí como días con lluvia, se realizó para cada estación y para cada año de datos de dicha estación y se calcularon los parámetros de la distribución normal para cada uno. Luego se calcularon los valores de los límites de media con dos desviaciones estándar para verificar las propiedades de la distribución normal esperando cubrir con este intervalo el 95% de los valores. Este análisis muestra que cada una de las series de datos anuales de lluvia no tienen una distribución normal; el 7.4% de ellas dan valores negativos en el límite inferior del intervalo. 1.3.5. Lluvias de un día. En las estaciones del desierto no es raro que en un sólo día llueva arriba del 50% de la precipitación total anual. Para mostrar esto se ha seleccionado, con conocimiento amplio de las estaciones y las zonas geográficas, la estación de San Felipe, B. C. que está en el desierto y es un caso interesante. En la Tabla 1 se muestran las precipitaciones diarias más grandes registradas en dicha estación comparadas con la precipitación anual de ella en la década 1971 a 1980, donde por cierto los datos de 1973 no están disponibles; en un sólo día puede llover más que la precipitación acumulada en varios años. 240 milímetros en un día es equivalente a la precipitación de 7 años en esa misma década. Tabla 1 Precipitación diaria comparada con la anual de una estación de desierto. Puede observarse que no es raro que en un día llueva el 50% de la precipitación total anual; y más que la precipitación anual acumulada en varios años. San Felipe, Baja California Total anual Tres lluvias diarias mayores Tres lluvias diarias mayores en porcentaje 1971 54 24 10 8.5 44% 19% 16% 1972 278.7 240 22 4.5 86% 8% 2% 1974 17 8 7 2 47% 41% 12% 1975 17 11.5 3 2.5 68% 18% 15% 1976 37.8 8 6 5.7 21% 16% 15% 1977 57.4 15.2 10 9 26% 17% 16% 1978 57.5 12 11 6 21% 19% 10% 1979 114 63 16.5 12 55% 14% 11% 1980 40.5 12 10.5 6 30% 26% 15% mediana 54 12 10 6 44% 18% 15% promedio 74.9 43.7 10.7 6.2 44% 20% 12% Pero no solamente eso, aun en zonas de alta precipitación tales como la sierra del norte de Puebla, Veracruz o en Chiapas y bajo fenómenos como ciclón, en unos pocos días o pocas horas puede precipitar un porcentaje importante del valor acumulado anual. Para mostrar esto se ha seleccionado otro caso interesante en los datos, la estación Tierra Blanca, Ver. El fenómeno de alta precipitación en un día en zonas donde llueve mucho puede verse en la Fig. 3, donde se muestra la precipitación diaria para la estación Tierra Blanca en el estado de Veracruz. Se muestran las precipitaciones diarias de junio y julio de 1971, donde se puede ver 8 que 16 días precipitó arriba de los 50 mm y 3 de las precipitaciones fueron mayores a trescientos mm. En estos dos meses la precipitación acumulada fue de 3,377 mm, que representa poco más de la mitad de la precipitación anual de 7,216 mm, para ese año en esa sola estación. Este fenómeno no es regular por lo que matemáticamente puede afectar ampliamente los valores medios. Hay ubicaciones donde típicamente la lluvia de un día es del mismo tamaño casi todos los días que llueve. Esto puede verse en la gráfica de la Fig. 4. donde la estación de la Presa Taxhimay en el estado de Hidalgo, para 1971 tuvo una precipitación anual acumulada de 791 mm y todas sus precipitaciones estuvieron por debajo de los 50 mm. Fig. 3 Precipitación diaria de la estación Tierra Blanca Veracruz de junio a julio de 1971. En este periodo ha caído la mitad (3,377 mm) de la precipitación anual acumulada (7,216 mm). Se tomó la intensidad de lluvia diaria y se comparó de manera porcentual al total anual de precipitación, tal y como se muestra en la gráfica de la Fig. 5. Casi una cuarta parte de las precipitaciones son de uno a diez milímetros, casi otra cuarta parte la conforman precipitaciones que van entre los 10.1 y 20 mm. En esta gráfica se muestra la distribución ocurrida categorizada en intervalos de 10 mm y es interesante observar que la forma de la distribución es similar a la distribución chi-cuadrada. En la Fig. 5 también se muestra la lluvia acumulada por frecuencia de intensidad de lluvia y su gráfica acumulada. Es sumamente interesante descubrir que el 90% de las veces que precipita en México caen de 1 a 60 mm de agua. Es muy interesante que la estadística realizada nos muestre que en promedio la magnitud de la precipitación en las zonas áridas sea muy similar a la magnitud de la precipitación en las zonas pluviosas, sólo que en éstas precipita de 30 a 100 días más que en aquellas. La precipitación máxima registrada en un día a nivel mundial es de 1,825 mm, en enero de 1966, enfrente de Madagascar [4]. 9 Fig. 4 En rojo la precipitación diaria de la estación de la Presa Taxhimay en el estado de Hidalgo, en azul la acumulada. Datos de junio y julio de 1971. Es una zona donde las precipitaciones los días del año que sí hay lluvia se parecen en magnitud mucho unas a otras, todas debajo de los 50 mm. Ahora, si comparamos la intensidad de lluvia media respecto al porcentaje anual se tiene la gráfica azul de la Fig. 5, donde se ve que la precipitación media entre 1 y 10 mm representa más del 20% de la precipitación anual. Algo similar sucede con las precipitaciones medias entre 10.1 y 20 mm. Al hacer la misma comparación pero con el máximo valor de precipitación para una estación, se tiene la línea rosa, donde se ve claramente que el peso relativo de la precipitación máxima fácilmente puede representar arriba del 70% de la precipitación que habrá todo el año, toda esa agua en un sólo día. Es muy interesante observar que la lluvia máxima típicamente representará arriba del 10% de la precipitación anual para el máximo valor de una estación. Fig. 5 Gráfica de la intensidad de lluvia diaria respecto al porcentaje de la lluvia anual total. Porcentaje de precipitación anual respecto a la anual por estación de los datos de la década 1971 a 1980. En rojo se muestra la lluvia acumulada por frecuencia de intensidad de lluvia y en azul el valor acumulado. 10 Fig. 6 Gráfica de la intensidad de lluvia media diaria y de la máxima lluvia diaria registrada en ese intervalo comparadas contra el porcentaje de lluvia anual. 1.3.6. Distribución espacial de las estaciones. La ubicación de las estaciones se muestra en el mapa de la Fig. 7. Las estaciones no tienen una distribución espacial uniforme. Hay lugares entre el estado de Chihuahua y el de Coahuila donde una estación dista más de 400 km de otra más cercana, y en el Distrito Federal se encuentran algunas a menos de dos kilómetros. La red pluviométrica de acopio diario de la CONAGUA está determinada por diversos factores. Mayormente la medición se hace en las zonas pobladas, donde hay personal para tomar y transmitir la lectura diaria en las estaciones convencionales con registro manual en papel. La transmisión de la información se hace por radio, teléfono, o electrónicamente por correo o captura en el sistema. La importancia económica de la región, cantidad de población, incidencia de fenómenos que provocan afectaciones a las poblaciones son otros de los criterios utilizados. 1.3.7. Altitud de las estaciones. Dado que en general las estaciones pluviométricas están ubicadas donde hay población, no hay muchas estaciones en las zonas altas del país. En general están por debajo de los tres mil metros de altura. 11 Fig. 7 Distribución espacial de las estaciones de medición pluviométrica. Las zonas urbanas y económicamente muy activas son las que tienen mayor densidad de estaciones. 1.3.8. Generalidades. Los mapas que se producen en este trabajo son interpolaciones de precipitación en el espacio; se parte de un conjunto de mediciones y la interpolación indica cuánto precipitó en aquellos lugares donde no se midió. Hacen el cálculo de una superficie de precipitación sobre un territorio dado, estableciendo la precipitación areal y sus isoyetas. Los mapas materia de este trabajo no son predicciones de clima ni son modelos termodinámicos que describen el comportamiento atmosférico. Dado lo anterior, nuestros objetivos en el presente trabajo son los siguientes. 1.4. Hipótesis Al tomarse en cuenta la orografía para el cálculo de la superficie de precipitación se puede obtener mejor representación de la realidad de precipitación sobre el territorio nacional. Además esto puede hacerse en tiempos de cálculo del orden de minutos, de manera que estos nuevos algoritmos puedan ser utilizados interactivamente gracias a una implementación de cómputo en paralelo. 12 1.5. Objetivos 1. Desarrollar un nuevo modelo de interpolación de precipitación con orografía de manera que se tome en cuenta el efecto de variación de la precipitación por la altitud para obtener una mejor estimación del valor promedio de precipitación por área. 2. Implementar en paralelo el algoritmo de generación de isoyetas con orografía en una computadora personal con tarjeta de video programable. 1.6. Metodología Para la realización de este trabajo primero se hizo una recopilación de información donde se buscaron los métodos matemáticos de interpolación que se utilizan para resolver problemas con conjuntos de datos dispersos y se investigó cuáles eran los métodos utilizados por los programas comerciales. Se seleccionó el modelo geométrico del algoritmo a utilizar, que fue la 1-Vecindad de Voronoi. Este modelo utiliza toda la información disponible de las estaciones cercanas y las que están a la vista del punto de interés. Se hicieron las implementaciones de los algoritmos de comparación y los propuestos; algunos cálculos se paralelizaron por hilos para múltiples CPUs, como el del trazado de las isoyetas, otros mediante el uso de funciones de la biblioteca de rutinas para resolver ecuaciones de álgebra lineal programada en el lenguaje Fortran llamada LAPACK® logrando por diferentes implementaciones el uso paralelizado por GPUs. Una vez establecidas las herramientas de generación de mapas se hicieron los experimentos para encontrar los parámetros adecuados de radio de competencia de las estaciones o resolución de la malla. Luego se corrieron experimentos para comparar el desempeño de los algoritmos. Finalmente se hizo un análisis de la información y se compararon los desempeños de los algoritmos mediante varias técnicas. También se tuvieron comparativas de procesamiento serial contra paralelo. Los datos que se utilizaron fueron los siguientes: la base de datos climatológica de México llamada SIH que le pertenece a la CONAGUA, cartografía digital y el modelo digital de elevación global llamado Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) de la Agencia Espacial Estadounidense (NASA). Las herramientas de programación utilizadas fueron los lenguajes de programación C++ y Delphi©. Además se usaron dos de las implementaciones de LAPACK®: primero la nativa de ACML©, que es la biblioteca de matemáticas esenciales de la compañía de procesadores AMD© y la segunda una implementación propietaria llamada CULA© para la arquitectura de procesamiento en paralelo de la compañía de procesadores gráficos nVIDIA© llamada CUDA©. El desarrollo fue hecho para el sistema operativo Windows© y corre sobre una computadora personal con o sin tarjeta gráfica programable. Como patrón de comparación del desempeño de los algoritmos propuestos se usaron los mapas de isoyetas de la CONAGUA, los valores medidos en las estaciones a través del método de verificación cruzada, los valores que ofrece el método de cálculo de precipitación areal Thiessen 13 a través de mapas de error y finalmente imágenes de satélite de Google Earth© para hacer una mera comparación cualitativa. Para la realización de este trabajo se emplearon las siguientes materias o áreas del conocimiento de este posgrado: Geometría Computacional, Cómputo Científico en Paralelo, Visualización Científica y Graficación por Computadora. También se emplearon unos pocos conocimientos de Hidrología, Geofísica y Geografía. 1.7. Las partes que conforman a este documento Este trabajo está organizado de la siguiente forma: 1. Introducción. Este capítulo en el cual se explican el contexto y la motivación de la investigación. Además se incluye un breve análisis estadístico de los datos de precipitación en México. 2. Conceptos preliminares. Definiciones de precipitación areal, isoyeta y superficie de precipitación. Métodos de interpolación utilizados actualmente para la generación de mapas: Kriging y Spline de curvatura continua sujeto a tensión. 3. Algoritmos. Detalle de los algoritmos propuestos. Consideraciones de resolución. 4. Implementación. Las partes más importantes de la implementación, preprocesamiento y la parte de paralelización y uso de GPUs para algunos de los cálculos. 5. Resultados. Explicación de los experimentos realizados y detalle de los resultados. Verificación y cuantificación de los métodos propuestos. 6. Conclusiones. Conclusiones y trabajo futuro. 15 2. Conceptos preliminares En este capítulo se hace una breve revisión de los conceptos y terminología alrededor de la precipitación areal. Se revisan las ideas que sirvieron de base para la integración de la orografía a la interpolación de precipitación. Finalmente se describen los métodos de interpolación utilizados en los programas utilizados por los profesionales del área: Kriging y Spline de curvatura continua sujeto a tensión. 2.1. Definiciones Como ya hemos establecido con anterioridad y con la finalidad de tener en esta sección concentradas las principales definiciones, la precipitación areal es una cantidad escalar que dice cuánto llovió en un territorio en un periodo determinado expresándolo como la altura en milímetros de una lámina uniforme de agua acumulada en el periodo y sobre todo el territorio. En esta tesis se hacen mapas de precipitación areal. Dado un conjunto de mediciones dispersas en el territorio los mapas interpolan la precipitación entre las estaciones, en los lugares donde no hay medición. Mapas de isoyetas es la denominación común de los mapas donde se muestra la precipitación areal a través de unas curvas de nivel que representan la precipitación en cada punto del territorio llamadas isoyetas. Estas curvas son una manera de representar la superficie de precipitación subyacente, son cortes de equivalores en ella. Línea de contorno e isolínea son otros nombres que se le dan a las curvas de nivel o curvas de isovalor que van describiendo los lugares de la superficie donde la altura es la misma, la indicada por dicha curva. Las líneas de contorno pueden pensarse como una función escalar con dominio bidimensional, con planos paralelos al dominio, esto se debe a que la superficie de precipitación es una función escalar definida sobre un territorio, un dominio bidimensional de latitud y longitud. Como un ejemplo puede pensarse que es similar a la superficie terrestre, donde cada punto del dominio, posición latitud y longitud, tiene asociada una altitud. Para el caso de la superficie de 16 precipitación la cantidad de precipitación hace las veces de altitud. Es así que se puede definir una función, la isolínea, que corresponde a los puntos en dicha superficie que tienen exactamente ese valor y los puntos que pertenecen a la isolínea se obtienen exactamente cortando a la superficie de precipitación con un plano paralelo al dominio a una distancia definida por el valor de la isolínea, por ejemplo, 500 mm. La precipitación se mide en estaciones pluviométricas o meteorológicas mediante un pluviómetro. En las llamadas estaciones convencionales el pluviómetro puede ser un tanque abierto de área conocida y donde se mide la lámina de agua cada día. En estaciones automáticas el registro puede ser cada 10 minutos. Precipitaciones normales, anuales, mensuales y de otras periodicidades son las que pueden encontrarse comúnmente para el manejo de esta información en los servicios meteorológicos, hidrológicos y estudios geológicos e ingenieriles. El establecimiento, mantenimiento y operación de las estaciones pluviométricas es caro, además de ser físicamente imposible medir en cada punto del territorio. Esto nos deja en la siguiente situación: se tiene un conjunto de mediciones dispersas sobre el territorio, no siempre en la misma posición, y definitivamente no están hechas de manera que sigan una malla regular. Interpolación es un proceso mediante el que se asigna valor a una ubicación en medio de mediciones conocidas. Extrapolación es cuando el punto a asignarle el valor no está entre los puntos conocidos sino fuera de ellos. Dependiendo del área del conocimiento de la que se trate, al mismo proceso matemático de interpolación se le conoce también como predicción, en estadística, o como reconstrucción en imágenes y señales. Como ya se ha mencionado, dada la naturaleza discreta y local de la precipitación, sucede que cuando se genera un mapa de la superficie de la precipitación se está haciendo una interpolación de la precipitación. El mapa con datos de un día o de una semana va a tener imprecisiones grandes. El análisis de la precipitación areal en muy corto plazo no tiene sentido o representación de la realidad por la fuerte localidad de la precipitación. Sin embargo, cuando el periodo es largo, de varios meses en adelante, entonces el efecto de precipitar aquí a veces y allá otras se promedia y comienza a representar la realidad. Según [5], los algoritmos se pueden dividir en 2 familias, aquellos que utilizan una malla fija para determinar los valores interpolados y algoritmos de malla variable de subdivisión recursiva. Todos los métodos de interpolación presentan divergencias matemáticas al extrapolar, situación que se presenta en el mapa de isoyetas dado que no se puede medir la precipitación en toda la superficie del país. A continuación se presenta una lista de algoritmos de las dos familias mencionadas. Familia 1: Algoritmos que usan una malla para generar los valores interpolados.  Tipo 1: Interpolación Spline en mallas regulares. 1. Interpolación polinomial. 2. Interpolación Spline bilineal. 3. Interpolación Spline bicuadráticos. 4. Interpolación Spline bicúbica. 5. Interpolación Spline birraccionales. 17  Tipo 2: Interpolación Spline para puntos distribuidos arbitrariamente. 1. Métodos globales sin triangulación. 2. Triangulaciones. 3. Interpolantes Spline lineales sobre triangulaciones. 4. La aproximación de la primera derivada parcial. 5. Interpolantes Spline cuadráticos sobre triangulaciones. 6. Interpolación Spline cúbica sobre triangulaciones. 7. Interpolación Spline C1 de quinto grado sobre triangulaciones. Familia 2: Algoritmos de subdivisión recursiva.  Algoritmo general de subdivisión interpolatoria sobre triangulaciones arbitrarias. Se denominará superficie de precipitación a la interpolación producida a partir del conjunto de mediciones de las estaciones pluviométricas, que por no estar distribuidas espacialmente de manera regular, pertenecen a la categoría de datos dispersos; es decir, el problema que aquí se trata es el de una interpolación con datos dispersos. 2.2. Mapas actuales En el Servicio Meteorológico Nacional de México (SMN) los mapas son hechos con el programa de graficación llamado Surfer© con un conjunto fijo y predeterminado de estaciones, y es interpolado con el método Kriging. Un ejemplo de estos mapas puede verse en la Fig. 8 donde se muestra la precipitación media anual obtenida con datos entre los años 1941 y 2005. Este mapa fue tomado de la página Web del SMN. Este mapa representa el promedio de la precipitación de todo un año, representada como una lámina acumulada de agua cuya altura se mide en milímetros, el mapa indica que la lámina media anual es de 773.5 mm para todo el territorio nacional. La escala de valores y colores en el mapa se encuentra del lado izquierdo y no es lineal. En realidad la precipitación se mide en sólo unos puntos y luego se estima cuánto precipitó en donde no se midió con algún método de interpolación, Kriging en este caso. Algunas características que pueden observarse en esta interpolación son la continuidad de la superficie a pesar del Golfo de California. Es diferente la cantidad de precipitación en la Península de Baja California que en la parte continental, en los estados de Sonora y Sinaloa. Esta metodología tiene características desventajosas, como decir que la precipitación en la península es un reflejo de la del continente y decir que llueve mucho en zonas áridas y extender demasiado los eventos de precipitación medidos. Estas desviaciones se presentan también en el interior del continente por la interpolación generada. 18 Fig. 8 Mapa oficial del Servicio Meteorológico Nacional, CONAGUA (tomado de su página Web), el cual muestra la precipitación media anual. Está indicado que la lámina media anual es de 773.5 mm. Puede observarse en el Golfo de California que el método matemático de interpolación asume continuidad del terreno y en la Península de Yucatán menor precipitación en la costa. 2.3. De la orografía y la precipitación La orografía y la precipitación se relacionan mediante el fenómeno conocido como enfriamiento adiabático de la atmósfera. Para describir el comportamiento atmosférico se usan unas ecuaciones de termodinámica de la atmósfera. Este fenómeno es estudiado por la geofísica y descrito a través de modelos termodinámicos de la atmósfera. La Universidad de Iowa tiene disponible en línea un simulador interactivo de la termodinámica de la atmósfera [6]. Este simulador es útil para ver y estudiar cómo se comportan las masas de humedad al encontrarse con montañas altas permitiendo estudiar el fenómeno termodinámico ya mencionado. En la Fig. 9 se muestra una simulación de una masa de aire que viaja de izquierda a derecha y se encuentra con una montaña como puede verse en la parte superior izquierda de dicha figura. Este encuentro se da bajo las condiciones iniciales de temperatura, presión barométrica, punto de rocío, tamaño del elemento de simulación atmosférica expresado en una cantidad de altitud, la pendiente de la montaña y la altitud de la que parte la masa de aire. Las condiciones iniciales de la simulación son una montaña de un poco más de 2 km de altura con condiciones iniciales de temperatura 20º C, presión de vapor de agua 16.4 mb y con punto de rocío de 14.4º C. Del lado derecho puede verse el resultado de la simulación. La masa de aire se mueve de izquierda a derecha, choca con la montaña subiendo cada vez más hasta que 19 pasa del otro lado. Primero, en la gráfica de la montaña puede verse que comenzó a llover a una altura cercana al kilómetro y dejó de llover antes de llegar a la punta. Junto a la montaña pueden verse dos gráficas, la primera de izquierda a derecha describe el comportamiento de la presión de vapor respecto a la temperatura; en ella se ve cómo desciende la temperatura (enfriamiento adiabático) y cómo va cambiando la presión de vapor con este proceso. Una vez que la masa de aire supera la montaña y conforme se va alejando de ella, se ve cómo se recupera la temperatura. La segunda gráfica también describe el comportamiento adiabático de la masa de aire en movimiento pero ahora comparándolo respecto a la altitud y también se ve cómo se recupera la temperatura cuando supera la montaña. El punto de rocío es aquel en el que se logra la condensación del agua de la atmósfera y está en función de la presión y la temperatura; las cuales son condiciones atmosféricas. Los trazos azules en las gráficas de la Fig. 9 describen el fenómeno simulado de principio a fin. En la gráfica de temperatura contra la presión atmosférica se muestra cómo se va enfriando la masa de aire. Ese enfriamiento provoca precipitación. Al mismo tiempo en la gráfica de la derecha se describe cómo va enfriándose la masa de aire conforme va aumentando la altitud. Esto nos permite concluir que la precipitación se ve fuertemente modificada por la orografía; al haber montañas altas se favorece la precipitación, sobre todo cuando hay masas de aire caliente provenientes del mar y que chochan con ellas. Este fenómeno sucede en el territorio mexicano por ambos litorales, el del Golfo de México y el del Pacífico, que están flanqueados por la Sierra Madre Oriental y la Sierra Madre Occidental respectivamente, y que son cadenas montañosas altas. Además se sabe que alrededor de una montaña alta hay una gran cantidad de microclimas y gran diversidad biológica derivada de ella. Así, la montaña constituye una contribución importante al agua que recibe la zona y se concluye que conforme sube la masa de nubes aumenta la precipitación. Fig. 9 Se muestra el simulador de la termodinámica de la atmósfera que muestra el efecto del enfriamiento adiabático de la atmósfera. El programa pertenece a la Universidad de Iowa, EUA. En México no toda la precipitación proviene de humedad del mar, hay formaciones de nubosidades locales y dentro del continente. 20 Algo importante que sucede en nuestra atmósfera es que con pequeñas variaciones de altitud se obtienen cambios notorios de la presión atmosférica. Recordemos que conforme aumenta la altitud disminuye la presión atmosférica debido a que la cantidad de atmósfera que se tiene encima es menor. Entonces a nivel del mar se miden 1000 mb mientras que en el D. F., que está a 2,200 metros sobre el nivel del mar (msnm) se miden solamente 760 mb. Se observa una variación del 25% en sólo 2,200 m de altitud. Entonces también se observa que la presión barométrica varía mucho aún con cambios de altitud pequeños, teniendo como consecuencia el que se tenga una termodinámica de la atmósfera muy variada. 2.4. Interpolación por Kriging y Spline CCST Por experiencia personal hemos observado que programas de cómputo como Matlab©, Surfer© y GMT® son las herramientas más comunes que usan los hidrólogos y meteorólogos para la producción de mapas de precipitación a nivel mundial. Por ello nos dimos a la tarea de investigar cuáles son los métodos matemáticos que usan para la interpolación. Dado que son programas que se protegen por propiedad intelectual no se pudieron conseguir los algoritmos exactamente, pero se encontró que los dos métodos de interpolación llamados Kriging y Spline de curvatura continua sujeto a tensión son los estándares más conocidos para la generación de superficies a partir de datos dispersos. Antes de que se extendiera el uso de las computadoras los mapas de precipitación areal se hacían a mano con el método de Thiessen, que es un diagrama de Voronoi. Además se usaba criterio experto que tomaba en cuenta la orografía y condiciones conocidas del clima en las diferentes zonas. Por ejemplo, para las isoyetas y el cálculo de la precipitación areal de México, con datos de 1961 a 1990, se contrató a una empresa la que con un equipo de personas trabajó durante varios meses procesando los datos para obtener los mapas finales. Estos métodos manuales no han sido reproducibles. A continuación se dará una breve explicación de los métodos Kriging y Spline de curvatura continua sujeto a tensión. En esta breve explicación de los métodos Kriging no dará el desarrollo matemático ni las ecuaciones que muestran cómo se interpola la función debido a lo vasto de la teoría que habría que abordar y que no es el objetivo principal de este trabajo, además de ser un tema ampliamente documentado en libros de estadística multivariable y métodos estadísticos para ciencias de la atmósfera como [8]. Sólo se abordará un pequeño aspecto, el variograma, que diferencia los modelos de interpolación seleccionados y las ecuaciones mostradas en la correspondiente explicación de los métodos Kriging son sólo para explicar el variograma. En la breve explicación del método Spline de curvatura continua sujeto a tensión sí se establece la ecuación de la función de interpolación donde a la superficie interpolada se le denomina I , sin embargo la definición formal de I se deja para el siguiente capítulo. Se hace de esta manera para tener separado el trabajo de otras personas que fue tomado como antecedente y la teoría desarrollada en esta tesis. 21 A pesar de que estos métodos se tienen disponibles en los programas comerciales ya mencionados y con la finalidad de poder hacer comparaciones objetivas se programaron estos métodos. Así se pueden tener salidas comparables en escala, color y datos de entrada equiparables. 2.4.1. Kriging Kriging es una familia de métodos de interpolación, también conocida como interpolación óptima, donde se asume que la superficie a interpolar se comporta como un campo aleatorio multivariable estacionario, el cual se conforma por la superposición lineal de las funciones Gaussianas que representan el comportamiento de cada una de las mediciones. Esta familia de métodos nació del trabajo de maestría del sudafricano Daniel Gerhardus Krige en el ámbito de la minería, donde se quiere saber la ubicación de las betas a partir de unas cuantas mediciones. Por el tipo de funciones usadas, la superficie producida por este método es continua y diferenciable. Una función Gaussiana queda especificada al conocer los valores de media y varianza. Para cada función que representa el comportamiento en cada estación pluviométrica habrá que definirlas de algún modo. Las diferentes formas de determinar los parámetros de las funciones Gaussinas determinan las distintas modalidades de Kriging. Los más comunes son los siguientes: 1. Kriging simple. Con media igual a cero. 2. Kriging ordinario. Con media constante y conocida. a. Variograma Gaussiano. b. Variograma Exponencial. c. Variograma Esférico. 3. Kriging universal. La media se especifica con un modelo lineal conocido. Kriging Ordinario estima el valor desconocido de la media (tendencia central) usando una combinación lineal ponderada de los datos disponibles. El Kriging simple se acostumbra seleccionar cuando se desea estimar la distribución de valores dentro de un área y el universal cuando se quiere estimar la media espacial y además los datos tienen una tendencia fuerte que puede ser modelada por funciones simples. Kriging Ordinario fue seleccionado como el método para comparar porque es un método de estimación espacial donde se minimiza el error de varianza, que está basado en cómo están distribuidos los datos y cómo varían unos respecto a los otros. La interpolación por Kriging supone que se tiene un fenómeno que puede ser modelado como si fuera un campo aleatorio multivariable estacionario sobre un espacio M con un vector de posición m . Sea ( )I m la función de interpolación, entonces el variograma teórico puede ser definido como           21 1 2 2 h Var I m I m h E I m I m h             , (2.1) donde , n m h , es decir que el valor promedio de la diferencia al cuadrado de los valores de la propiedad en dos puntos separados por el valor absoluto de la distancia h . Por lo que  es independiente de la ubicación m pero dependiente del tamaño y dirección del vector h . El variograma es una herramienta que cuantifica la correlación espacial. 22 Para obtener interpolaciones por Kriging Ordinario se hace uso de un modelo teórico de variograma, que es el que indica cómo cambian espacialmente los valores. El modelo del variograma teórico es el que da apellido Gaussiano, exponencial o esférico a estos Kriging ordinarios. Para el Kriging ordinario la varianza puede estar especificada de diferentes formas. Una manera de hacerlo es definiendo el variograma que ya definimos conceptualmente y del cual puntualizaremos algunas cosas más. El variograma  h es una estadística descriptiva cuantitativa que puede ser representada gráficamente (ver Fig. 10) en una manera que se caracteriza la continuidad espacial (rugosidad) de un conjunto de datos. 0C también es conocido como el efecto de pedazo y es el valor del que parte la variación, que aunque debería ser cero para valores de 0h , debido a variaciones en las mediciones normalmente se tiene un valor diferente. Intervalo a es la distancia a la que la variación se ha estabilizado. 10 CC  es la meseta o valor final donde la variación se estabiliza. La distancia h es diferencia entre el valor estimado y una observación. Fig. 10 Variograma el cual describe la rugosidad espacial de un conjunto de datos. El variograma experimental puede ser encontrado calculando la varianza de cada punto en el conjunto respecto a cada uno de los otros puntos y graficando las varianzas respecto de la distancia h entre los puntos. El variograma experimental indica cómo varían cada uno de los grados de libertad respecto a los otros, de esta manera, para cada par de puntos del conjunto ,i j se obtiene un elemento del variograma ,i j calculado como la diferencia entre la varianza global 2̂ y la covarianza entre sendos elementos , ˆ i jC , así 2 ˆˆ ij ijC   , (2.2) donde la varianza estimada es 2 0 1 ˆ n i i i w      (2.3) siendo iw los valores medidos de la función a interpolar, 0i la variación respecto al origen o punto de referencia y  la media aritmética del valor de la función entre todos los puntos del conjunto. Una vez que se ha calculado el variograma experimental se selecciona el modelo del variograma, que es una función matemática que modela la tendencia del variograma 23 experimental. Las tres versiones de variograma: Gaussiano, Exponencial y Esférico fueron las seleccionadas para hacer las comparaciones en esta tesis y son las que se muestran a continuación. 2.4.1.1. Kriging Gaussiano. 0 1 2 0 1 , 0 3 1 exp , 0 ij C C si h h C C si h a                     (2.4) El modelo Gaussiano del variograma es el que se muestra en la ecuación (2.4), donde ij es el variograma. La interpolación que produce este método se muestra en la Fig. 11 y en la Fig. 12 en una vista en perspectiva, donde puede observarse que la superficie es muy plana, indicando que precipita de igual manera en la Península de Yucatán que en Durango. Fig. 11 Superficie de precipitación obtenida mediante el método Kriging Gaussiano. El color indica cuál es la precipitación media anual con los datos en la década 1971 a 1980. El color verde significa que llueve igual en la península de Baja California que en la Cuenca Lerma e igual que en Guerrero. Las estaciones tienen una influencia espacial muy pequeña y se ven como picos y depresiones en la lámina de la Fig. 12. En estos mapas se utilizaron datos anuales de la década 1971 a 1980 de las estaciones convencionales. De las 1400 estaciones consideradas se tamizaron a través de un filtro de calidad de datos para que cumpla las normas de la OMM para le generación de datos normales, donde sólo se consideran las series de datos con arriba del 85% de datos existentes y no se permite que haya más de 3 días seguidos sin dato. 24 Fig. 12 Mismos datos que la figura anterior pero vista en perspectiva. Así se observa que el radio de competencia de cada estación es demasiado pequeño. Los parámetros con los que se interpola el método de Kriging son críticos para la calidad del resultado. 2.4.1.2. Kriging Exponencial. 0 1 0 1 , 0 3 1 exp , 0 ij C C h h C C h a                   (2.5) Cuando el variograma experimental se ajusta a un modelo Exponencial, entonces se usa la expresión establecida en (2.5). El valor a es la distancia de variación a la que ya se ha alcanzado la meseta 10 CC  , como puede verse en la Fig. 10. La interpolación de Kriging con el modelo de variograma Exponencial puede verse en ambas, la Fig. 13 y en perspectiva en la Fig. 14. La superficie de precipitación producida por el variograma Exponencial tiene mayor suavidad que las otras y es un mapa mejor comportado, de manera que tiene más sentido el resultado que se produce. En estos mapas la escala a la que se está representando la precipitación areal puede verse del lado inferior izquierdo, donde precipitación cero está representado por el color blanco y el máximo de la escala está representándose con el color azul. La escala va de 0 a 5,000 y tiene una resolución de 50 mm, es decir, hay un cambio de color cada 50 mm de precipitación en el mapa. El rótulo encima de la escala es un identificador único que corresponde al nombre del archivo (png) producido por el programa y que tiene de manera abreviada el método por el cual fue producido, la resolución y el conjunto de datos utilizados. El espacio cartografiado en el mapa va de los 86º a los 118º de longitud oeste y de los 14º a los 34º de latitud norte. La división de paralelos está puesta cada cinco grados. 25 Fig. 13 Interpolación Kriging de tipo exponencial. Este tipo de interpolación por el método de Kriging es el que da resultados con más sentido. Puede verse perfectamente en el mapa las zonas secas y húmedas. Fenómenos que se observan en esta interpolación para este conjunto de datos es la identificación clara de las mesetas centrales del norte que son áridas y desérticas y la distribución de la precipitación en la Península de Baja California y en Sonora. Aún así pueden ser observados problemas, sobre todo en Oaxaca donde aparece como de baja precipitación cuando no es así en la realidad. Fig. 14 Acercamiento y vista en perspectiva de la Interpolación de Kriging tipo Exponencial de la figura anterior en donde se aprecia la continuidad de la superficie y cómo cada estación parece un pico o depresión en ella. 2.4.1.3. Kriging Esférico. El modelo Esférico del variograma es el que se muestra en la ecuación (2.6). La interpolación producida para el mismo conjunto de datos puede verse en la Fig. 15 y un acercamiento e una visión 3D de la superficie puede verse en la Fig. 16. 26 3 0 1 0 1 3 0.5 , , 2 , . ij h h C C h a a a C C h a                   (2.6) Esta interpolación es suave y se ve claramente la influencia esférica de las estaciones y se ve cómo amolda la superficie producida. Tiene como ventajas hacer superficies muy suaves pero tiene la desventaja de producir desviaciones grandes respecto a los valores originales rápidamente. Las vistas 3D de estos mapas tienen otra escala de color porque fueron obtenidas con un programa de visualización de datos cartográficos llamado Global Mapper®. En la Fig. 16 puede apreciarse uno de los problemas que produce la interpolación de Kriging, que bajo ciertas condiciones la superficie de precipitación se vuelve negativa. Todas las zonas azules obscuro en el continente así lo establecen, en medio de la Península de Baja California, al fondo del Mar de Cortés, por la zona del desierto de Altar y en una zona de Chihuahua. Fig. 15 Superficie de precipitación interpolada por el método de Kriging Esférico. Pueden apreciarse detalles como que en medio de la selva Lacandona no hay medición y en todos los mapas se cae la superficie interpolada en esa zona. Físicamente no es posible tener precipitación negativa y si los mapas van a ser utilizados para la toma de decisiones entonces no pueden ser usados de esta manera. Aun cuando no haya suficiente medición tiene que tomarse un criterio experto que tenga sentido con la realidad. 27 Fig. 16 Recorte de la vista en perspectiva de la superficie de precipitación interpolada por Kriging Esférico. Se aprecia claramente cómo las estaciones van dando forma a la interpolación y cómo aberraciones matemáticas producen zonas de precipitación negativa. 2.4.2. Spline de curvatura continua sujeto a tensión Aquí se describe el método GMT® de Spline de curvatura continua sujeta a tensión: Métodos Biarmónico, Laplaciano y Tensión Biarmónico-Laplaciano. En 1990 W. Smith y P. Wessel publicaron en la revista de geofísica la manera de generar las mallas usando Spline de curvatura continua sujeto a tensión [7], que es la base matemática para el programa, también de su autoría, GMT: The Generic Mapping Tools, que por su accesibilidad y utilidad se ha convertido en uno de los programas más usados en el área de geofísica y meteorología a nivel mundial. El método de Spline de curvatura continua sujeto a tensión es precisamente la superficie de curvatura mínima y tiene analogía con una lámina flexible que se hace pasar por todos los puntos. Este método interpola los datos usando una superficie que tiene segundas y cuartas derivadas continuas y minimizando la curvatura cuadrada sobre la superficie total. La superficie producida por este método es continua y diferenciable. La superficie producida se denominará I , y las coordenadas de ubicación de longitud y latitud ,x y respectivamente. De esta manera la notación queda congruente con la usada al explicar la teoría desarrollada en esta tesis. Sin embargo aquí se usan así y no se definen previamente. Las definiciones formales se dejan para el siguiente capítulo para tener separada la teoría desarrollada por otras personas y tomada como antecedente en este trabajo. En el método de curvatura mínima se construye un interpolante con segundas derivadas continuas tal que la integral de la curvatura cuadrada sobre la superficie total es minimizada. Estas ecuaciones pueden solucionarse con un método iterativo de diferencias finitas. Los métodos de generación de malla de curvatura mínima utilizan la norma  22 C I dxdy  . (2.7) Que es una aproximación válida a la curvatura total de I cuando I es pequeña. De esta ecuación se puede derivar la ecuación diferencial 28    2 2 , i i i i I f x x y y     . (2.8) Donde ( , , )i i ix y  son los datos medidos. Las funciones if deben ser seleccionadas de manera que i  conforme    ii yxyx ,,  . Y las condiciones de frontera son: 2 0 I n    (2.9) y  2 0I n     (2.10) a lo largo de las aristas, donde n indica una derivada normal a una arista y 2 0 I x y     (2.11) en las esquinas. Las ecuaciones de las restricciones (2.9), (2.10) y (2.11) son llamadas condiciones de arista libre, y con estas condiciones la ecuación (2.8) tiene una solución única con segundas derivadas continuas y es llamada Spline bicúbica natural. Pequeños desplazamientos de  de una lámina delgada elástica con rigidez constante D, sujeta a un estrés normal vertical q y fuerzas horizontales constantes por unidad de longitud vertical de xxT , xyT y yyT , satisfacen aproximadamente   2 2 2 2 2 2 2 2xx xy yy I I I D I T T T q x x y y                . (2.12) La ecuación de generación de malla de curvatura mínima descrita en la ecuación (2.8) es un caso especial de la descrita en la ecuación (2.12), cuando las fuerzas horizontales son cero y las condiciones de frontera representan momento de doblado cero en las aristas de la ecuación (2.9), cero estrés vertical de corte en las aristas, ecuación (2.10) y cero momento de torsión en las esquinas, ecuación (2.11). Físicamente las funciones if representan las fuerzas Dq de los puntos de carga sobre la lámina elástica. Matemáticamente son los coeficientes en una solución compuesta de una combinación lineal de funciones de Green para la placa elástica debida a las cargas puntuales unitarias. La energía de tensión total almacenada en la placa flexible es aproximadamente proporcional a la curvatura. Para todas las superficies doblemente diferenciables interpolando los datos la superficie de curvatura mínima almacena la mínima energía de tensión. Si uno se imagina doblando una placa delgada se requerirá trabajo extra para obtener la superficie de curvatura mínima. La solución de curvatura mínima puede parecer la manera natural de estimar una malla de valores débilmente restringidos, sin embargo está demostrado que puede crear oscilaciones extrañas. Derivamos nuestra ecuación para investigar la función del estrés de tensión isotrópico uniforme T en la ecuación (2.12). Suponga que xxT = yyT y que xyT = 0, entonces la ecuación (2.12) se vuelve  2 2 2 A I T I q     . (2.13) Cuando 0T  entonces la ecuación (2.13) se convierte en la (2.8), pero para una T arbitrariamente grande la solución es dominada por el segundo término. Aquí T tiene unidades de fuerza por unidad de longitud y T requiere un ajuste de la solución escalada con la constante de rigidez flexible A y con el valor escalar q . Se puede evitar esto reescribiendo de la siguiente manera 29      2 2 21 , I I i i i i T I T I f x x y y        . (2.14) La ecuación de interpolación por el método de Spline de curvatura continua de los valores I a partir de un conjunto de i puntos  , , i i i x y  , puede expresarse como está indicado en la ecuación (2.15) para cualquier punto  , ,x y  . De esta manera se tienen dos términos conjugados por el factor de tensión T , uno es el Laplaciano 2 I y el otro el Biarmónico  2 2I  .      2 2 21 ,i i iT I T I f x x y y        (2.15) Así, cuando 1T la superficie de precipitación interpolada sólo contiene el comportamiento del término Laplaciano y cuando 0T sólo contiene el comportamiento Biarmónico. Para variables geofísicas los autores recomiendan un factor de tensión 7.0T . Respecto a las características de la interpolación producida bajo dichas diferentes condiciones tenemos que para 1T el mapa resultante es el de la Fig. 17, donde se puede apreciar cómo la lámina se deprime uniformemente hacia el noreste y es medianamente alta en el sureste. Además, las estaciones parecieran tener poco peso individual, cosa notoria es las estaciones con registro de abundante precipitación en el sur del país donde se ven como pequeñas protuberancias. Ahora, con 0T y sólo el comportamiento Biarmónico se muestra en la Fig. 18, donde pueden verse una zona de precipitación negativa entre los estados de Oaxaca y Guerrero, en una zona donde no hay medición. Físicamente no tienen sentido las precipitaciones negativas y si el mapa se quiere para la toma de decisiones, como cuando la CONAGUA determina si los municipios van a ser sujetos a compensaciones del gobierno federal por eventos como la sequía o las inundaciones, entonces este mapa con desviaciones tan evidentes no puede ser utilizado. Se han mostrado estos dos mapas con condición extrema del factor de tensión para hacer evidente los comportamientos que siempre presenta este método, pero que no es tan fácil identificar. Cuando se utiliza el valor recomendado del valor de tensión 7.0T produce una interpolación mejor comportada, se ve más suave y tiene más sentido respecto a la experiencia personal, tal como se muestra en Fig. 19. 30 Fig. 17 Spline de Curvatura Continua Sujeto a Tensión con valor T=0, llamado Laplaciano. Se aprecia la lámina uniformemente deprimida en el noroeste y medianamente alta en el sureste. Las estaciones parecieran tener poco peso individual y más como conjunto hablando espacialmente. Fig. 18 Método de interpolación de precipitación Spline de Curvatura Continua Sujeto a Tensión con valor T=1, llamado Biarmónico. Presenta una zona de precipitaciones negativas entre Oaxaca y Guerrero, en una zona sin medición. 31 Fig. 19 El mejor desempeño del método de Spline de Curvatura Continua Sujeto a Tensión se da con el parámetro recomendado de T = 0.75. 33 3. Algoritmos En este capítulo se hace la definición matemática formal del problema a resolver, se precisan y describen partes y términos utilizados. Luego se describen los nueve algoritmos de interpolación que se están proponiendo. Finalmente se describe la manera en que se hizo el cálculo de la constante de orografía. Los algoritmos que se describen en este capítulo son: 1. Vecinos Voronoi. 2. Vecinos Voronoi con Orografía. 3. Cuatro Estaciones Cercanas por el inverso de la distancia. 4. Dieciséis Estaciones Cercanas por el inverso de la distancia. 5. Cuatro Estaciones Cercanas por el inverso de la distancia al cuadrado. 6. Dieciséis Estaciones Cercanas por el inverso de la distancia al cuadrado. 7. Cuatro Estaciones Cercanas con Orografía. 8. Dieciséis Estaciones Cercanas con Orografía. 9. Thiessen con Orografía. 10. Trazado de isoyetas con marching triangles método recursivo y paralelizable. 3.1. Definición del problema Sea U el universo de todas las mediciones diarias de precipitación disponibles donde , ( , , , , )i i i i i i id U d x y z p t  (3.1) donde Xxi  donde X es el conjunto de longitudes (posición geográfica),  180, 180X   (3.2) Yyi  donde Y es el conjunto de las latitudes (posición geográfica),  90, 90Y   (3.3) 34 Zzi  donde Z es el conjunto de las altitudes (posición geográfica), el dominio de Z son los números reales Z (3.4) se considera altitud cero el nivel medio del mar. Las ubicaciones que están por encima de ese nivel son positivas, y negativas las que están por debajo. Ppi   , 0, i p P P   (3.5) es el conjunto de posibles valores de precipitación y Tti  , donde T es la fecha de la precipitación diaria de manera que / /it dd mm aaaa (3.6) Sea  el intervalo a consideración definido como un subconjunto de T , T , ,ini finT fecha fecha       (3.7) el intervalo de tiempo sobre el cual se tomarán los datos del conjunto D del universo de datos disponibles U. De esta manera definimos al conjunto D como  ( , ) i ii D U D d d U y t     (3.8) es decir, el conjunto de m datos dispersos  mddddD ,,,, 321  donde ( , , , , )i i i i i id x y z p t (3.9) llamaremos a D conjunto de datos en el intervalo a consideración. Ahora, sobre T se definen los siguientes subdominios a los que llamaremos periodos  y para definirlos se usarán la variables mudas k como elementos de numeración y por lo tanto subíndices para indicar un elemento válido del conjunto representando por ejemplo kaaaa un año cualquiera entre los años válidos. Para la definición de estos periodos se requiere establecer una fecha inicial del periodo y una fecha final, para ello la fecha se representará en el formato numérico / /dd mm aaaa donde dd es un número de dos posiciones que representa el día, mm es un número de dos posiciones que representa el mes y aaaa es un número de cuatro posiciones que representa el año, de manera que el primero de junio del 2011 se representaría 01/06/2011. Además, se usará la variable findemes para representar el último día del mes correspondiente, ya que es diferente dependiendo del mes de que se trate y si es o no bisiesto el año. Entonces, usando la definición de la fecha de precipitación diaria i t en (3.6) un periodo  se puede definir como un intervalo continuo de n fechas que inicia en un i t dado y termina en El día del mes se representará en las siguientes definiciones como una cifra de dos posiciones de manera que , así: El periodo  anual es  01/ 01/ ,12 / 31/ anual k k aaaa aaaa  kaaaa (3.10) El periodo  mensual es 01/ / , / /mensual k j k jmm aaaa findemes mm aaaa     jk aaaamm , (3.11) El periodo  de secas es  101/12 / ,30 / 04 / secas k k aaaa aaaa  kaaaa (3.12) El periodo  de aguas es 35  01/ 05 / ,30 /11/ aguas k k aaaa aaaa  kaaaa (3.13) El periodo  de intervalo arbitrario es , con 1 añoarbitrario inicial final arbitrariot t     (3.14) donde  define el tamaño del periodo. Sea  el periodo sobre el que se normalizarán los datos. Sea c la constante de calidad de la que se habló en la sección que describe la lluvia del capítulo de introducción 1.3 y los criterios de la OMM para la definición de los valores normales, donde 0.85c  (3.15) Sea f la función de calidad definida sobre D donde:  , , i i i el número total de t existentes sobre el periodo f D c E d d D y c             (3.16) por lo que E es un conjunto de n datos dispersos  1 2 3, , , , n E d d d d (3.17) donde ),,,,( iiiiii tpzyxd  . Llamaremos a E el conjunto de datos de buena calidad. YX  es el dominio de las posiciones geográficas viendo a la tierra como una esfera, plana en su superficie. En el dominio definido por YX  se pueden ubicar las estaciones de medición de precipitación o pluviómetros. Por otro lado tenemos que ZYX  es el dominio de las posiciones geográficas considerando la orografía de la superficie terrestre. El dominio definido por ZYX  también es un dominio de ubicación de las estaciones. Cada id representa una ubicación con la medida que se hizo en ella en un periodo de tiempo determinado. Todos los m datos del conjunto D corresponden a un mismo intervalo de tiempo a consideración  , y por la manera que está definido, también todos los n datos del conjunto E . Sea  la regla de acumulación para un periodo  definido para datos de precipitación como la suma de todos los datos disponibles en el periodo  . Sea   ( , , ) , , ,i i i i i ig E B b b x y z p     (3.18) la función de acumulación de datos se define como E bajo la regla de acumulación  sobre el periodo  , y donde ip es el valor de precipitación acumulado bajo la regla de acumulación  en el periodo  , por lo que B es un conjunto de n datos dispersos  nbbbbB ,,,, 321  . Sea  11 21 1 12 22 2 1 2, , , , , , , , , , , , r r s s rs v v v v v v v v v  (3.19) un territorio definido por un conjunto de s polígonos cerrados definidos cada uno por jr vértices, donde ),( ijijij yxv  pertenecen al dominio YX  y definen áreas no necesariamente contiguas sobre él. Sea z una teselación de  a la que llamaremos zonas del territorio, donde  11 21 1 12 22 2 1 2, , , , , , , , , , , , z u u v v uw v v v v v v v v v  (3.20) 36 es un conjunto de w polígonos cerrados definidos cada uno por jw vértices, donde ( , )ij ij ijv x y (3.21) pertenecen al dominio YX  y definen áreas no necesariamente contiguas sobre él. Sea M la malla de resultados definida como un lattice regular rectangular de resolución  0.1, ,0.9,0.01, ,0.09,0.001, ,0.009res  (3.22) indicando la correspondiente fracción de grado (posición geográfica) a la que están definidos los puntos a partir del rectángulo definido por los puntos extremos   iniminiminim yxp , y  , m fin m fin m fin p x y (3.23)  , ,mini m finM M res p p (3.24) donde cada punto de cálculo ji, sobre la malla M tiene coordenadas  resjyresix iniminim  , . Para México  00.14,00.114inimp y  00.33,00.86finmp . Una malla de baja resolución para representar el mapa de México podría ser de un décimo de grado, lo que nos da una malla de 321 columnas por 191 renglones, 61,311 puntos y representaría cada celda aproximadamente 10 km por 10 km. Una de alta resolución podría ser de un centésimo de grado 3,201 columnas por 1,901 renglones, 6,085,101 puntos en la malla y representaría un kilómetro cuadrado cada punto. Una vez definido el tamaño de la malla se ha de calcular el valor de la función a interpolar en cada uno de los puntos que define. Sea ijC la celda correspondiente a las posiciones del punto de cálculo ji, en la malla M donde  0, /m fin minii x x res    (3.25) y  0, /m fin minij y y res    (3.26) Cada ijC tiene una correspondencia a un área geográfica sobre el dominio YX  . El problema a resolver es interpolar una superficie de precipitación sobre el territorio  a partir de B . Las funciones de interpolación  ,I B  (3.27) producen un valor escalar sobre todos los puntos del territorio  definidos sobre una malla M de resultados, estimando cuánto precipitó en puntos donde no se midió.   ( ) , 999 ij ij I B cuando C I B dato no existente cuando C        (3.28) 37 Sea R el radio de influencia de cualquier estación ib , definido como la distancia geodésica Dis máxima a la que será tomada en cuenta el valor ip para las funciones  BI que así lo determinen. Sea Dis la distancia geodésica entre dos puntos 1b y 2b definida como  1 2 1 2 1 2 1 2 2 2 2 2 2 2 ( , ) arccos cos cos cos sin sin 360 360 360 360 360 360 Dis b b r y y x x y y                                       (3.29) donde el radio 249145.6278r que es el semieje mayor del elipsoide Clarke 80 y suponiendo una tierra esférica. En este trabajo se proponen nueve funciones de interpolación  BI y se comparan con las tres de uso más extendido (Kriging, Spline CCST y Thiessen). 3.1.1. Definición de precipitación diaria La precipitación diaria se define como el acumulado de precipitación en un periodo de 24 horas. Para estaciones convencionales oficialmente según criterios de la OMM, el corte de este periodo se hace a las 8:00 de la mañana hora local. Para estaciones automáticas, que toman varias mediciones al día, por ejemplo cada 10 minutos o una hora, se puede tomar otro criterio de acumulación, por ejemplo que el corte se haga a las cero horas. 3.1.2. Datos dispersos Se considera que los datos son dispersos si no se encuentran ubicados espacialmente siguiendo el patrón regular definido por un lattice. Datos que siguen un patrón regular son por ejemplo los que se adquieren al tomar una fotografía digital. Aunque es cierto que para ciertos tipos de fotografías y dependiendo de las condiciones del lente utilizado pueden existir aberraciones que provoquen que el muestreo sobre la escena no haya sido regular, la representación de las intensidades de luz registradas en renglones y columnas se manejan comúnmente como que siguen un patrón regular. Esto es tan cierto que cuando se requiere precisión en las fotografías, como cuando son tomadas desde un avión, se someten a un proceso llamado ortorrectificación que corrige las desviaciones por movimiento, lentes y algunas de estas aberraciones conocidas para finalmente obtener algo más cercano a un muestreo regular sobre la malla. Como ya se dijo las estaciones de medición están ubicadas de manera arbitraria, principalmente en las poblaciones donde personal puede tomar lecturas y dar mantenimiento. En las zonas donde la actividad económica es mayor, hay más medición. Entonces los datos con los que se trabajan en esta tesis son del tipo disperso. 38 3.1.3. Cálculo de malla En el proceso de generación de la solución sobre la malla M con la función  ,BI primero se hace el siguiente proceso llamado de bloque mínimo: Si dentro del área geográfica que define la celda ijC queda más de una estación lb , entonces se sustituyen en B las k estaciones que quedan dentro de ijC por 1 1 1 1 1 1 1 1 , , , k k k k l l l l l l l l b x y z p k k k k                (3.30) De esta manera nos estamos ajustando al criterio Bessel y Smith. La manera de determinar a cual celda pertenece lb es mini l mini x x i x         (3.31) mini l mini y y j y         (3.32) 3.1.4. Puntos de cálculo Para la malla M con celdas ijC , son todos los ijC que quedan dentro del territorio  . La determinación de los ijC que están fuera y dentro del territorio  se hace como un preproceso que depende de la malla M , y se almacena en un archivo externo con extensión (*.pim). 3.1.5. Distancia geodésica La distancia geodésica está definida como la distancia entre dos puntos sobre la superficie terrestre, suponiendo que la tierra es una esfera plana, por lo que afectan al cálculo la latitud de los puntos entre los que se quiere medir la distancia. No afecta la longitud porque está definida como diámetros de la tierra que pasan por los polos, a diferencia de la latitud que está definida como rebanadas paralelas al ecuador, por lo que se hacen más pequeñas conforme se acercan a los polos. De hecho en las latitudes extremas 90 y -90 la rebanada mide cero, es el punto extremo de la esfera. La distancia geodésica entre dos puntos, definida en la ecuación (3.29), fue calculado con el radio 249145.6278r que es el semieje mayor del elipsoide Clarke 80 y suponiendo una tierra esférica. 39 3.1.6. Determinación de interioridad de un punto a un polígono cerrado Para cuando los algoritmos presentados necesitan determinar si un punto está dentro de un polígono de Voronoi, Región Hidrológica o en general cualquier territorio  una parte de él denominada anteriormente como zona z , se utilizó al algoritmo llamado trazado de rayo. El punto del cual se quiere determinar la interioridad al polígono puede ser un punto de la malla ijC o una estación ib . El rayo se traza a partir del punto de interés y luego se cuenta el número de veces que interseca al polígono cerrado kV que representa a dicho  o z en cuestión. Si la cantidad de intersecciones es non, la celda es interior. Para facilidad de cálculo se seleccionó el rayo de x constante. Si interseca con un extremo de segmento se cuenta una vez. Si el punto está sobre la línea del polígono entonces es interior. Recordemos que el territorio  está definido como en la ecuación (3.19), es decir como un conjunto de s polígonos cerrados definidos cada uno por jr vértices, donde ),( ijijij yxv  pertenecen al dominio YX  y definen áreas no necesariamente contiguas sobre él. Las zonas interiores al territorio están definidas como z una teselación de  a la que llamaremos zonas del territorio, donde (3.20) es un conjunto de w polígonos cerrados definidos cada uno por jw vértices, donde ),( ijijij yxv  pertenecen al dominio YX  y definen áreas no necesariamente contiguas sobre él. 3.1.7. Puntos interiores El Algoritmo - 1 es el que determina cuándo un punto es interior. ALGORITMO PUNTOS INTERIORES Algoritmo - 1 ENTRADAS: PUNTO  ,P x y , Y POLÍGONO CERRADO kV . SALIDAS: VALORES DE VERDAD RESPECTO INTERIORIDAD. 1. Cada polígono tiene unos valores extremos mínimo en x y máximo en x . También se tienen los valores mínimo y máximo en y . 2. Ahora, si el punto es mayor que el máximo o menor que el mínimo entonces está fuera. 3. Si la latitud del punto es mayor que la latitud del máximo o la latitud del punto es menor que latitud mínima, entonces está fuera. 4. Para cada segmento del polígono otra vez le aplicamos el criterio, pero ya no sobre puntos extremos del polígono sino 40 por el rectángulo definido por el segmento. 5. Se estableció que el rayo para determinar la interioridad va hacia arriba. 6. Si la latitud de Cij es más grande que la latitud máxima del segmento ks entonces quedó fuera, no analizar. 7. Si la latitud de Cij es más pequeña que la latitud mínima del segmento ks entonces quedó fuera, no analizar. 8. Si la longitud de Cij es menor que la longitud mínima de ks no analizar. 9. Si la longitud de Cij es mayor que la longitud máxima de ks no analizar. 10. Como el rayo va hacia abajo y parte de Cij , entonces si la latitud Cij es menor que la latitud mínima de ks queda fuera. Solamente que la latitud mínima de ks sea menor de la latitud de Cij debo analizar. Entonces veo si corta en x. Si la longitud de Cij es menor que la longitud máxima ks , y la longitud de Cij >= que la longitud mínima ks , entonces cruza. Dos cruces, número par, está fuera. Cero cruces, está fuera. Sólo en número non hace cruce. 11. Una vez que cumplió en quedar en el rectángulo de ks se verifica si realmente cruza utilizando las ecuaciones paramétricas de la recta. Si realmente cruza se acumula el número de cruce. Si te dice que está exactamente sobre el segmento ks entonces es interior y ya no se hacen pruebas. 12. Si la longitud mínima de ks es menor o igual que la longitud de Cij y es menor que la longitud máxima de ks entonces interseca. 3.1.8. Criterio de cercanía Para tener certidumbre y convergencia en el criterio de igualdad en los algoritmos de interioridad, que puede verse afectada por la variabilidad intrínseca de la representación de números flotantes en la computadora, se estableció un criterio de cercanía cumplido el cual establece la igualdad de dos elementos. Se estableció una tolerancia, lo que comúnmente se conoce como delta diferente para la longitud y la latitud. De esta manera si la delta de longitud 41 entre la celda y el elemento del polígono, entonces el vértice es mayor que la distancia máxima de cobertura de la estación por lo tanto se considera fuera de él. 3.1.9. Cálculo de la distancia Y se busca que sea menor que el radio máximo de competencia. Cuando la distancia es adecuada se inserta en la lista de estaciones cercanas. El cálculo de la distancia es geodésica, no euclidiana. Se calcula cuando hay al menos una estación cercana. Cuando hay sólo una estación cercana estará funcionando localmente como Thiessen. 3.1.10. Determinación de la diferencia de altitud A partir de los archivos DEM de la orografía de radar satelital obtenidas de la NASA, se generó un modelo orográfico en un archivo de extensión *.eip con resolución espacial de cada segundo y resolución de altitud de a metro. Sea   iji Cb , la diferencia de altitud entre la estación ib y el punto en la malla ijC . Como cada celda en la malla representa a un punto de cálculo y ese punto de cálculo es una ubicación geográfica real, entonces de los archivos preprocesados (*.eip) se lee la altitud iz correspondiente. 3.2. Vecinos Voronoi Este es un método geométrico. Cada punto en la malla cae en un polígono de Voronoi del diagrama de Voronoi definido por las estaciones. Cada polígono de Voronoi tiene una 1-vecindad de polígonos de Voronoi. En este método el valor asociado en cada punto en la malla es la contribución ponderada por el inverso de la distancia de las estaciones que definen la 1- vecindad más la del propio polígono. Dado el conjunto de estaciones B puede calcularse el diagrama de Voronoi )(BVor donde cada polígono Voronoi iPV define el área más cercana a cada estación. Para cada iPV se definen a ks como las k aristas de dicho polígono. Para cada iPV se definen como sus primeros vecinos al conjunto de polígonos Voronoi jPV que comparten una arista con iPV . 42 Fig. 20 Para el punto en la malla representado por el cruce de las líneas rojas en blanco se muestra el polígono de Voronoi al que pertenece y en gris se muestra la primera vecindad o también llamada 1-vecindad de Voronoi. Cada punto de cálculo en la malla ijC es interior a un polígono de Voronoi iPV . A cada iPV le corresponde una 1-vecindad definida por el conjunto de los jPV que comparten una arista con iPV . Este conjunto de estaciones correspondientes son las que sí son tomadas en cuenta para el cálculo de  BI VoronoiVecinos . Sea BVV  donde ( ) i j ij i i i b es una estación que define a un polígono de PV VV VV C b ó b es la estación que define a PV             . (3.33) Sea VVn 2 el número de estaciones ib en VV . Para México en particular, nosotros encontramos vecindades de hasta 30 estaciones. Entonces    2 2 1 1 1 ( , ) 1 ( , ) n k k k ij Vecinos Voronoi ij n k k ij p D b C I B C D b C       (3.34) , para VVbk  . 43 ALGORITMO VECINOS VORONOI Algoritmo - 2 ENTRADAS: DATOS B, Y MALLA M. SALIDAS: VALORES DE LA SUPERFICIE DE PRECIPITACIÓN PARA CADA ijC . Se calcula el diagrama de Voronoi de B. Para cada ijC se determina a cuál iPV pertenece. Se determina jPV la primera vecindad de iPV . Se determina el conjunto de estaciones VV . Se determina el valor para ijC con la ecuación (3.34) . 3.3. Vecinos Voronoi con Orografía Este es un método geométrico. Cada punto en la malla cae en un polígono de Voronoi del diagrama de Voronoi definido por las estaciones. Cada polígono de Voronoi tiene una 1-vecindad de polígonos de Voronoi. En este método el valor asociado en cada punto en la malla es la contribución ponderada primero por el inverso de la distancia de las estaciones que definen la 1- vecindad más la del propio polígono y ponderado también por el factor de diferencia de altitud entre el punto en la malla que se está calculando y la estación correspondiente al valor a ponderar. Se determinan iPV , jPV y VV como se ha determinado en Vecinos Voronoi con la ecuación (3.33). Entonces si  , 1k ijb C            2 2 1 1 1 1 , ( , ) 1 1 , ( , ) n k k ij k k ij VV Orografía ij n k ij k k ij p b C D b C I B C b C D b C             (3.35) para VVbk  . Cuando  , 1k ijb C   se hace  , 0.9999k ijb C   y se usa (3.35). ALGORITMO VECINOS VORONOI CON OROGRAFÍA Algoritmo - 3 ENTRADAS: DATOS B, Y MALLA M. SALIDAS: VALORES DE LA SUPERFICIE DE PRECIPITACIÓN PARA CADA ijC . 1. Se calcula el diagrama de Voronoi de B. 2. Para cada ijC se determina a cuál iPV pertenece. 3. Se determina jPV la primera vecindad de iPV . 4. Se determina el conjunto de estaciones VV . 5. Se determina el valor para ijC con la ecuación (3.35). 44 La determinación de los valores es diferente a Thiessen. 3.4. 4 / 16 Estaciones Cercanas por el inverso de la distancia y por el inverso de la distancia al cuadrado. En estos cuatro métodos el valor de cada punto en la malla se calcula sumando una contribución de las cuatro o dieciseis estaciones más cercanas pesada por el inverso de la distancia entre el punto en la malla y la estación correspondiente o pesada por el inverso de cuadrado de la distancia. Este método produce una superficie continua y suave, con pendientes menores que el método del inverso del cuadrado de la distancia. Para abreviar Estaciones Cercanas se utilizará en su lugar la palabra Cercanos. Sea  1 4 |16n  el número de estaciones ib en consideración por el criterio. Ahora, sea 1 1 1 1 1 , 2 , _ ( ) _ n l l o ij l ij l b B y b es el más cercano más cercano n Cercanos C n C b a C n ésimo más cercano                (3.36) Entonces    1 1 1 1 _ 1 1 ( , ) 1 ( , ) n i i i ij n Cercanos INV distancia ij n i i ij p D b C I B C D b C       (3.37) donde Cnbi 1 . ALGORITMO N1_CERCANOS INV DISTANCIA Algoritmo - 4 ENTRADAS: DATOS B Y MALLA M. SALIDAS: VALORES DE LA SUPERFICIE DE PRECIPITACIÓN PARA CADA ijC . Para cada ijC se determina el conjunto Cn1 . Para cada elemento en Cn1 se calcula el valor con la fórmula (3.37). 45 También    1 2 1 1 1 2 1 1 ( , ) 1 ( , ) n i i i ij n Cercanos INV distancia cuadrada ij n i i ij p D b C I B C D b C       (3.38) donde Cnbi 1 . ALGORITMO N1-CERCANOS INV DISTANCIA CUADRADA Algoritmo - 5 ENTRADAS: DATOS B Y MALLA M. SALIDAS: VALORES DE LA SUPERFICIE DE PRECIPITACIÓN PARA CADA ijC . Para cada ijC se determina el conjunto Cn1 . Para cada elemento en Cn1 se calcula el valor con la fórmula (3.38). Estas funciones de interpolación  BI definen un subconjunto de estaciones más cercanas al punto de cálculo mediante la ecuación (3.36). Así, aunque la expresión está definida de manera general aquí se particulariza a sólo dos casos: cuatro y dieciséis estaciones más cercanas. Conforme se va cambiando de punto de cálculo dicho subconjunto de estaciones también varía. Entonces se toman los valores de las estaciones seleccionadas por (3.36) y para cada punto en la malla se hace el cálculo del valor con la ecuación (3.37) cuando se trata de que esté ponderado por el inverso de la distancia y con la ecuación (3.38) cuando se trata de que esté ponderado por el inverso de la distancia al cuadrado. La diferencia en ambos tipos de ponderaciones consiste en el grado de concentración de la influencia a partir del punto de medición, de manera que cuando es ponderado por el inverso de la distancia la superficie interpolada tiene pendientes mayores, pareciéndose más a las superficies producidas por el método Thiessen. Al ponderar por el inverso de la distancia al cuadrado la superficie que produce la interpolación es más suave. 3.5. 4 / 16 Estaciones Cercanas con Orografía En estos dos métodos el valor de cada punto en la malla se calcula sumando una contribución de las cuatro o 16 estaciones más cercanas pesada por el inverso de la distancia entre el punto en la malla y la estación correspondiente y ponderada también por un factor de altitud. Se determina el conjunto ijC como lo determina (3.37) en Algoritmo - 4. Sea   iji Cb , la diferencia de altitud entre la estación ib y el punto en la malla ijC . Como cada celda en la malla representa a un punto de cálculo y ese punto de cálculo es una ubicación geográfica real, 46 entonces de los archivos preprocesados (*.eip) se lee la altitud iz correspondiente. Sea  el factor de precipitación por altitud que describe cuánto se modifica la precipitación en función de la diferencia de altitud por unidad de altitud. Entonces          1 1 1 1 1 1 1 , ( , ) 1 1 , ( , ) n i i ij i i ij n Cercanos Orografía ij n i ij i i ij p b C D b C I B C b C D b C             (3.39) donde Cnbi 1 . 3.6. Thiessen con Orografía En este método a cada punto en la malla se le asigna el valor de la estación más cercana ponderado por el inverso de la distancia del punto en la malla a la estación y ponderado también por un factor de altitud. El diagrama de Voronoi de las estaciones define el polígono de cada una, los puntos en la malla que están dentro de dicho polígono son los que tienen el valor ponderado de la estación que lo definió. Se determinan iPV , jPV para identificar el polígono de Voronoi al que pertenece cada punto de cálculo en la malla. Entonces    2 1 2 1 1 ( , ) 1 ( , ) n k k k ij Thiessen Orografía ij n k k ij p D b C I B C D b C         (3.40) ALGORITMO THIESSEN CON OROGRAFÍA Algoritmo - 6 ENTRADAS: DATOS B, Y MALLA M. SALIDAS: VALORES DE LA SUPERFICIE DE PRECIPITACIÓN PARA CADA ijC . 1. Se calcula el diagrama de Voronoi de B. 2. Para cada ijC se determina a cuál iPV pertenece. 3. Se determina el valor de ijC con la ecuación (3.40). 47 3.7. Cálculos areales Como ya se mencionó la precipitación areal es la cantidad precipitada acumulada como una lámina uniforme de agua cuya altura se mide en milímetros sobre un territorio  o una parte del territorio z . La nombraremos PA ó z PA . Se calcula acumulando el valor de las celdas ijC que son interiores al territorio  ó z . ALGORITMO DE CÁLCULOS AREALES Algoritmo - 7 ENTRADAS: LA MALLA CON LOS VALORES DE LA SUPERFICIE INTERPOLADA Y LOS POLÍGONOS DEL TERRITORIO. SALIDAS: CANTIDAD ESCALAR DE LA PRECIPITACIÓN AREAL CALCULADA. Para cada ijC y si ijC en el caso de cálculos areales en todo el territorio o zijC  en caso de cálculos areales en las zonas del territorio: 1. Si es interior se incrementa el contador de celdas para  ó z . 2. Se acumula el valor de la celda por el coseno de la latitud  ii yp cos . 3. Se acumulan los cosenos de las latitudes  iycos . 4. La precipitación areal en  ó z que se define como PA ó z PA va a ser igual a: Para cada elemento de la malla. Si es interior a nuestra área de interés (mapa) determinar para cada zona del mapa de cálculo lo siguiente: 1. Si es interior a la zona (cuenca, polígonos de cuenca o de estado o de Organismo de Cuenca), incrementamos el contador de celdas en la zona. 2. Acumulamos el valor de la celda por el coseno de la latitud. 3. Acumulamos los cosenos de las latitudes. 4. El valor de la zona va a ser igual a la sumatoria de las Cij por el coseno de las latitudes de ijc entre el acumulado de los cosenos de las latitudesCij . Algunas acotaciones de este algoritmo son las siguientes:  El área de la zona es igual al número de elementos por dos cajitas latitud y longitud por el radio terrestre por  medialatitudcos 360 2 .  El área del polígono es el pseudodeterminante de   terrestreradiomedialatitudcos 2  .  El área se tiene en grados. Se tiene que convertir a kilómetros cuadrados. 48 3.8. Trazado de isoyetas Una vez calculada la superficie de precipitación por cualquiera de los métodos, valuada en una malla regular, entonces se procede a hacer el trazado de las isoyetas. Para ello se utilizó el método llamado marching triangles. Como ya se ha dicho la superficie calculada es ( )I B y la isoyeta es el conjunto de puntos de dicha superficie que tienen un valor específico al que llamaremos Iso , función que se puede definir como     | Iso I B B I B Iso  . Así cuando en un mapa se trazan las isoyetas cada 500 mm entonces se trazan las funciones      500 1000 1500, , ,I B I B I B hasta llegar al límite de escala de los mapas, en el caso particular de los que se presentan en esta tesis el valor máximo de la escala es 5000 mm por lo que en los mapas se trazan las isoyetas de valores quinientos, mil, mil quinientos y así hasta llegar a cinco mil:        500 1000 1500 5000, , ,I B I B I B I B . La implementación de marching triangles es una particularidad en la que se establecieron convenciones de orden y procesamiento; la siguiente descripción es de dicha implementación particular. Se tiene la malla regular de la Fig. 21 donde se considera el origen en la esquina inferior izquierda, el índice i se incrementa hacia la derecha hasta llegar al valor máximo in ; el índice j se incrementa hacia arriba hasta el valor jn . La malla se recorre una vez en este mismo orden. Fig. 21 Triangulación de la malla para el trazado de isoyetas. Sobre la malla de la Fig. 21 se hace una triangulación en donde cada conjunto de 4 puntos se establecen dos triángulos, el inferior (en verde) y el superior (en azul). El triángulo inferior tiene los vértices       1,1,,1,,  jijiji y sus aristas son nombradas 1, 2 y 3 en el sentido contrario a las manecillas del reloj partiendo del punto  ji, . El triángulo superior tiene los vértices       1,,1,1,,  jijiji y también siguiendo la convención de ir en sentido contrario a las manecillas del reloj, sus aristas son nombradas 1, 2 y 3 partiendo del punto  ji, . 49 Si se está buscando el trazo de la isoyeta con valor k, entonces se prueba cuándo k está entre los valores de dos vértices de cada triángulo. Al encontrarlo se ubica mediante interpolación lineal. Si se encuentra que la isoyeta cruza en un punto de una arista de cada triángulo entonces se sabe que para cada triángulo va a haber un punto de entrada y otro de salida de la isoyeta, excepto en los casos límite donde el valor k de la isoyeta sea exactamente igual al valor de al menos un vértice del triángulo. En esta implementación se repasan primero los triángulos inferiores y luego los superiores. Se lleva cuenta de cuáles ya han sido procesados. Y por el orden establecido se sabe que si la isoyeta entra por la arista número 1 entonces sólo puede salir por las aristas dos o tres. La arista dos de un triángulo inferior tiene una adyacencia única, la de un triángulo superior y por la arista 3. La arista tres de in triángulo inferior también tiene una adyacencia única, la de un triángulo superior y por la arista 1. ALGORITMO MARCHING TRIANGLES Algoritmo - 8 ENTRADAS: VALOR DE LA ISOYETA A TRAZAR Y MALLA DE DATOS CON LA SUPERFICIE DE PRECIPITACIÓN CALCULADA. SALIDAS: TRAZOS DE LAS ISOYETAS SOBRE LA IMAGEN EN MEMORIA. 1. Se repasan los triángulos en el orden establecido. 2. Se comparan los valores de cada vértice del triángulo con el valor de la Isoyeta a trazar para encontrar si cruza cada arista. 3. Cuando se encuentra un cruce de la isoyeta buscada en una arista se hace interpolación lineal para ubicar en qué punto de la arista del triángulo se ubica la isoyeta y se marca el triángulo como ya procesado. 4. Cuando se tienen determinados el punto de entrada y el punto de salida de la isoyeta para un triángulo, se traza el segmento en la imagen en memoria. 5. Además se llama recursivamente a esta función para que encuentre por dónde sigue la isoyeta en el triángulo adyacente. 6. Cuando no se encuentra cruce se sigue buscando y se marca cada triángulo donde ya se buscó como ya procesado. 3.9. Cálculo de la constante de orografía En este trabajo se propone poner un factor de aumento o decremento de precipitación en el punto de cálculo respecto al punto de medición en función directa de la diferencia de altitud, sin embargo el fenómeno físico es mucho más complejo y esta es sólo una simplificación. Se usaron ideas conocidas en el ambiente meteorológico e hidrológico, como las relaciones que hay entre las altitudes y las precipitaciones fueron también evaluadas mediante la aplicación del nivel de correlación de Spearsman [8]. 50 En realidad la cosa puede ser tan compleja que para que sea adecuado el uso de este algoritmo se requiere que las constantes de ajuste se determinen para cada una de las zonas de interés o cuencas. Para la determinación de el factor de precipitación por altitud  se analizaron los datos de las estaciones automáticas de altura. Dado que la medición se hace principalmente en las zonas habitadas, casi no existe medición en altura, pero existen dos valiosos casos en los que estamos basando nuestras estimaciones. Se usaron datos de estaciones automáticas los datos de cuatro años cada 10 minutos, convertidas a dato diario de las estaciones automáticas de altura de Cerro Catedral a 3,754 metros sobre el nivel del mar (msnm) y Nevado de Toluca a 4,139 msnm. Se compararon sus series de datos con estaciones del valle como San Juan de las Huertas a 2,282 msnm, tal y como se muestra en la Fig. 22. Con cálculos y criterio experto se determinó que se va a usar una  de 0.4 por cada 500 m de diferencia de altura. Fig. 22 Un año de las series de datos diarios de las estaciones Cerro Catedral, Nevado de Toluca y San Juan de las Huertas. La constante de orografía se determinó 0.4 por cada 500 m de diferencia de altura. Esta adyacencia interrumpe el orden inicial para continuar de manera recursiva el algoritmo siguiendo el trazo de la isoyeta. Este proceso recursivo puede ser interrumpido en cualquier momento y dado que está determinado de manera única el trazo de la isoyeta, los tramos siempre se unirán en la imagen en memoria, sin importar en qué orden fueron encontrados. Esto permite que este algoritmo sea procesable por hilos mediante el método de separación de dominio, siempre que se lleve un control unificado de cuáles triángulos han sido ya procesados. El flujo general del algoritmo se muestra en Algoritmo - 2. Está sujeto a restricciones de nivel máximo de llamadas recursivas como un parámetro, establecido a 1000 en nuestras pruebas, que depende directamente del tamaño de la pila que pueda manejar el lenguaje y la computadora. 51 4. Implementación En este capítulo se describen las partes más importantes de la implementación, preprocesamiento y la parte de paralelización y uso de GPUs para algunos de los cálculos. Los periodos de tiempo representados en los mapas para este trabajo son de un año, a menos que se indique otra cosa. Los mapas de precipitación areal de periodos muy cortos de tiempo no tienen sentido físico debido a la característica intrínseca de localidad de la precipitación. El efecto de las montañas en las precipitaciones areales es evidente cuando se trata de medidas de precipitación de periodos mayores a 6 horas [9]. 4.1. Datos utilizados Para la generación de las superficies de precipitación se utilizaron datos de las estaciones pluviométricas, cartografía digital de México, cartografía digital de las cuencas hidrológicas y el modelo digital de elevación de mayor resolución y más reciente disponible. La precisión de las ubicaciones geográficas, escala, proyecciones y transformación de proyecciones a que fueron sujetas son de una calidad tal que todos empatan a la perfección a las escalas que se manejan en este trabajo dond se usa una resolución máxima 0.001 de grado, aproximadamente. 4.1.1. Precipitación diaria La CONAGUA a través del generoso acuerdo del Gerente de Aguas Superficiales, el M.C. Horacio Rubio Gutiérrez nos extendió permiso de uso de la base de datos de precipitación diaria de Aguas Superficiales de 1950 a 2008. Esta base conforma toda la información disponible en el país para este tipo de datos. Son 1,446 estaciones de medición pluviométrica, cuya ubicación puede verse en la Fig. 23. Esta información fue convertida a una base de datos en MySQL© y los programas acceden a ella en este formato. 52 La información, dado que son datos reales, es cruda y tiene diferentes calidades debidas a las múltiples condiciones de operación que siempre son variables. Por ello se utilizó un criterio de calidad de datos, el descrito por la función de calidad de la ecuación (3.16). Fig. 23 Ubicación de las estaciones pluviométricas con registro diario de la CONAGUA. 4.1.2. Cartografía digital Se utilizó cartografía digital de México y de sus estados basados en cartografía del INEGI. La cartografía de cuencas hidrológicas y organismos de cuenca fue basada en información de la CONAGUA. Esta cartografía fue convertida a polilíneas y después cargada a la base de datos. Los programas la manejan de ahí directamente. La proyección cartográfica en la que se encuentran y manejan en los programas es en coordenadas geográficas. Es notable la calidad y la precisión de todos los conjuntos de datos utilizados, ya que a pesar de provenir de fuentes independientes, cazan perfectamente utilizando sólo las posiciones geográficas y sin necesidad alguna de ajuste de otro tipo. 4.1.3. Modelo digital de elevación En junio de 2009 la NASA puso a disposición a través de su página Web [10] un nuevo modelo digital de elevación sin precedentes en la historia. Con una sorprendente resolución de 30 m y una cobertura mundial, este modelo fue obtenido a través de imágenes de radar satelital ASTER con la colaboración del Ministro de Comercio Económico e Industria de Japón. 53 Fig. 24 Índice de archivos DEM del Modelo Digital de Elevación ASTER Global de la NASA. Mosaicos de 1º × 1º. Escala de color arbitraria. Para México los archivos DEM del tipo geo tiff en mosaicos de un grado por un grado suman 14 GB. Este volumen de información y debido al tipo y densidad de información en el DEM, no es procesable en una computadora personal todo al mismo tiempo en memoria RAM ni puede subirse a una base de datos sin necesidad de herramientas de súper cómputo para hacer búsquedas y consultas. El modelo original puede verse en la Fig. 24. Fig. 25 Malla preprocesada a 0.02 del modelo digital de elevación ASTER Global de la NASA visualizado mediante el programa desarrollado. Escala de color arbitraria. Por ello se generó una estrategia de preprocesamiento generando las mallas fijas en archivos binarios (*.eip) a las resoluciones 0.1, 0.01, 0.02, 0.05, 0.001, 0.002 y 0.005 de grado. Resoluciones intermedias son generadas a través de interpolación lineal basada en estas mallas preprocesadas. La malla de resolución 0.02, aproximadamente 2 km se muestra en la Fig. 25. 54 4.2. Cálculo por malla Para contener el cálculo de las superficies de precipitación y dado que no todos los métodos son funciones explícitas simples sino funciones definidas a pedazos o discontinuas y para fines de facilitación del manejo las Interpolaciones se calcularon sobre malla regular rectangular cuya resolución espacial es definible como un porcentaje de grado, que no es del mismo tamaño para la latitud y la longitud. Una malla de baja resolución para representar el mapa de México podría ser de un décimo de grado, lo que nos da una malla de 321 columnas por 191 renglones, 61,311 puntos y representaría cada celda aproximadamente 10 km por 10 km. Una de alta resolución podría ser de un centésimo de grado 3,201 columnas por 1,901 renglones, 6,085,101 puntos en la malla y representaría un kilómetro cuadrado cada punto. La malla donde cada punto está a aproximadamente 20 km es la de 0.02 de grado. Para los mapas de los experimentos donde se interpola todo el territorio mexicano, esta malla tiene 1600 puntos en el sentido de la longitud, y 950 puntos en el sentido de la latitud, lo cual hace una malla con 1’520,000 puntos. Esta es la resolución que fue seleccionada para hacer los experimentos, de manera que cada uno de los mapas mostrados tiene ese millón y medio de puntos. Sin embargo no todos los puntos son valuados porque se hizo la optimización de tener una máscara de validación de punto interior al área de interés. El cálculo efectivo se hace en alrededor de un millón de puntos. Algunos pocos puntos calculados son descartados al hacer el recorte final del perfil del país. En la siguiente figura se muestra una malla calculada y recortada con el perfil del país. Fig. 26 Cálculos resultantes en una malla de 20 km después de ser perfilado y la imagen es el estado previo a la rotulación. Son aproximadamente un millón de puntos calculados aunque la malla contiene en realidad millón y medio 1600 × 950. 55 Independientemente del método seleccionado, los valores de la superficie de precipitación interpolada por cada uno de los 16 métodos mostrados, tanto los propuestos como los comparativos, para fines de cálculo son valuados en una malla regular rectangular fija a una resolución. Esto para facilitar la comparación de los métodos: en cada mapa, independientemente del método, es valuado en los mismos puntos en la malla. Aunque la superficie de precipitación que entrega cada método sea continua y o tenga escalones y puntos de discontinuidad, para fines de cálculo todas son muestreadas en esta malla. Se seleccionó la malla de resolución 0.02 porque las áreas definidas con la suma de los rectángulos definidos con cada punto de la malla y sus contiguos dan una diferencia promedio de dos centésimas, al hacer la comparación con el área que se obtiene directamente al calcularlas con los polígonos definidos que constituyen el mapa base de ellas. Dos centésimas es pequeño comparado con las 81 centésimas de diferencia promedio que se obtienen con la malla a 0.05 km, distancia entre puntos en la malla separados aproximadamente 50 km. Esto se muestra en la siguiente tabla, donde además se incluyó la cantidad de puntos en la malla que están conformando a dicha región. En la tabla siguiente se observa una pequeña muestra de los cálculos para las 37 Regiones Hidrológicas y sus 159 subregiones. Tabla 2 Muestra del cálculo de las áreas con una malla de 0.02 y con una de 0.05. La malla de 0.02 da una diferencia de dos centésimas de punto porcentual mientras que la de 0.05 da ocho centésimas. Con malla 0.02 Con malla 0.05 Zona Área calculada por el polígono. Númer o de puntos. Área calculada por las celdas. Porcentaje de la diferencia entre áreas para 0.02. Área calculada por el polígono. Número de puntos. Área calculada por las celdas. Porcentaje de la diferencia entre áreas para 0.05. RH01 26757.83 6210 26460.48 1.11% 26757.83 1001 26657.55 0.37% RH01A 9274.2 2078 8935.58 3.65% 9274.2 335 9003.31 2.92% RH01B 9741.87 2306 9767.87 -0.27% 9741.87 371 9821.87 -0.82% RH01C 7651.45 1826 7663.31 -0.16% 7651.45 295 7737.8 -1.13% ... ... ... ... ... ... ... ... ... RH37C 7714.6 1710 7713.85 0.01% 7714.6 272 7668.74 0.59% RH37D 9265.72 2056 9290.79 -0.27% 9265.72 327 9235.44 0.33% RH37E 14662.15 3226 14672.9 -0.07% 14662.15 516 14668.35 -0.04% RH37F 12996.71 2836 12961.04 0.27% 12996.71 464 13253.54 -1.98% RH37G 11359.18 2487 11374.13 -0.13% 11359.18 392 11204.92 1.36% RH37H 11505.94 2521 11505.54 0.00% 11505.94 405 11552.32 -0.40% TOTALES 3937064.1 9 3936295.0 1 0.02% 3937064.1 9 3933837.5 8 0.08% Además para comprobación de la precisión de lo representado en las mallas, con dichos valores se hizo el cálculo de la superficie de México como se muestra en la siguiente tabla, donde se coteja de la precisión de los mapas utilizados y la adecuada selección de la resolución espacial para la malla. 56 Tabla 3 Comparación de las áreas del territorio mexicano calculado con nuestras mallas de 0.02 y 0.05 contra los datos oficiales. La diferencia es muy pequeña. Área total de México incluyendo el territorio insular Diferencia INEGI 1’964,375 km2 Con malla de 0.02 de grado. 1’968,532 km2 0.21% Con malla de 0.05 de grado. 1’968,147 km2 0.19% 4.3. Programas Se programaron rutinas en Delphi y en C++ para este proyecto. Además se utilizaron herramientas externas para la conversión de imágenes, el cálculo de los diagramas de Voronoi e implementaciones de LAPACK© con ACML© y CUDA©. 4.3.1. Programa de generación de mallas de preprocesamiento de orografía Lectura de los archivos de orografía de México que mide 14 GB y que contienen información grado por grado y la extracción y escritura de archivo (*.eip) que contiene exclusivamente los datos requeridos para un tamaño de malla prefijado. Se generaron mallas de preprocesamiento de orografía para las resoluciones 0.1, 0.05, 0.02, 0.01, 0.005, 0.002 y 0.001 de grado. Dada la naturaleza del formato en que se encuentra esta información, formato DEM, el volumen de datos y lo detallada de la información, malla de 30 m, se hubiera requerido súper cómputo para hacer uso de esta información si no hubiéramos seleccionado esta estrategia. 4.3.2. Rutina de consulta a la base de datos climatológica. La base de datos original proporcionada por la CONAGUA tiene varios millones de registros y mide más de 40 GB. Las consultas dinámicas que se pueden hacer permiten el manejo de cualquier subconjunto de datos proporcionando una entrada adecuada para cualquier experimento que se quiera realizar. 4.3.3. Rutinas de los 16 métodos de cálculo de la superficie de precipitación con salidas parciales en imagen, archivos de datos Se programaron los métodos que son estándar en el ambiente de las isoyetas y contra los que se contrastó en esta investigación: Thiessen, Kriging ordinario con variogramas de tres tipos Gausiano, exponencial y esférico; y la muy utilizada versión del Spline de curvatura continua sujeta a tensión. 57 Parámetros de comando de línea: -d SIH -o mexico -rg -118/-86/14/33.5 -g NAETISOL -fi 1971/01/01 -ff 1980/12/31 -fp AÑO -I 0/5000/500 -D 500 -dm SI -H 4 -im 0.02 -Mp PerfMexCG:MEXICO -Mb MEXRHCG -t 6 -a 0.4 -Ed SI -GM NO -Mo mexico -iA 500 -iC 250 -GPU AMD -Ej 5/40/192/3 -E1 "Precipitación Media: #media mm" -E2 "Estaciones en cálculo: #estac" -E1p 1750/100/48 -Le 20 Donde el significado de cada uno de estos parámetros se explica brevemente en la Tabla 4. Tabla 4 Parámetros de comando de línea. Parámetro Significado -d SIH Conexión a la base de datos. -o mexico Opción de mapa. -rg -118/-86/14/33.5 Paralelos y meridianos que delimitan el área geográfica. -g NAETISOL Grupo de estaciones a utilizar. -fi 1971/01/01 -ff 1980/12/31 Fechas inicial y final. -fp AÑO Acumulación anual. Los datos originales son diarios. -I 0/5000/500 Trazado de isoyetas, de cero a 5 mil cada 500. -D 500 -dm SI Escritura de archivos de mallas. -H 4 Número de hilos a utilizar en los procesos paralelos. -im 0.02 Resolución de la malla. -Mp PerfMexCG:MEXICO Perfil del territorio. -Mb MEXRHCG Perfiles de los subterritorios. -t 6 Tipo de cálculo -a 0.4 Factor de altitud. -Ed SI Escribe archivos para exportar datos de cálculo. -GM NO Generar mapa con técnica GMT. -Mo mexico Nombre del mapa con la orografía. -iA 500 Nomenclatura de las isoyetas sobre valores cada 500. -iC 250 Valor de isoyeta que no va nombrada. -GPU AMD Tipo de procesamiento paralelo en GPU. -Ej 5/40/192/3 Parámetros de las referencias finales para presentación estilo mapa. Marcas de paralelos y meridianos cada 5 grados, tamaño de los números de referencia, grueso de la línea y tamaño del círculo de las estaciones. -E1 "Precipitación Media: #media mm" Título. -E2 "Estaciones en cálculo: #estac" Subtítulo. -E1p 1750/100/48 Posición de los títulos y tamaño de letra. -Le 20 Tamaño de la letra para la escala de color. Ahora, el parámetro más importante es –t, el que determina el tipo de cálculo. Las 16 opciones de tipo de cálculo se muestran en la Tabla 5, donde además se describen los diferentes nombres con que pueden ser identificados dentro de los menús del programa en su interfaz gráfica, en la identificación breve de los mapas y en este documento. 58 Tabla 5 Parámetro de tipo de cálculo para el programa Etisol. -t Tipo de cálculo Nombre en menús. Nombre en mapas. Denominación en tesis. 0 Thiessen Thiessen. Thiessen 1 4EstInvDist Metodo: 4 cercanos / dist. 2 4EstInvDist2 4 cercanos / (dist*dist). 3 16EstInvDist 16 cercanos / dist. 4 16EstInvDist2 Metodo: 16 cercanos / (dist*dist). 5 Voronoi Vecinos Voronoi. 6 VoronoiDifAlt Voronoi con Orografia. Vecinos Voronoi con orografía. 7 ThiessenDifAlt Thiessen con Orografia. 8 4EstDifAlt 4 cercanos con Orografia. 9 16EstDifAlt 16 cercanos con Orografia. 10 KrigGauss Kriging Gausiano. 11 KrigExp Kriging Exponencial. 12 Laplaciano Laplaciano. 13 Biarmónico Biarmonico. 14 TensorLaplBiar m Tensor Lapl-Biarmonico. 15 KrigEsferico Kriging Esferico. Tabla 6 Archivos que se general al correr Etisol. Extensión de los archivos generados .asc Malla calculada que puede abrirse en Global Maper® o ArcGIS©. .bat Archivo auxiliar para convertir el bmp a png .etv La rutina externa de convex hull produce este archivo y contiene los puntos y los polígonos de Voronoi producidos. .mll nombre corto Archivo donde se ve la malla tal cual. Se hizo para depuración de isolíneas. .mll nombre largo Archivo donde se ve la malla tal cual. Contiene los mismos valores que la otra malla. Se hizo para depuración de isolíneas. .png *.png Malla escala *I.png Resultado final Isolíneas solitas *Mll.png Malla *MllIs.png Resultado final Malla isolíneas *MllM.png Resultado final Malla mapa *Vor.png Diagrama de Voronoi *VorM.png Diagrama de Voronoi con el perfil encima Orografía. Sale en interactivo o descomentando y recompilando. .txt *Tiempos.txt Rastreo de tiempos de proceso .trg Máscara de validez .xls Resultados Precipitaciones areales por regiones hidrológicas. Precipitación areal nacional. Datos originales. Estación, total, coordenadas, fechas inicial y final. Comparación: Estación, valor interpolado, valor original, coordenadas .xy Malla de depuración .xyz Malla 3D punto a punto 59 Finalmente en la Tabla 6 se enlistan los tipos de archivos producidos por el programa. Se da el tipo de extensión, la forma del nombre para algunos tipos, el tipo de información que contiene y si es o no un archivo que contiene un resultado. Algunos archivos son para rastreo, depuración y muestra de los pasos del cálculo. 4.3.4. Procesamiento paralelo El programa tiene unas partes que corren en paralelo. La inversión de matrices que para los conjuntos de datos aquí manejados es de tamaño típico mil por mil, pero es variable en función de las estaciones que hayan pasado el control de calidad para el cálculo, se hace mediante GPUs con ACML© y con CUDA© a través de funciones de LAPACK© y exclusivamente para el método Kriging. Además hay otras partes que corren en procesos paralelos utilizando hilos, como el trazado de las isoyetas y el método Kriging, para lo cual se puso la opción para el uso de múltiples CPUs con ACML© . Este procesamiento en paralelo es controlable desde los parámetros del programa. El programa puede correr perfectamente en computadoras sin GPUs programables y con un sólo procesador, pero también se pueden aprovechar los mejores recursos que se tengan disponibles. 4.3.5. Lista de rutinas y otros programas. Además del programa principal fueron creados otros programas menores y rutinas que realizan tareas diversas, aquí se enlistan las consideradas más importantes.  Rutina de trazado de isoyetas: Implementación recursiva altamente paralelizable mediante la técnica de partición de domino.  Rutina para verificación cruzada: Generación de estadísticas para subconjuntos determinados por un grupo de estaciones para hacer la verificación cruzada y medir la calidad de la predicción.  Rutinas de generación del mapa: Los resultados se escriben en una imagen tipo PNG. Son las rutinas de etiquetado, escala, graduación espacial de los resultados.  Programa de cálculo de mapas de error: La verificación con medición de distancia usando datos artificiales se hizo en mapas con desviación positiva y negativa usando como estándar del cálculo de precipitación areal por Thiessen. 60 4.4. Radio máximo de competencia Para todos y cada uno de los métodos propuestos se incorporó una restricción de radio máximo de competencia del valor de una estación. Se hicieron pruebas para que los datos disponibles fueran suficientes para cubrir el territorio nacional sin dejar zonas sin dato. Se seleccionó 200 km como radio máximo de competencia. En las siguientes figuras se muestran resultados con diferentes radios. Fig. 27 Pruebas de cobertura con diferentes radios de competencia de las estaciones. Tal y como lo establece De la Cruz en [11], las funciones de interpolación aquí propuestas son variantes de funciones de base radial. 4.5. Procesamiento en paralelo Uno de los objetivos de este trabajo fue producir una herramienta que en minutos permitiera el cálculo de las superficies de precipitación, sus isoyetas y el cálculo de la precipitación areal en minutos en computadoras personales. 61 Al principio se veían muchos retos para procesamiento en paralelo, pero al hacer la implementación se fueron resolviendo de otra manera, especialmente el preprocesamiento de las mallas de orografía, el preprocesamiento de la interioridad de las estaciones a los polígonos de las Regiones Hidrológicas, Estados y Organismos de Cuenca, resultados que se almacenan en la base de datos. Los dos únicos retos de procesamiento en paralelo que se conservó fue el cálculo completo de Kriging, no una aproximación numérica, y el trazado de isoyetas por un método recursivo y paralelizable. El cálculo de Kriging implica la inversión de la matriz de varianzas y covarianzas. Así que para un conjunto de datos de 900 estaciones, entonces se tiene una matriz de varianzas y covarianzas de tamaño 900 × 900 que hay que invertir. Un reto interesante para una computadora personal. Esto se logró a través de las funciones de álgebra lineal de LAPACK©, biblioteca de procesamiento de álgebra lineal. Se utilizó sorprendentemente rápido y transparente ACML© de AMD©. ACML© permitió la implementación para múltiples CPUs en una misma computadora y la implementación para GPUs AMD©. Estas herramientas son completamente gratuitas. También se hizo la implementación para uso de GPUs con CUDA© de nVIDIA©, mediante la biblioteca que implementa LAPACK© llamada CULA©, que se compró una licencia de un año para una computadora expresamente para este proyecto. Nuestro programa tiene parámetros de entrada que permiten indicar si se quiere utilizar o no procesamiento en paralelo e indicar el número de hilos para múltiples CPUs que se quieren utilizar para cada corrida. 63 5. Resultados En este capítulo se presentan los resultados de los experimentos hechos para esta tesis y se valida la calidad de las interpolaciones por diferentes métodos. Se generaron 43 GB de archivos de resultados que incluyen mapas, mallas de procesamiento y rastreo, archivos de tiempos, tablas de datos y resúmenes. De todo esto se presenta una selección, primordialmente resultados de datos anuales de la década 1971 a 1980 con una malla de resolución de 2 km a menos que se indique otra cosa. 5.1. Resultados de los experimentos realizados 5.1.1. Experimentos realizados Se seleccionó la década 1971-1980 para generar los mapas de comparación de cada uno, con los datos de esta década se presentan mapas de precipitación areal anual para México. Este conjunto de datos nos da una adecuada visión del comportamiento de los métodos. La variación del resultado por el conjunto de datos se muestra en la Fig. 57, donde se muestra la diferencia del cálculo de precipitación areal tomando datos de cuatro décadas diferentes, se aprecia que los resultados son similares a nivel nacional. Ahora, cuando se toman zonas pequeñas, por ejemplo las 37 Regiones Hidrológicas, las variaciones también son pequeñas pero mayores que a nivel nacional, tal y como lo muestra la gráfica de diferencias de precipitación areal por Regiones Hidrológicas con datos de 4 décadas diferentes de la Fig. 58. Se generaron mapas de México con cartografía y datos reales. El perfil de la República Mexicana, los mapas por Regiones Hidrológicas, Estados y Organismos de Cuenca fueron convertidos a proyección rectangular de igual área para hacer los cálculos e imágenes. Éstos empataron perfectamente con los datos de orografía de la NASA y con las coordenadas de las estaciones pluviométricas. Los datos utilizados fueron datos diarios de precipitación de la CONAGUA, base de datos de Aguas Superficiales de 1950 a 2008. Los datos utilizados, dado que son reales están sujetos a las complejidades y los problemas propios de cualquier conjunto de mediciones. Por ello, para poder ser tomados en cuenta, conjunto definido por la ecuación (3.17) en estos mapas, fueron sujetos a la función de calidad de datos descrita en (3.16). 64 Se puede considerar que las mediciones tienen en sí mismas un error de cuantificación respecto al valor real de la función de precipitación. Los modelos Kriging y Spline de Curvatura Continua Sujeto a Tensión que consideran cierta variabilidad de las mediciones lo hacen suponiendo un modelo lineal, el decir, suponen que la función real es la función medida más una función de error. Esta superposición lineal de funciones no cambia en sí misma ningún comportamiento de los métodos de interpolación, sino más bien es una manera de interpretar los datos y la función obtenida. Dado que el objetivo principal de esta tesis es la evaluación de los métodos de interpolación para precipitación y no el desempeño y calidad de los instrumentos de medición, en esta tesis se consideran las mediciones como mediciones directas de la función de precipitación real. Se generaron experimentos a resoluciones de 0.1, 0.05, 0.02, 0.01 de grado. Las diferencias de precisión debidas al cambio de resolución en el área geográfica de cálculo a nivel nacional pueden verse en la Fig. 43 y a nivel local, por Regiones Hidrológicas en la Fig. 45. Las diferencias en la precipitación areal a nivel nacional puede verse en la Fig. 42 y por Regiones Hidrológicas en la Fig. 44. Por los resultados obtenidos se considera la malla de 2 km una excelente resolución para mapas a nivel nacional. Se generaron 16 tipos de mapas con los algoritmos de comparación y los nuevos. Tabla 7 Tabla de métodos de interpolación propuestos y contra los que se compara. Métodos propuestos Métodos contra los que se compara 1. Vecinos Voronoi. 2. Vecinos Voronoi con Orografía. 3. Cuatro Estaciones Cercanas por el inverso de la distancia. 4. Dieciséis Estaciones Cercanas por el inverso de la distancia. 5. Cuatro Estaciones Cercanas por el inverso de la distancia al cuadrado. 6. Dieciséis Estaciones Cercanas por el inverso de la distancia al cuadrado. 7. Cuatro Estaciones Cercanas con Orografía. 8. Dieciséis Estaciones Cercanas con Orografía. 9. Thiessen con Orografía. 1. Kriging Gaussiano. 2. Kriging Exponencial. 3. Kriging Esférico. 4. Spline de Curvatura Continua sujeto a Tensión con T=1, Laplaciano. 5. Spline de Curvatura Continua sujeto a Tensión con T=0, Biarmónico. 6. Spline de Curvatura Continua sujeto a Tensión con T=0.75, Tensor Biarmónico-Laplaciano. 7. Se calcula el error contra Thiessen. Las imágenes producidas son archivos individuales en formato png de unos 850 KB a una resolución de 3,500 píxeles de ancho por 2078 píxeles de alto en color de 24 bits. De acuerdo a la Tabla 9 se han numerado los siguientes mapas resultantes de los métodos de interpolación propuestos. Todos los métodos aquí presentados respetan el valor de las mediciones en los puntos medidos. En todos los mapas que se muestran a continuación se muestra la superficie de precipitación interpolada por algún método. La precipitación representada es el valor anual de precipitación nacional tomando los datos de la década 1971 a 1980. La escala numérica y en representación de color es la misma para cada mapa, de manera que se facilita hacer comparaciones visuales directas. El color blanco es representa el valor de 0 mm, y el color morado representa 5,000 mm. Se usaron exactamente el mismo conjunto de datos los mismos criterios de selección 65 y calidad de datos para todos los mapas presentados. En el rótulo que está encima de la escala del lado inferior izquierdo de los mapas está una identificación abreviada de la información presentada. La resolución a la que están hechos estos mapas es de 0.02 de grado, cada punto en la malla está separado aproximadamente por dos kilómetros. Esto es una resolución alta para un mapa a nivel nacional. Son más de un millón de puntos calculados para cada mapa. La malla es de 1600 puntos en el sentido de la longitud, y 950 puntos en el sentido de la latitud. El color con que está rellena la superficie del mapa corresponde exactamente al color de la escala de cada punto en la malla, con la técnica llamada generar el píxel, donde se rellena del color correspondiente al área entre píxeles asignándole la de uno en específico. La imagen escrita a disco tiene una resolución de 3,500 x 2,078 píxeles y está en formato png. En la Fig. 28 se muestra el resultado de la interpolación de la precipitación anual por el método de Vecinos Voronoi. La superficie de precipitación que produce este método es discontinua. Cada polígono de Voronoi tiene asignado un mismo valor el cual fue calculado con base a la ponderación de la 1-Vecindad de Voronoi. Este método nunca muestra valores más allá de los extremos medidos. Es un método adaptable a la cantidad de mediciones disponible; donde hay más mediciones, se utilizan para obtener el valor ponderado sobre ellas, cuando hay poca medición se utiliza la disponible, siempre hasta llegar al límite del radio de competencia máxima que es de 200km para estos mapas. El resultado de este método puede pensarse como una función discontinua conformada por un conjunto de escalones a diferentes niveles donde cada escalón tiene la forma de un polígono de Voronoi. Fig. 28 Interpolación por Vecinos Voronoi. En la Fig. 29 se presenta la interpolación de la precipitación anual nacional interpolada por el método de Vecinos Voronoi con Orografía. Esta interpolación tiene características similares a Vecinos Voronoi, dado que es una modificación directa de dicho método, sólo que además toma en cuenta cómo varía la precipitación en función de la altitud con el parámetro calculado previamente, 0.4 por cada 500 m de diferencia de altitud entre el punto en la malla y el medido. 66 Esta también es una superficie discontinua no diferenciable entre los bordes de los polígonos de Voronoi, pero es mucho más suave que el método sin orografía. Aquí los polígonos de Voronoi no son superficies planas paralelas al dominio sino cada polígono tiene la forma de la orografía. Sorprendentemente esto sigue de una manera más fiel la vegetación, ya que en las cañadas típicamente hay más humedad. Fig. 29 Interpolación por Vecinos Voronoi con Orografía. En la Fig. 30 se muestra la superficie de precipitación anual interpolada por el método de 4 Estaciones Cercanas por el inverso de la distancia. La superficie producida Es muy suave en casi todos los puntos pero puede presentar discontinuidades. Los valores interpolados se parecen al método de Thiessen pero están más suavizados, ya que están tomadas en cuenta las cuatro estaciones más cercanas al punto en la malla que se está calculando. En la Fig. 31 se muestra la interpolación de precipitación anual por el métodos de 16 Estaciones Cercanas por el inverso de la distancia. Sus características son similares al anterior, diferenciándose en que la superficie producida es aún más suave por el efecto de tomar más estaciones en cuenta. Es de notarse la localidad y adaptabilidad de estos métodos. En las Fig. 32 y Fig. 33 se presentan las superficies de precipitación anual interpoladas por los métodos 4 y 16 Estaciones Cercanas por el inverso de la distancia al cuadrado. Esta variación hace que el decaimiento debido a la distancia de la medición sea más lento por lo que produce superficies más suaves, aunque aún pueden presentar discontinuidades. Estas discontinuidades no constituyen ningún problema, dado que se busca tener mejor precisión en la interpolación y la continuidad no tiene sentido práctico para este tipo de mapas. Alta resolución y mayor precisión son características que ayudan a la lectura directa del valor a partir de ellos. 67 Fig. 30 Interpolación por 4 Estaciones Cercanas por el inverso de la distancia. Fig. 31 Interpolación por 16 Estaciones Cercanas por el inverso de la distancia. 68 Fig. 32 Interpolación por 4 Estaciones Cercanas por el inverso de la distancia al cuadrado. Fig. 33 Interpolación por 16 Estaciones Cercanas por el inverso de la distancia al cuadrado. En las Fig. 34 y Fig. 35 se muestran las interpolaciones de precipitación anual con los métodos de 4 y 16 Estaciones Cercanas con Orografía, respectivamente. En este método la ponderación se hace por el inverso de la distancia además del factor por la diferencia de altitud entre el punto en la malla y el punto medido a considerar. Estas superficies sí pueden producir valores mayores al máximo medido y menores al mínimo medido, cosa que es clara en el punto morado de la Fig. 34 y en la zona azul de la costa de Chiapas en la Fig. 35, pero estos valores que 69 superan a los extremos no están demasiado desviados de los medidos. Jamás producen precipitaciones negativas por las restricciones del algoritmo. Fig. 34 Interpolación por 4 Estaciones Cercanas con Orografía. Fig. 35 Interpolación por 16 Estaciones Cercanas con Orografía. En la Fig. 36 se presenta el mapa de la precipitación anual interpolado por el método de Thiessen con Orografía. Este método equivaldría al de Vecinos Voronoi con Orografía salvo que en lugar de 1-vecindad se emplearía 0-vecindad, es decir, sólo el valor del polígono de Voronoi donde cae el punto en la malla. Por ello la discontinuidad es mayor en la superficie producida. 70 En la Fig. 36 pueden verse con claridad los límites de los polígonos de Voronoi y es muy interesante cómo quedan afectados por el efecto de la orografía. Fig. 36 Interpolación por Thiessen con Orografía. Para cualquiera de los métodos con orografía, las diferencias respecto a los mapas producidos por el método equivalente sin orografía son sólo patentes en las zonas montañosas. En las zonas planas como la Península de Yucatán, los resultados son casi iguales. Esta adaptabilidad en función de la altitud de las montañas es un resultado importante para la precisión de la interpolación. A continuación se muestran ejemplos de otros mapas que se generaron para cada uno de los experimentos. El mapa de isoyetas se muestra en la Fig. 37. El mapa de la interpolación con isoyetas se muestra en la Fig. 38. En la Fig. 39 se hace un acercamiento para ver el detalle del trazo de las isoyetas. El trazo que se muestra en todos los casos es el obtenido directamente por el algoritmo de marching triangles y no tiene ningún procedimiento posterior de regularización ni simplificación. Para algunos casos es más detalle del que se acostumbra a manejar en este tipo de mapas, pero es valioso si se quieren hacer lecturas a nivel local. 71 Fig. 37 Mapa de isoyetas de Vecinos Voronoi con Orografía a cada 500 mm. Fig. 38 Mapa de interpolación por Vecinos Voronoi con Orografía con isoyetas a cada 500 mm. 72 Fig. 39 Trazo de las isoyetas por el método de marching triangles. En todos los mapas están indicadas las estaciones en su ubicación original como se muestra en el acercamiento de la Fig. 40 Fig. 40 Acercamiento que muestra la ubicación de las estaciones en los mapas. 73 Para cada mapa también se produjo un archivo en formato ASCII compatible con programas genéricos de visualización de datos cartográficos como el Global Mapper® y el ArcGIS©. Una visualización en Global Mapper® puede vers en la Fig. 41. Fig. 41 Interpolación por Vecinos Voronoi con Orografía mostrada en Global Mapper®. Para cada experimento se generó también un archivo de rastreo de tiempos de proceso. En Listado - 1 se muestra el del proceso de la Interpolación por Vecinos Voronoi con Orografía. -------------------------- Rastreo de tiempos de proceso con las opciones: Cálculos áreales: 6 - Vecinos Voronoi Dif. alturas Tamaño de malla: 0.02 Con GPU AMD / ACML . Con 4 para proceso en paralelo Carga de datos: 5906 ms Preparación de datos para visualización: 219 ms Carga de alturas en ZC: 578 ms Escritura de archivo xyz para visualización datxyz.txt: 0 ms Inicialización de arreglos MVv y MVz: 32 ms Determinación de altura en malla MVz: 2156 ms Inicialización de arreglos artr, artt para triángulos: 140 ms Cálculo de polígonos Voronoi: 4954 ms Cálculos de malla por hilos: 34093 ms Escritura de malla de datos MVv: 3766 ms Tiempo de proceso en paralelo: 45719 ms Determinación de triángulos artr: 94 ms Escritura de malla arpv: 115625 ms Escritura de malla arpv1 tipo .asc: 5093 ms Contrucción de gráfica 3D (pego): 157 ms Escritura de triángulos artr: 6828 ms Preparación de datos. Cálculo de malla: 174172 ms Cálcula de Precipitaciones areales por Grupos en Mapa: 16750 ms AjustaArx 16750 ms Tiempo de proceso ------------------------------: 196828 ms Exportado de datos: 672 ms Dibujo de malla: 12781 ms Trazado de isolíneas: 5032 ms Listado - 1 74 Ahora se presentarán unos análisis de la precipitación areal calculada bajo diferentes condiciones con la finalidad de ver la variabilidad por los factores que se irán mencionando. En la Fig. 42 se muestra cómo varía la precipitación areal nacional cuando se calcula por el métodos de Vecinos Voronoi con Orografía y con mallas de diferentes resoluciones: malla de 1 km, 2 km, 5 km y 10 km. Para las mallas probadas no hay diferencia significativa entre ellos. Fig. 42 Comparativo de la precipitación areal nacional calculada con mallas a 1 km, 2 km, 5 km y 10 km. No hay diferencia significativa entre ellas. En la Fig. 43 se muestra cómo varían las áreas totales calculadas por malla a diferentes resoluciones: 1 km, 2 km, 5 km y 10 km, y comparadas con la cartografía digital original definida en polígonos a los cuales se les ha calculado el área. Aunque sí hay pequeñas diferencias, no hay cambios significativos. En la Fig. 44 se muestra cómo varía la precipitación areal a nivel de Regiones Hidrológicas, que son 37 en el país, cuando se varía la resolución de la malla a la que se está calculando. A diferencia de las comparaciones a nivel nacional, sí hay variaciones más grandes para algunas Regiones Hidrológicas, sobre todo a las de menor área. Sin embargo son todavía pequeñas. Tal y como puede verse en la Fig. 45 donde se muestran las variaciones de las áreas definidas por las mallas a diferentes resoluciones y comparadas con las áreas originales de las Regiones Hidrológicas definidas por los polígonos de la cartografía original, la variación es todavía pequeña. 75 Fig. 43 Comparativo de las áreas totales conformadas por los elementos de las celdas de la malla a resoluciones de 1 km, 2 km, 5 km y 10 km comparadas con las áreas de los polígonos originales. Fig. 44 Comparativo del cálculo de la precipitación areal por Regiones Hidrológicas calculada por mallas a diferentes resoluciones: 1 km, 2 km, 5 km y 10 km. 76 Fig. 45 Comparativo de las áreas de las Regiones Hidrológicas calculadas por cada elemento de la malla a diferentes resoluciones: 1 km, 2 km, 5 km y 10 km y comparada con las áreas de los polígonos originales. Se generaron también mapas por épocas del año haciendo la comparación de época de lluvias de mayo a noviembre con la época de secas de diciembre a abril. 77 Fig. 46 Mapas de las épocas de secas y de lluvias de cuatro décadas diferentes. Las diferencias se aprecian principalmente en las épocas de lluvias. En la Fig. 46 se muestran de manera comparativa los promedios anuales de las épocas de lluvias de las décadas 1951 a 1960, 1961 a 1970, 1971 a 1980 y 1981 a 1990. Las diferencias se notan en la cantidad de precipitación, pero también en la cobertura de las estaciones, en el mapas de 1951 a 1960 el radio de 200 km no alcanza a cubrir todo el territorio en la zona norte del país. Ahora, respecto a los tiempos de proceso y aceleración por paralelización con múltiples CPUs, con GPU y con CPUs más GPU se tiene lo siguiente. Con base en los archivos testigo donde se rastrean los tiempos de cada paso del proceso y que se generaron para cada experimento, como el que se muestra en Listado - 1 se obtuvieron las gráficas siguientes. En Fig. 47 se muestran los tiempos de proceso con un sólo CPU para los métodos propuestos y Kriging Esférico. En la gráfica están diferenciados los tiempos de la lectura y preparación de los datos, que como se ve es idéntico para todos los procesos excepto para Kriging. También en rojo están los tiempos de cálculo especiales para cada método, aplicable tan sólo para el cálculo del diagrama de Voronoi en los métodos de Interpolación Voronoi y Voronoi con Orografía, y para los métodos Spline de curvatura continua sujeto a tensión como se aprecia en las Fig. 50 y Fig. 51. En amarillo están los tiempos asociados a los cálculos para la valuación de la función en los puntos de la malla y los cálculos por zonas, Regiones Hidrológicas, Estados u Organismos de Cuenca, según sea el caso. Este es mayormente el tiempo de procesamiento de cada método. Puede verse cómo son prácticamente iguales los tiempos de procesamiento para todos los métodos excepto para Kriging. Finalmente en azul está el tiempo de la escritura de resultados, que incluye la escritura de archivos de texto con las mallas de datos, el archivo de rastreo y la imagen del mapa producido. 78 Fig. 47 Gráfica de tiempos separada por procesos para cada uno de los métodos de interpolación cuando se procesa en serie a un sólo CPU. No fueron ejecutados para los experimentos en serie todos los métodos pesados en cálculo por las horas de procesamiento que implicaba. Pero todos los métodos fueron ejecutados en las versiones en paralelo. En la Fig. 48 se aprecia la comparación de los tiempos de procesamiento entre 1 CPU y 4 CPUs para cada uno de los diferentes métodos propuestos. El tiempo de proceso cae a la tercera parte. Fig. 48 Gráfica que muestra la aceleración producida al procesar con 4 CPUs comparada con 1 CPU para todos los métodos de cálculo propuestos. 79 Fig. 49 Gráfica de aceleración del procesamiento del cálculo de los algoritmos que incluye un Kriging. Fig. 50 Gráfica de tiempos de proceso diferenciada por etapas. Se omitieron los Kriging para tener una visión clara de comparación. 80 Fig. 51 Similar a la figura anterior se muestra la Gráfica de Tiempos de proceso diferenciado por etapas que incluye los métodos Kriging. Fig. 52 Precipitación areal anual nacional con datos en la década 1971 a 1980 calculada por los diferentes métodos. 81 -500 0 500 1000 1500 2000 2500 3000 3500 RH01 RH02 RH03 RH04 RH05 RH06 RH07 RH08 RH09 RH10 RH11 RH12 RH13 RH14 RH15 RH16 RH17 RH18 RH19 RH20 RH21 RH22 RH23 RH24 RH25 RH26 RH27 RH28 RH29 RH30 RH31 RH32 RH33 RH34 RH35 RH36 RH37 M ed ia m ap a Región Hidrológica P re ci p . M ed ia A re al Thiessen Thies-Orog Voronoi Voron-Orog 4 Cerc/d 4 Cerc/d2 4 Cerc-Orog 16Cerc/d 16Cerc/d2 16Cerc-Orog Krig-Gaus Krig-Expn Krig-Esf Laplaciano Biarmónico Lapl-Biarm Fig. 53 Precipitación areal anual por Regiones Hidrológicas comparada entre los diferentes métodos de cálculo. - 500.0 1,000.0 1,500.0 2,000.0 2,500.0 OCAVM OCB OCCCN OCFS OCGC OCGN OCLSP OCNW OCPBC OCPN OCPS OCPY OCRB Método de Cálculo P re ci p it ac ió n m ed ia Thiessen Thies-Orog Voronoi Voron-Orog 4 Cerc/d 4 Cerc/d2 4 Cerc-Orog 16Cerc/d 16 Cerc/d2 16Cerc-Orog Krig-Esf Lapl-Biarm Fig. 54 Precipitación areal anual por Organismos de Cuenca comparada entre los diferentes métodos. 82 - 500.0 1,000.0 1,500.0 2,000.0 2,500.0 AGS BC BCS CAMP CHIH CHIS COAH COL DF DGO GRO GTO HGO JAL MEX MICH Método de Cálculo P re ci p it ac ió n m ed ia a n u al ( m m ) Thiessen Thies-Orog Voronoi Voron-Orog 4 Cerc/d 4 Cerc/d2 4 Cerc-Orog 16Cerc/d 16 Cerc/d2 16Cerc-Orog Krig-Esf Lapl-Biarm Fig. 55 Precipitación areal anual por Estados parte 1 comparada entre los diferentes métodos. - 500.0 1,000.0 1,500.0 2,000.0 2,500.0 3,000.0 MOR NAY NL OAX PUE QRO QROO SIN SLP SON TAB TAMS TLAX VER YUC ZAC Método de Cálculo P re ci p it ac ió n m ed ia a n u al ( m m ) Thiessen Thies-Orog Voronoi Voron-Orog 4 Cerc/d 4 Cerc/d2 4 Cerc-Orog 16Cerc/d 16 Cerc/d2 16Cerc-Orog Krig-Esf Lapl-Biarm Fig. 56 Precipitación areal anual por Estados parte 2 comparada entre los diferentes métodos. 83 680.0 700.0 720.0 740.0 760.0 780.0 800.0 820.0 1971-1980 1981-1990 1991-2000 1998-2007 Decada P re ci p it ac ió n m ed ia a re al ( m m ) Fig. 57 Precipitación areal anual nacional calculada con datos de diferentes décadas y calculada por el método de Vecinos Voronoi. - 500.0 1,000.0 1,500.0 2,000.0 2,500.0 3,000.0 3,500.0 RH01 RH02 RH03 RH04 RH05 RH06 RH07 RH08 RH09 RH10 RH11 RH12 RH13 RH14 RH15 RH16 RH17 RH18 RH19 RH20 RH21 RH22 RH23 RH24 RH25 RH26 RH27 RH28 RH29 RH30 RH31 RH32 RH33 RH34 RH35 RH36 RH37 Región Hidrológica P re ci p it ac ió n M ed ia A re al ( m m ) 1971-1980 1981-1990 1991-2000 1998-2007 Fig. 58 Comparación de la precipitación areal anual tomando datos de diferentes décadas y calculado por el método de Vecinos Voronoi con Orografía. 84 - 100.0 200.0 300.0 400.0 500.0 600.0 700.0 800.0 900.0 1,000.0 19 71 19 72 19 73 19 74 19 75 19 76 19 77 19 78 19 79 19 80 19 81 19 82 19 83 19 84 19 85 19 86 19 87 19 88 19 89 19 90 19 91 19 92 19 93 19 94 19 95 19 96 19 97 19 98 19 99 20 00 20 01 20 02 20 03 20 04 20 05 20 06 20 07 Año P re ci p it ac ió n M ed ia ( m m ) Fig. 59 Precipitación areal nacional anual por años calculada por Vecinos Voronoi con Orografía. 5.2. De los métodos de validación Cuando se propone un nuevo método de interpolación es deseable tener una medida de calificación del nuevo método. Específicamente en el ámbito de validación de métodos de interpolación, han sido utilizados diferentes enfoques. A saber: 5.2.1. Comparación contra datos artificiales Técnicas modernas como la obtención de superficies sin ruido a partir de datos dispersos y a través de interpolaciones hechas con redes neuronales [12] hacen la verificación de la calidad de sus modelos calculando el error de interpolación contrastando con datos artificiales, sin embargo esto tiene dos enfoques que son diametralmente opuestos al tipo de funciones que describen la precipitación areal. Primero está la asunción subyacente de que la función a ser interpolada es suave; por lo que el ruido se considera cualquier elemento que hace a la función no suave. Los datos artificiales que son usados para hacer la medida del error de interpolación son precisamente de una superficie matemáticamente suave. De hecho el promediar efectuado por cualquier método tiene como efecto suavizar la superficie. El conocimiento empírico y estadístico de la precipitación nos indica que no es necesariamente suave. Por ejemplo se pueden encontrar bosques de menos de 30 85 kilómetros cuadrados rodeados de zonas desérticas tal como el Parque Nacional El Chico que está en la montaña arriba de la cuidad de Pachuca en el estado de Hidalgo, México. El segundo está en la presunción de que se conoce a priori la función (que en nuestro caso sería una función de precipitación), y es con base en ello, o a menos un modelo de aproximación, que se generan los datos artificiales con los que se ha de medir el error de estimación. Este tipo de medición de error de estimación, basados en la medición de distancia entre un conjunto de datos artificiales conocidos a priori, luego muestreados para dar como entradas al modelo de interpolación a valorar y la interpolación así obtenida, sólo miden precisamente eso, qué tanto se parce la interpolación dada a la función original. Sin embargo no están considerados los siguientes aspectos:  Dado un conjunto de datos artificiales se escogen unas muestras para hacer con ellas la interpolación. Así, no se toma en cuenta la variabilidad de la interpolación cuando se cambia el conjunto de muestras; qué pasa si están concentradas en una zona, o agrupadas alrededor de unos polos de atracción, o espacialmente distribuidas siguiendo algún patrón. Es previsible que ninguna interpolación va a dar buenos resultados cuando no se tienen suficientes datos, aunque es difícil determinar cuántos son suficientes por la distribución espacial tridimensional de ellos y el desconocimiento de las frecuencias máximas de la función de precipitación en la naturaleza.  La sensibilidad a la variabilidad de cantidad de muestras dispersas utilizadas para hacer la interpolación.  La sensibilidad a la varianza de los valores de las muestras, es decir, cómo se comporta cuando los valores de las muestras son muy similares o no.  La resolución espacial máxima que puede tener la interpolación dada la resolución espacial específica de las muestras dispersas. Esto no debería de ser aplicable porque no se está haciendo un muestreo periódico, sin embargo al calcular sobre una malla regular se está haciendo un muestreo periódico de la superficie de precipitación. Otra manera de escribir esto es la siguiente: El proceso que se estudia es el siguiente: Dado un conjunto de datos artificiales se escogen unas muestras dispersas para hacer con ellas la interpolación. Luego se mide la distancia entre la interpolación y los datos artificiales originales conocidos a priori para dar una medida de calidad de la interpolación. Éste método de valoración de la interpolación no toma en cuenta la variabilidad debida a:  La distribución espacial de las muestras.  La cantidad de muestras utilizadas.  La similitud del valor de las muestras.  La precisión con la que se quiere obtener la interpolación. 5.2.2. Comparación por el método de validación cruzada de datos Modernos perfeccionamientos a la familia de interpolaciones conocida como Kriging tal como la utilización de la técnica llamada Máquina de Vectores de Soporte de Mínimos Cuadrados para la Interpolación de Kriging [13], que permite no tener que seleccionar un modelo de variograma 86 que describa con precisión la distribución espacial de los datos y que busca la adaptabilidad a las características diferentes de los conjuntos de datos, hacen la valoración de la calidad de su método a través de la técnica de validación cruzada. La validación cruzada consiste en tener un conjunto de datos originales, por ejemplo mediciones reales de un fenómeno, seleccionar un subconjunto y medir la distancia de la predicción contra el complemento de los datos. La comparación del desempeño de los métodos se hace para un conjunto de datos originales o medidos, con un mismo subconjunto se generan las interpolaciones con cada uno de los métodos y se comparan las distancias de cada uno al complemento del conjunto. Para obtener un índice global de desempeño del algoritmo independiente del número de muestras y de su dimensión se utilizan medidas globales como error medio, error medio absoluto o error medio cuadrado. Aún más, se puede hacer una experimentación sobre varios conjuntos de datos y obtener el promedio de las mediadas globales. Esta técnica, al igual que la anterior, tampoco toma en cuenta la variabilidad debida a:  La distribución espacial de las muestras.  La cantidad de muestras utilizadas.  La similitud del valor de las muestras.  La precisión con la que se quiere obtener la interpolación. 5.3. Validación comparando contra con datos artificiales El modelo utilizado para generar los datos artificiales fue el estándar mundial de cálculo de precipitación areal llamado Thiessen. Se realizaron los mapas de diferencia de los métodos contra éste y se muestran a continuación. En estos mapas se muestra la desviación positiva y negativa. Cuando el valor interpolado es igual al valor que da originalmente el método de Thiessen, entonces el color en el mapa es blanco. Los colores cercanos al blanco, aunque sean rojos o azules, quieren decir que la cercanía es mucha. Ahora, si la diferencia es negativa, es decir se está subestimando el valor entonces en el mapa se representa con color azul. Los valores mostrados son porcentuales debido a que se requiere una normalización para hacer una comparación objetiva. No es lo mismo tener una variación de 10 mm cuando el valor debería se 8 mm, que tener una variación de igualmente 10 mm cuando el valor debería de ser 4,950 mm. Por ello se usa una variación porcentual, está normalizado de manera local al valor en cada punto en la malla. Entonces cuando la variación es del 50% o mayor, entonces en el mapa se muestra en color azul fuerte. Para variaciones positivas del 50% o mayores, es decir, sobreestimaciones, el mapa se muestra en color rojo. En la Fig. 60 se muestra el error de Vecinos Voronoi contra Thiessen. En todos los mapas de error mostrados se ven subyacentes los polígonos de Voronoi del método de Thiessen y esto es inherente, ya que el error es sólo una resta entre dos funciones. 87 Fig. 60 Error del método Vecinos Voronoi contra Thiessen. Fig. 61 Error del método Vecinos Voronoi con Orografía contra Thiessen. En la Fig. 61 se muestra el error del método Vecinos Voronoi con Orografía contra el método Thiessen. Se ve el efecto suavizado de la orografía. 88 Fig. 62 Error del método 4 Estaciones Cercanas por el inverso de la distancia comparado contra Thiessen. Fig. 63 Error del método 16 Estaciones Cercanas por el inverso de la distancia comparado contra Thiessen. En la Fig. 62 se muestra el error del método de 4 Estaciones Cercanas por el inverso de la distancia contra Thiessen y se ve claramente que es muy parecido este método, que difiere más si se pondera por el inverso de la distancia al cuadrado, mostrado en la Fig. 64. 89 Fig. 64 Error del método 4 Estaciones Cercanas por el inverso de la distancia al cuadrado comparado contra Thiessen. Fig. 65 Error del método 16 Estaciones Cercanas por el inverso de la distancia al cuadrado comparado contra Thiessen. En la Fig. 63 se muestra el error del método de 16 Estaciones Cercanas por el inverso de la distancia y en la Fig. 65 se muestra el error de la versión por el inverso de la distancia al cuadrado. Es muy claro el suavizamiento por la velocidad de decaimiento de la función. 90 Fig. 66 Error del método 4 Estaciones Cercanas con Orografía comparado contra Thiessen. Fig. 67 Error del método 16 Estaciones Cercanas con Orografía comparado contra Thiessen. En la Fig. 66 se muestra el error del método 4 Estaciones Cercanas con Orografía y en la Fig. 67 se muestra el de 16 Estaciones Cercanas con Orografía. En ambos se aprecia cómo las diferencias aumentan en las zonas montañosas y son casi cero en las zonas planas. Esto mismo puede verse en la Fig. 68 donde se compara el método contra su versión con orografía. 91 Fig. 68 Error del método Thiessen con Orografía comparado contra Thiessen. Fig. 69 Error del método Kriging Gaussiano comparado contra Thiessen. Es claro en la Fig. 69 que el método de Kriging es capaz de producir grandes desviaciones de los valores. Siendo un método global los errores se extienden en todo el territorio. Al comparar con el otro método global, el mostrado en la Fig. 70 es claro por qué se prefiere el método Spline de Curvatura Continua sujeto a Tensión. 92 Fig. 70 Error del método Spline de Curvatura Continua Sujeto a Tensión con T=0.75 comparado contra Thiessen. 5.4. Validación usando referencia cruzada Evaluación del desempeño de los métodos de interpolación usando referencia cruzada. Tabla 8 Tabla de evaluación del desempeño de los métodos de interpolación usando referencia cruzada. Método Máxima diferencia porcentual en el valor interpolado. Mínima diferencia porcentual en el valor interpolado. Promedio de la diferencia porcentual del valor interpolado. Varianza de la diferencia porcentual del valor interpolado. Thiessen 455.0% 0.1% 26.0% 24.4% 4EstInvDist 442.5% 0.1% 23.8% 21.8% 4EstInvDist2 638.5% 0.3% 37.7% 80.4% 16EstInvDist 421.1% 0.1% 26.6% 19.5% 16EstInvDist 2 426.1% 0.1% 24.4% 19.7% Voronoi 381.8% 0.3% 24.2% 17.4% VoronoiDifAlt 365.4% 0.1% 37.0% 19.6% ThiessenDifA lt 607.2% 0.2% 41.4% 38.9% 4EstDifAlt 454.4% 0.2% 37.4% 24.2% 16EstDifAlt 350.1% 0.3% 40.6% 21.0% KrigGauss 2518.8% 0.3% 79.6% 497.7% 93 KrigExp 2541.8% 0.2% 83.4% 508.2% Laplaciano 518.0% 0.1% 45.8% 28.2% Biarmónico 363.5% 0.2% 23.5% 15.9% TensorLaplBi arm 338.4% 0.1% 24.1% 14.6% KrigEsferico 1173.7% 0.2% 67.0% 152.3% Para obtener los datos de la Tabla 8 se generaron para cada renglón seis mapas, tres extrayendo el 20% de las estaciones de manera aleatoria y luego midiendo la diferencia de la interpolación contra el dato verdadero. Otros tres extrayendo el 30% de las estaciones. Fig. 71 Tres experimentos con el 20% de las estaciones omitidas y 3 con el 30% para medir los valores interpolados contra las mediciones reales para el método Vecinos Voronoi con Orografía. 94 5.5. Validación comparando contra con isoyetas manuales históricas Dentro de la CONAGUA el último proyecto de generación de isoyetas a nivel nacional de datos históricos fue el realizado para los datos de 1930 a 1990, del cual se conservan una libreta de mapas grande en papel. En este documento se tienen mapas por estados a diferentes escalas y no se tiene mapa nacional, aunque los cálculos nacionales sí están plasmados en él. En la Fig. 72 puede verse cómo tanto la precipitación areal calculada para cada una de las Regiones Hidrológicas se parece mucho a la calculada por el método Vecinos Voronoi con Orografía. Los datos y detalle de esta gráfica pueden consultarse en la Tabla 9. Una medida de verificación de equiparabilidad entre los datos nos la da la comparación de los porcentajes de área de las diferentes Regiones Hidrológicas. Casi todas son muy parecidas, excepto unas cuantas como la 4 y la 8 que han sido históricamente redefinidas por diferentes criterios hidrológicos y administrativos. El detalle de estos datos puede verse en la Tabla 9. Comparación de datos manuales CONAGUA 1930-1990 contra Vecinos Voronoi con Orografía 0% 50% 100% 150% 200% 250% RH01 RH03 RH05 RH07 RH09 RH11 RH13 RH15 RH17 RH19 RH21 RH23 RH25 RH27 RH29 RH31 RH33 RH35 RH37 Regiones Hidrológicas V al o r re la ti vo a l m an u al % valor % área Fig. 72 Comparación de los datos de precipitación y área por Regiones Hidrológicas calculada por el método de Vecinos Voronoi contra los datos manuales CONAGUA 1930 – 1990. 100% quiere decir que son iguales. 95 Tabla 9 Tabla de datos donde se comparan la precipitación areal por Región Hidrológica y las áreas de las mismas del método Vecinos Voronoi contra los datos manuales CONAGUA 1930 – 1990. 100% quiere decir que son iguales. DATOS MANUALES CONAGUA 1931-1990 Comparación contra resultados Vecinos Voronoi con orografía 1971-1980 malla 2km con 1 CPU Comparación Región Hidrológica Área Lámina (mm) zona valor area Región Hidrológica % valor % área RH01 26,785 197 RH01 309.78 26757.83 RH01 157% 100% RH02 42,990 102 RH02 151.25 43020.63 RH02 148% 100% RH03 29,563 166 RH03 171.4 29072.76 RH03 103% 98% RH04 6,741 134 RH04 207.35 15223.6 RH04 155% 226% RH05 12,898 101 RH05 156.4 13022.99 RH05 155% 101% RH06 11,590 226 RH06 187.99 12139.68 RH06 83% 105% RH07 15,963 131 RH07 274.22 10776.88 RH07 209% 68% RH08 56,383 227 RH08 332.75 56487.66 RH08 147% 100% RH09 144,881 481 RH09 524.84 139023.76 RH09 109% 96% RH10 104,790 790 RH10 844.97 104744.99 RH10 107% 100% RH11 51,837 841 RH11 921.34 53059.95 RH11 110% 102% RH12 137,144 695 RH12 748.58 135869.97 RH12 108% 99% RH13 5,072 1354 RH13 1199.14 5558.53 RH13 89% 110% RH14 11,664 1052 RH14 892.56 12222.97 RH14 85% 105% RH15 13,567 1093 RH15 1385.38 13756.58 RH15 127% 101% RH16 17,863 875 RH16 1021.96 17869.3 RH16 117% 100% RH17 9,196 1008 RH17 902.42 9222.7 RH17 90% 100% RH18 115,895 954 RH18 1077.99 118519.75 RH18 113% 102% RH19 13,090 1250 RH19 1074.32 13166.48 RH19 86% 101% RH20 39,994 1411 RH20 1039.58 40313.82 RH20 74% 101% RH21 10,148 1309 RH21 1003.62 10656.31 RH21 77% 105% RH22 16,443 1020 RH22 960.65 15936.53 RH22 94% 97% RH23 12,347 2382 RH23 2758.1 12099.23 RH23 116% 98% RH24 247,163 424 RH24 443.51 231966.06 RH24 105% 94% RH25 55,060 668 RH25 746.16 55259.62 RH25 112% 100% RH26 94,073 905 RH26 822.51 98534.06 RH26 91% 105% RH27 25,789 1577 RH27 1236.16 26304.23 RH27 78% 102% RH28 61,189 1818 RH28 1550.39 58149.73 RH28 85% 95% RH29 30,187 2292 RH29 2303.8 30140.42 RH29 101% 100% RH30 106,210 1903 RH30 1872.09 102346.15 RH30 98% 96% RH31 16,544 1199 RH31 1167.44 24541.85 RH31 97% 148% RH32 65,143 1073 RH32 1067.68 57710.94 RH32 100% 89% RH33 35,431 1245 RH33 1222.33 39207.29 RH33 98% 111% RH34 84,185 348 RH34 412.91 90178.37 RH34 119% 107% RH35 56,883 304 RH35 313.17 64307.27 RH35 103% 113% RH36 89,239 387 RH36 384.67 92082.51 RH36 99% 103% RH37 94,243 434 RH37 429.25 88489.18 RH37 99% 94% De estas comparaciones se concluye que el método Vecinos Voronoi con Orografía es un método que aproxima los mapas que antes se lograban mediante cálculos manuales y criterio experto y que hasta ahora no eran reproducibles por ningún método computacional automatizado. 96 5.6. Validación cualitativa contra imágenes de satélite En inspección visual al usar el programa Google Earth© y sobreponerle el mapa interpolado por Vecinos Voronoi con Orografía como se muestra en Fig. 73 a las imágenes de satélite en espectro de luz visible que este programa provee encontramos primero que hay una buena correlación entre los bordes de nuestros mapas con los del programa. También, como se puede apreciar en Fig. 74 existe una gran correlación entre las intensidades de precipitación definidas por la interpolación de Vecinos Voronoi con Orografía con la vegetación subyacente. Hacer un análisis no visual sino determinístico, calculado, contra la vegetación del país está fuera del alcance de este trabajo debido a la falta de acceso a mapas de vegetación detallados a la suficiente resolución. Fig. 73 Interpolación por Vecinos Voronoi con Orografía sobrepuesta a las imágenes de satélite de Google Earth©. Fig. 74 Acercamiento que muestra la enorme correlación entre la intensidad de precipitación y la vegetación subyacente en las imágenes satelitales. 97 6. Conclusiones Los métodos globales de interpolación de precipitación difieren más de la precipitación real debido a que tratan de compensar sobre toda la superficie fenómenos de clima local. La precipitación es mayormente un fenómeno local y múltiple, por ejemplo la debida a procesos de convección atmosférica. Fenómenos relativamente globales o que abarcar mayor territorio se presentan sólo bajo condiciones de huracán. Los métodos que interpolan a través de una superposición lineal de funciones, tales como Kriging y Biarmónico-Laplaciano no toman en cuenta las condiciones locales de cuencas hidrológicas determinadas por la orografía. Sobre grandes extensiones de terreno y para periodos largos de consideración, desde anuales a normales (OMM, 30 años), la precipitación areal de todo el territorio puede calcularse casi por cualquier método de los aquí expuestos, ya que los fenómenos de compensación de valores promedian los errores y todos convergen a casi un mismo valor. El mejor método para calcular la superficie de precipitación es el de Vecinos Voronoi con Orografía porque es capáz de reproducir con cercanía los mapas que se hacían manualmente con criterio experto. El trazo de las isoyetas puede hacerse perfectamente y sin ambigüedades, y en tiempo lineal con marching triangles. La implementación a través de una función recursiva es altamente eficiente. La utilización de las bibliotecas BLAS y LAPACK® permite hacer uso del cómputo en paralelo, para este tipo de problemas, con muy poco esfuerzo de programación, de manera rápida y obteniendo grandes beneficios. La aceleración de la inversión de la matriz requerida por Kriging fue de 200 veces con GPUs, lo cual permite que un cálculo que tarda días tarde sólo segundos. Las pruebas demuestran que la tecnología de paralelización múltiple y optimización del lenguaje ACML de la compañía AMD es muy superior a la de CUDA de la compañía nVIDIA©. 98 Trabajo futuro La incorporación de más información geofísica al modelo de interpolación pareciera que podría aumentar la precisión de la estimación a nivel local. Modelado de microclimas incorporando mediciones climatológicas podría establecer funciones locales de ajuste de interpolación. Mapas de vegetación podrían también establecer un parámetro para la interpolación, sin embargo no existen para todo el territorio mexicano con una precisión y resolución adecuada para hacer una consideración seria. 99 Anexo - Siglas, acrónimos, nombres y marcas MPI Protocolo de paso de mensajes que sirve para hacer cómputo en paralelo. IDE Ambiente de desarrollo integrado. C/C++ Dos lenguajes de programación de computadoras. Lenguaje C. Lenguaje C++. DeinoMPI Implementación del marco de trabajo MPICH 2 para Windows. Windows Sistema operativo de computadoras. Le pertenece a la compañía Microsoft. OpenGL Biblioteca de Gráficos Abierta. CONACYT Consejo Nacional de Ciencia y Tecnología. CUDA Arquitectura de cómputo en paralelo de la compañía nVIDIA©. GNU General Public License Acuerdo de licenciamiento de programas de cómputo en el que se permite el estudio, utilización y modificación de los programas y sus códigos fuentes siempre que se exprese en todo momento el acuerdo y se permitan los mismos derechos a terceros. GCC Colección de compiladores GNU. Dev-C++ IDE para programar en los lenguajes C y C++. Software libre. Bloodshed Dev- C++ IDE para programar en los lenguajes C y C++. Software libre. GPU Unidad de Procesamiento Gráfica. IIMAS Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas. Isolínea Línea que describe un isovalor a través de una superficie. Isoyeta Isolínea sobre una superficie de precipitación. Isocurva Curva que describe un solo valor específica a través de una superficie. Curva de nivel Línea que describe un sólo valor específico a través de una superficie. También es conocida como Isolínea. Se estila llamar curva de nivel a la isolínea sobre un mapa de altitudes. Win32 Sistema operativo Windows de 32 bits. CONACYT Consejo Nacional de Ciencia y Tecnología. México. UNAM Universidad Nacional Autónoma de México. MCC Maestría en Ciencia e Ingeniería de la Computación. MPICH 2 Implementación de los estándares 1 y 2 de pasado de mensajes. GLUT Graphics Library Utility Toolkit. http://www.opengl.org/resources/libraries/glut.html NSF National Science Foundation, USA. ACML AMD Core Math Library. Biblioteca de matemáticas esenciales de AMD. AMD Compañía que produce procesadores de computadora y otras cosas. Intel Compañía que produce procesadores de computadora y otras cosas. Tesselación Partición de la totalidad de un espacio con objetos que no se traslapan ni dejan espacios entre sí. Triangulación Teselación de un espacio en triángulos. Open manifold Un manifold no compacto sin frontera. Un manifold es un espacio topológico que a una escala suficientemente pequeña semeja un espacio euclidiano de dimensión específica menor al del espacio topológico. OpenCL Open Computing Language. PCIC Posgrado en Ciencia e Ingeniería de la Computación. Surfer© Programa de visualización de datos científicos de la empresa Golden Software. http://www.goldensoftware.com/ 101 Referencias [1] V. Ponce, Engineering Hydrology, EUA: Prentice Hall, 1989. [2] W.C.M.V. Beers and J.P.C. Kleijnen, “Kriging Interpolation in Simulation: A Survey,” Proceedings of the 2004 Winter Simulation Conference, EUA, 2004. [3] P. Bedient and W. Huber, Hidrology and Floodplain Analysis, EUA: Addison-Wesley, 1992. [4] WMO, “World Weather / Climate Extremes Archive http://wmo.asu.edu,” 2010. [5] H. Späth, Two Dimensional Spline Interpolation Algorithms, Natick, EUA: A. K. Peters, Ltd., 1995. [6] I.S. University, “Mountain Simulation http://www.pals.iastate.edu/simulations/Mtnsim/index.html,” 2010. [7] W.H.F. Smith and P. Wessel, “Gridding with Continuous Curvature Splines in Tension,” Geophysics, vol. 55, 1990, p. 293. [8] D.S. Wilkis, Statistical Methods in the Atmospheric Sciences, Academic Press Professional, Inc., 2006. [9] Z. Sokol and B. Bliznak, “Areal Distribution and Precipitation-Altitude Relationship for Heavy Short-Term Precipitation in the Czech Republic in the Warm Part of the Year,” Atmospheric Research, vol. 94, 2009, pp. 652-662. [10] METI, “Modelo Digital de Elevación Mundial ASTER NASA http://asterweb.jpl.nasa.gov,” 2009. [11] L.M. De La Cruz, “MQ-RBF Meshless Method for solving CFD problems using an Object-Oriented Approach,” SCAT Mobility Grant Report, vol. 1, 2008, pp. 1-41. [12] X. Liu, W. Xu, J. Xu, and Y. Guan, “Scattered Points Denoising of TC-Bézier Surface Fitting,” Proceedings of the 2010 Third International Conference on Knowledge Discovery and Data Mining, 2010. [13] W. Huizan and Z. Ren, “Improved Kriging Interpolation Based on Support Vector Machine and Its Application in Oceanic Missing Data Recovery,” Proceedings of the 2008 International Conference on Computer Science and Software Engineering, EUA: IEEE Computer Society, 2008, pp. 726-729. [14] V.B. Anand, Computer Graphics and Geometric Modeling for Engineers, EUA: Willey, 1993. [15] J. Boissonnat and D. Cohen, Effective Computational Geometry for Curves and Surfaces, Berlin, Springer Verlang, 2008. 102 [16] E.J. Candes and M.B. Wakin, “An Introduction to Compressive Sampling,” IEEE Signal Processing Magazine, vol. 25, 2008, pp. 21-30. [17] S. Cunningham, Programming in OpenGL for Visual Communication Graphics, EUA: Pearson. Prentice Hall, 2007. [18] F.C. Curriero, “On the Use of Non-Euclidean Distance Measures in Geostatistics,” Mathematical Geology, vol. 38, 2006, pp. 907-926. [19] A. Davies and P. Samuels, An Introduction to Computational Geometry for Curves and Surfaces, EUA: Oxford University Press, Inc., 1996. [20] H. Edelsbrunner, Geometry and Topology of Mesh Generation, EUA: Cambridge University Press, 2001. [21] K.E. Kerry and K.A. Hawick, Kriging Interpolation on High-Performance Computers Technical Report DHPC-035, Australia: 1998. [22] E. Langetepe and G. Zachmann, Geometric Data Structures for Computer Graphics, Alemania, A. K. Peters, Ltd., 2006. [23] NVIDIA©, “CUDA - Compute Unified Device Architecture http://www.nVIDIA©.com/cuda,” 2010. [24] H. Nyquist, “Certain Topics in Telegraph Transmission Theory,” Proceedings of the IEEE, vol. 90, EUA, 2002, pp. 280-305. [25] J. Reif and S. Rajasekaran, Handbook of Parallel Computing: Models, Algorithms, and Applications, EUA: Chapman & Hall / CRC, 2006. [26] J. Romberg, Teorema de Nyquist http://cnx.org/content/m12971/1.2/, 2005. [27] S. Sakata, “An Efficient Algorithm for Kriging Approximation and Optimization with Large-Scale Sampling Data,” Computer Methods in Applied Mechanics and Engineering, vol. 193, EUA, 2004, pp. 385-404. [28] C.E. Shannon, “A Mathematical Theory of Communication,” The Bell System Technical Journal, vol. 27, EUA, 1948, pp. 379-423, 623-656. [29] B. Su and D. Liu, Computational Geometry: Curve and Surface Modeling, EUA: Academic Press Professional, Inc., 1989. [30] W.G. Tuller, Theoretical Limitations on the Rate of Transmission of Information, Technical Report No. 114 MIT, EUA: 1949. [31] X. Weixiang and W. Liuqiang, “Quadratic TC-Bézier Curves with Shape Parameter,” Advanced Materials Research, vol. 179-180, EUA, 2011, pp. 1187-1192. 103 Índice 16 Estaciones Cercanas con Orografía, 71 16 Estaciones Cercanas por el inverso de la distancia, 69 16 Estaciones Cercanas por el inverso de la distancia al cuadrado, 70 1-vecindad, 42 4 Estaciones Cercanas con Orografía, 71 4 Estaciones Cercanas por el inverso de la distancia, 69 4 Estaciones Cercanas por el inverso de la distancia al cuadrado, 70 aceleración, 80, 81 ACML©, 58, 63 algoritmos, 16 altitudes, 34 AMD©, 63 área del polígono, 47 áreas totales, 77 base de datos climatológica, 58 Bessel, 38 Biarmónico, 29 bloque mínimo, 38 Cálculo de la distancia, 41 Cálculo de malla, 38 Cálculo por malla, 56 Cálculos areales, 47 cartografía digital, 54 celda, 36 Cercanos, 44 cobertura de la medición, 4 CONAGUA, 18, 58 conjunto de datos de buena calidad, 35 conjunto de datos en el intervalo a consideración, 34 constante de calidad, 35 constante de orografía, 49 cosenos, 47 CPUs, 63 Criterio de cercanía, 40 CUDA©, 58, 63 curvatura mínima, 27 datos artificiales, 86, 88 datos dispersos, 17 Datos dispersos, 37 de modelos termodinámicos de la atmósfera, 18 décadas, 85 Definición del problema, 33 Delphi, 58 DEM, 55 desempeño de los métodos, 94 diferencia de altitud, 41, 45 distancia geodésica, 38 dominio de las posiciones geográficas, 35 dominio de las posiciones geográficas considerando la orografía, 35 104 elipsoide Clarke 80, 37 época de lluvias, 78 Escala de medición, 4 Esférico, 25 estaciones, 10 Estaciones Cercanas, 44 Estaciones Cercanas con Orografía, 45 estaciones convencionales, 16 estaciones del desierto, 7 estaciones pluviométricas, 7, 53 Estados, 84 Experimentos, 65 Exponencial, 24 Extrapolación, 16 factor de precipitación por altitud, 46 función de acumulación, 35 función de calidad, 35 funciones de base radial, 62 Gaussiano, 23 GPUs, 53, 63 gráfica de contornos, 3 Hipótesis, 11 imágenes de satélite, 98 intensidad de lluvia media, 9 interioridad, 39 Interpolación, 16 intervalo a consideración, 34 isolínea, 15 isoyetas, 2, 58, 61, 73 isoyetas manuales, 96 Japón, 54 Kriging, 2, 20, 21, 37, 63 Kriging ordinario, 21 Kriging simple, 21 Kriging universal, 21 LAPACK©, 58, 63 Laplaciano, 29 latitudes, 33 lattice, 36 Línea de contorno, 15 lluvia, 3 longitudes, 33 malla, 56 malla de resultados, 36 Mapas actuales, 17 Mapas de isoyetas, 15 mapas de precipitación areal, 20 marching triangles, 48 matriz de varianzas y covarianzas, 63 Metodología, 12 Métodos contra los que se compara, 66 métodos de cálculo, 58 105 métodos de interpolación, 66 métodos de validación, 86 Métodos propuestos, 66 Modelo digital de elevación, 54 msnm, 50 NASA, 54 nVidia©, 63 Objetivos, 12 Organismos de Cuenca, 83 orografía, 18, 58 Parámetros, 59 periodo, 16 periodos, 34 polígono de Voronoi, 42 posición geográfica, 33 precipitación, 18 precipitación acumulada, 35 precipitación areal, 2, 3, 15, 77 Precipitación areal anual, 83, 86 precipitación areal nacional, 76 precipitación diaria, 37, 53 precipitación máxima, 8 precipitaciones areales, 53 predicción, 16 preprocesamiento, 55, 58 primera vecindad, 42 primeros vecinos, 41 Procesamiento, 61 Procesamiento en paralelo, 62 Puntos de cálculo, 38 Puntos interiores, 39 radio de influencia, 37 radios de competencia, 62 Rastreo de tiempos, 75 reconstrucción, 16 referencia cruzada, 94 Regiones Hidrológicas, 57, 78, 83 regla de acumulación, 35 resolución, 36, 56, 57 resoluciones, 77, 78 Resultados, 65 rugosidad, 22 secas, 78 segmento, 39 serie, 80 Servicio Meteorológico Nacional, 18 Smith, 38 Spline CCST, 37 spline de curvatura continua sujeta a tensión, 27 Spline de curvatura continua sujeto a tensión, 20 Spline de Curvatura Continua Sujeto a Tensión, 2 superficie de precipitación, 17, 36 106 territorio, 35 Thiessen, 2, 37 Thiessen con Orografía, 46, 72 tiempos de proceso, 81 Tiempos de proceso, 82 tiff, 55 trazado de las isoyetas, 48 ubicación de las estaciones, 74 validación cruzada de datos, 87 valor de la celda, 47 valor normal de precipitación, 4 valores normales, 3 Vecinos Voronoi, 41, 67 Vecinos Voronoi con Orografía, 43, 68 zonas del territorio, 35