Estimación de la radiación solar mediante técnicas de geoestadística utilizando una red de sensores distribuidos

i Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Estimación de la radiación solar mediante técnicas de geoestadística utilizando una red de sensores distribuidos Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Autor: Mario Fernández Rangel Tutor: Amparo Núñez Reyes Sevilla, 2019 ii iii Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Estimación de la radiación solar mediante técnicas de geoestadística utilizando una red de sensores distribuidos Autor: Mario Fernández Rangel Tutora: Amparo Núñez Reyes Profesora titular Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 iv v Trabajo Fin de Grado: Estimación de la radiación solar mediante técnicas de geoestadística utilizando una red de sensores distribuidos Autor: Mario Fernández Rangel Tutora: Amparo Núñez Reyes El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2019 El Secretario del Tribunal vi vii A mi familia viii ix Agradecimientos Quiero darle las gracias a mi tutora, Amparo Núñez Reyes, sin su ayuda este Trabajo de Fin de Grado no hubiera podido salir adelante. También estoy muy agradecido a todas las personas que me han ayudado a lo largo de este grado, sus consejos han sido muy valiosos para llegar hasta este momento. Darle las gracias a mis amigos que siempre han estado conmigo durante este tiempo y ayudarme a aliviar cuando la carga era demasiado excesiva para soportarla. Finalmente, agradecer a mi familia todo el ánimo que he recibido. Quiero acabar con una cita de mi escritor favorito: “The most important step a person can take is always the next one.” Brandon Sanderson Mario Fernández Rangel Alumno del Grado en Ingeniería en Tecnologías Industriales Sevilla, 2019 xvi xvii ÍNDICE DE FIGURAS Figura 1. Distribución de la red de sensores (1). 3 Figura 2. Arquitectura del protocolo de comunicación de una WSN (1). 5 Figura 3. Ejemplo de variograma experimental (7). 7 Figura 4. Variograma experimental con modelo lineal. 9 Figura 5. Mapa del desierto de Tabernas (11). 11 Figura 6. Ubicación de los nodos de la red de sensores simulada. 12 Figura 7. Interfaz de usuario del CAMS Radiation Service. 13 Figura 8. Visualización de los datos descargados. 14 Figura 9. Nube de puntos. 15 Figura 10. Gráfica generada mediante Desc del mes de junio. 16 Figura 11. Gráfica q-q para el mes de junio. 17 Figura 12. Distribución de los valores de los datos de los días de junio 17 Figura 13. Gráfica q-q para los datos de días de junio 18 Figura 14. Histograma del mes de junio 19 Figura 15. Mapa con los valores de radiación de los nodos. 19 Figura 16. Zonas UTM de España (20). 20 Figura 17. Variograma experimental con la muestra de datos del mes de junio. 21 Figura 18. Direcciones para el variograma del mes de junio. 22 Figura 19. Variograma con modelo teórico del mes de junio. 23 Figura 20. Variogramas direccionales junto al modelo del variograma del mes de junio. 23 Figura 21. Variograma con el modelo de b) 24 Figura 22. Variograma con el modelo de c) 25 Figura 23. Variograma con el modelo de d) 25 Figura 24. Variogramas de los distintos meses del año 2015. 26 Figura 25. Predicción de las horas con el modelo de variograma del mes de junio. 27 Figura 26. Desviación típica de la predicción de horas con el modelo de variograma de junio. 28 Figura 27. Resultado de la validación cruzada con el modelo de junio del variograma. 29 Figura 28. Comparación de las estimaciones con datos normalizados y sin normalizar. 30 Figura 29. Comparación de las predicciones de la hora de los distintos modelos de variograma. 31 Figura 30. Puntos de validación sobre la predicción. 32 xviii xix Notación WSN Wireless Sensor Network (Red inalámbrica de sensores) ZEPA Zona Especial de Protección para las Aves PVGIS Photovoltaic Geographical Information System TOA Radiación en plano horizontal por encima de la atmósfera GHI Radiación global en plano horizontal a nivel del suelo BHI Radiación directa en plano horizontal a nivel del suelo DHI Radiación difusa en plano horizontal a nivel del suelo BNI Radiación directa en plano móvil perpendicular a la incidencia solar CPV CRS Concentrador fotovoltaico Coordinates Reference System UTM Coordenadas universales transversal de Mercator EPSG European Petroleum Survey Group RMSE Root Mean Square Error MAE Mean Absolute Error R2 Coeficiente de determinación. IDW Inverse Distance Weighting 1 1 INTRODUCCIÓN a creciente preocupación por las graves consecuencias producidas por el calentamiento global, está conduciendo a la sociedad, a un progresivo endurecimiento sobre el control normativo de las emisiones de los gases de efecto invernadero a la atmósfera. Surgiendo cada año nuevas normas locales, nacionales, europeas e internacionales, como el protocolo de Kyoto que comprende dos períodos. El primero desde 2008 hasta 2012 y el segundo y actual desde 2013 hasta 2020. Siendo el principal objetivo de dicho acuerdo internacional, la reducción de las emisiones de seis gases de efecto invernadero relacionados directamente como causantes del calentamiento global, uno de los mayores problemas a los que se enfrenta el mundo globalizado actualmente. Por todo ello se fomenta y justifica el uso de las energías renovables. Dentro de este contexto, se encuentra el presente trabajo, donde se estudia la variabilidad espacial de la radiación solar de una determinada zona del desierto de Tabernas para posteriormente realizar la interpolación espacial, estimando la radiación solar en los lugares donde no se ubiquen sensores, ya sea por motivos geográficos o por costes económicos. El escenario hipotético en el que se realiza el estudio consiste en una planta solar ubicada en el desierto de Almería. Dicha planta puede ser de cualquier tipo y dispone de una red de sensores distribuidos irregularmente por todo el terreno. La planta presenta colectores solares que pueden ser controlados automáticamente mediante un controlador avanzado de procesos. Para ello es primordial conocer la radiación solar por toda la zona en la que se encuentran los paneles solares. Los datos se han obtenido de satélites disponibles de forma libre y se ha simulado una red de sensores de 49 nodos distribuidos irregularmente. Con el uso de técnicas de geoestadística, el análisis estructural de los datos ha permitido obtener información sobre la variabilidad espacial de la radiación solar, utilizándose para la estimación posterior de la radiación solar en diferentes periodos de tiempo. No se pretende entrar en profundidad en los desarrollos teóricos de los métodos y técnicas de la geoestadística, el principal objetivo consiste en lograr aplicar las técnicas mediante programas informáticos realizados mediante los entornos de programación R y Matlab, así como estudiar los resultados obtenidos para comprobar la viabilidad que puedan tener en entornos reales. • Organización de la memoria. El documento se ha dividido en las siguientes partes: • Redes de sensores inalámbricos. En el primer capítulo, se hace una descripción de esta tecnología y las complejidades de su diseño, los problemas que presenta actualmente y las investigaciones actuales en curso para su mejora con el fin de resolverlos. • Geoestadística. Se explica los conceptos importantes de la geoestadística necesarios para comprender su funcionamiento y realizar las estimaciones mediante kriging, el método elegido para estimar la radiación solar. • Origen y análisis de los datos. Donde se explica en que consiste nuestro trabajo, elegimos el lugar y describimos la búsqueda de los datos de radiación solar, realizamos el análisis necesario de los datos y en último lugar, explicamos en que consiste la reproyección del sistema de coordenadas. • Variabilidad espacial y estimación de la radiación solar. Dentro de este capítulo se muestran los resultados de las investigaciones realizadas en el trabajo que incluyen: Representación del variograma experimental, ajuste de los modelos teóricos para el variograma, la obtención de las estimaciones de los valores de radiación solar mediante kriging y los cálculos de los errores cometidos. Validaremos el resultado de nuestro modelo mediante validación cruzada y el grado de ajuste que tenemos a nuestros datos. L Introducción 2 2 • Conclusiones. Aquí se exponen las conclusiones más importantes obtenidas sobre el trabajo realizado y proponemos líneas de investigación futuras para desarrollar este trabajo y lograr nuevos objetivos o mejoras posibles. Finalmente, incluimos un anexo con todo el código desarrollado para el trabajo. Disclaimer: Contains modified Copernicus Atmosphere Monitoring Service Information 2019. Neither the European Union nor ECMWF is responsible for the use of the information 3 2 REDES DE SENSORES INALÁMBRICAS as redes de sensores inalámbricas es un tema de actualidad en la investigación científica, todavía falta un tiempo para que su uso se extienda en diversas aplicaciones, pero probablemente sea una tecnología que tendrá mucha incidencia en el futuro. 2.1. Redes de sensores inalámbricas Una red de sensores inalámbricos, conocida en inglés como Wireless sensor network (WSN), se define como un grupo de nodos distribuidos espacialmente conectados de forma inalámbrica, normalmente cuentan con un terminal común encargado de recoger toda la información que obtienen los sensores colocados en los nodos y trabajan de forma cooperativa para mandar la información obtenida. Su función suele ser monitorizar y/o controlar una variable ambiental como presión, temperatura, contaminantes, movimiento, niveles de ruido, etc. Es una tecnología reciente y en investigación, por lo que sus aplicaciones todavía están por desarrollar, en (1 págs. 393-422), se repasan sus posibles aplicaciones, que pueden incluir: • Aplicaciones militares. • Aplicaciones en el medio ambiente: detección de incendios, mapeado del medio, control de especies animales, detección de inundaciones. • Aplicaciones en la salud: Monitorización de datos fisiológicos, rastreo de pacientes y médicos. • Aplicaciones en hogares: Monitorización de variables ambientales, automatización del hogar. • Otros usos comerciales: Automatización de fábricas, control de inventarios En (2 págs. 2292-2320) se clasifican las WSN en dos tipos, estructuradas y desestructuradas. Las desestructuradas contienen muchos nodos de sensores y se pueden colocar de manera aleatoria, se usan para funciones de monitorización, aunque su mantenimiento es costoso. Las redes estructuradas permiten colocar menos nodos y tienen menor coste de mantenimiento, aunque su colocación puede ser más costosa. Para ilustrar de forma gráfica la WSN, podemos visualizar este esquema: Figura 1. Distribución de la red de sensores (1). L Redes de sensores inalámbricas 4 4 El diseño de la WSN depende de factores como la tolerancia de la tolerancia, la escalabilidad, los costes de producción, las limitaciones del hardware, la topología de la red, el medio, el canal de transmisión y el consumo de energía. En la actualidad, una de las mayores restricciones de las WSN puede ser el consumo de energía debido a que el método más usual de suministro ha sido mediante baterías que se agotan en el tiempo, por eso es necesario reducir el consumo en todos los pasos del diseño. Si fuera posible suministrar una energía renovable como la energía solar, no existiría esta limitación, por eso la instalación en una planta fotovoltaica, que además pueda suministrar energía a la WSN que permita mejorar la eficiencia en la obtención de energía sería ideal. Otro de los problemas es que los nodos de sensores sean capaces de organizarse a sí mismos en una red y ser capaz de manejarse de forma eficiente (2). Normalmente, las WSN se colocan en la región objetivo sin una infraestructura de apoyo (3), los propios nodos toman la información que captan los sensores y se mandan al nodo “getaway” o “sink” en la que se vuelca toda la información y que suele estar conectado a Internet, mediante este nodo llegan los datos al terminal conectado mediante Internet, llamado estación base, desde la que un usuario puede controlar o monitorizar la información generada. Debido a que los WSN suelen ser colocados sin una infraestructura de soporte, necesitan tener un protocolo de comunicación para transmitir la información de forma inalámbrica que suele incluir señales y detección de error como los protocolos usados habitualmente por los equipos informáticos, por ejemplo, el protocolo TCP. Este protocolo debe ser diseñado con el objetivo de reducir el consumo de energía de los nodos, transmitir de forma eficiente la información y promover el esfuerzo cooperativo de los nodos (1). Las capas del protocolo se describen de la siguiente forma (1) y (2): • Capa de aplicación: Dependiendo de la tarea, se pueden usar diferentes tipos de aplicaciones software, según las necesidades o requerimiento del uso que le demos a la red de sensores inalámbricos. • Capa de transporte: Asegura la fiabilidad y calidad de los datos en el proceso. Podemos ver una lista de los protocolos de comunicación de la capa de transporte en (2). • Capa de red: Enruta los datos que manda la capa de transporte a través de la red desde su fuente hasta el destino. • Capa de enlace de datos: Maneja el acceso al canal de datos para evitar errores. Dado que el medio es inalámbrico, es necesario un control del acceso al medio. • Capa física: Lleva a cabo las necesidades de métodos de modulación, recepción y envío. Además, los planos de energía, movilidad y de gestión de tareas coordinan las tareas y reducen el consumo de energía de la red. En la figura 2 se observa una representación del modelo de comunicación de las redes de sensores: 5 Figura 2. Arquitectura del protocolo de comunicación de una WSN (1). Otra de las áreas de investigación actual es la eficiencia de la red en el diseño, además de tener en cuenta el coste y los requerimientos de la tarea a realizar, para ello hay que optimizar el hardware y el software utilizado. Esto incluye la elección de los sensores en la parte de hardware y el software se refiere a la seguridad, robustez, tolerancia de fallos, organización, la red utilizada para el intercambio de información (2). Se dividen los problemas actuales de la tecnología en tres categorías que son (2): • La plataforma de las distintas capas y el sistema operativo. • La arquitectura del protocolo de comunicación. • Los servicios de red, aprovisionamiento y problemas de despliegue. También se analizan y comparan las posibles soluciones y las investigaciones actuales en cada categoría. En este trabajo se va a suponer que disponemos de una red de sensores inalámbricos de tipo desestructurado, ya que permiten una colocación más sencilla en el espacio y resulta conveniente para nuestro método de predicción mediante kriging. Origen y Análisis de los datos. Simulación de una red de sensores 12 12 Figura 6. Ubicación de los nodos de la red de sensores simulada. 4.2 Origen de los datos Una vez elegidos los puntos, el siguiente paso es buscar una base de datos disponible que se ajuste a las necesidades de nuestro problema. No hay datos disponibles para puntos específicos de la zona que se ha elegido, ya sean mediante estaciones meteorológicas o sensores colocados, por lo que se hace necesario usar datos obtenidos mediante imágenes de satélite y métodos, como el método Heliosat-2, que permite convertir observaciones de satélites en estimaciones de la radiación solar en el suelo (14). Se ha optado por usar el CAMS Radiation Service, que tiene disponible en el portal SoDa, este servicio se encuadra dentro del programa Copernicus, el programa de observación de la Tierra de la Unión Europea con el objetivo de tener datos sobre el estado de nuestro planeta y el medioambiente. Otros servicios que se han encontrado y que podrían ser útiles son, el PVGIS, servicio de la Unión Europea para obtener información sobre radiación y herramientas de comprobación de eficacia de paneles fotovoltaicos o el CM SAF, la agencia de monitorización climática de la unión europea, dispone de muchos productos de datos de variables climatológicas, incluyendo series históricas de datos y a tiempo real. No se han utilizado porque son más complejos de descargar los datos y su posterior utilización, pero proporcionan alternativas aceptables, en caso de que no estuviera disponible el servicio utilizado. Toda la información generada dentro del servicio CAMS Radiation es libre y gratuita para todos los usuarios, si se especifica al usar los datos lo indicado en el acuerdo de licencia (15). La guía de usuario del servicio nos repasa los métodos utilizados con anterioridad para el cálculo de la radiación solar. Este servicio utilizó previamente el método Heliosat-2 con el que se construyeron las bases de datos HelioClim-3 y SOLEMI que proporcionan datos de radiación solar, usando distintas implementaciones para el cálculo de la radiación solar, es decir, las dos bases de datos se basan en modelos de cálculo de la radiación solar distintos. También se avisa de que el método Heliosat-2 tiene limitaciones, como que no tiene en cuenta la aparición de nieve repentina. (16). En una versión posterior, actualizaron el método al Heliosat-4, que es el usado actualmente en el servicio CAMS. Se basa en variables como las propiedades de aerosoles, vapor de agua total en columna o el ozono. Es capaz de 13 estimar la radiación a nivel de suelo en todas las condiciones meteorológicas del cielo, siendo posible calcular otras variables como los componentes directo y difuso o la radiación en plano horizontal (16). Está disponible el informe de validación por trimestre sobre las estimaciones, comparando con datos de estaciones a nivel del suelo, y podemos leer las limitaciones existentes para validar los datos. Se ha leído parte de los informes y se concluye que existe un error en las medidas, tanto por sobreestimación como por subestimación y muestran distintos valores para indicar el grado de error cometido por estación, como el RMSE, por lo que concluimos que los datos no son del todo exactos. Figura 7. Interfaz de usuario del CAMS Radiation Service. Hay datos disponibles desde febrero de 2004 hasta datos de hace dos días desde la fecha en que lo miremos, con resolución temporal de 1 min, 15 min, 1 hora y 1 día. Se ha descargado un periodo de cuatro años, por lo que se ha optado por resolución de 15 min, ya que el tamaño de los archivos generados provocaría problemas de espacio y las lecturas de datos serían muchos más lentos. Si queremos un determinado punto, el servicio interpola para obtener el valor de la variable en el punto que le indiquemos, lo que equivaldría a disponer en la realidad de un sensor en ese punto. Para descargar los datos, hay que ir al portal e introducir las coordenadas del punto deseado e introducir el rango de tiempo a calcular, y seleccionar la resolución temporal y el formato del documento, que puede ser “.nc” o “.cv”, se han descargado en formato “.cv” ya que el programa Excel puede manejarlos sin mucha dificultad. Una vez descargado el fichero que se genera, si lo abrimos con Excel (puede ser necesario cambiar algunas opciones de la configuración del programa para leer los datos correctamente) podemos ver que tenemos en intervalos de 15 min los datos de las siguientes variables: • TOA: Radiación en plano horizontal por encima de la atmósfera. [Wh/m2]. • GHI: Radiación Global en plano Horizontal a nivel del suelo [Wh/m2]. • BHI: Beam irradiance o Radiación directa en plano Horizontal a nivel del suelo [Wh/m2]. • DHI: Radiación difusa en plano Horizontal a nivel del suelo [Wh/m2]. • BNI: Radiación directa en plano móvil perpendicular a la incidencia solar. [Wh/m2]. También están estas variables con cielo despejado (clear sky), que son los datos que se tendrían si no hubiera nubosidad sobre nuestro punto. Reliability: Última columna, proporción fiable de datos (0-1). Mostramos una captura de los datos descargados mediante el programa Excel: Origen y Análisis de los datos. Simulación de una red de sensores 14 14 Figura 8. Visualización de los datos descargados. “Generated using Copernicus Atmosphere Monitoring Service Information 2019” En las horas de noche, se observa que la radiación que tenemos es cero. También podemos ver el efecto de las nubes en las horas de día, en las que tendremos una caída del valor de la radiación global en la zona afectada. Se trabaja con la GHI, o radiación global en plano horizontal a nivel del suelo, corresponde a la suma de las radiaciones que llegan por radiación difusa y directa ya que los paneles más usados absorben el total de la radiación que llega, aunque existe una excepción, los sistemas de energía solares de concentración o CPV, que usan paneles curvos que concentran la energía solar que llega a las células fotovoltaicas que utilizan principalmente la radiación directa, es una tecnología en desarrollo, en el futuro podrían sobrepasar la eficiencia de la tecnología actual. Se han descargado 49 datos para predecir los valores de radiación solar en toda nuestra zona elegida, mediante kriging, para ello es necesario estudiar la variabilidad espacial entre los puntos de forma previa y tenemos 10 datos de validación, que usaremos para ver la diferencia entre la estimación y el valor “real” de la variable. 4.3 Análisis exploratorio de datos para cuatro periodos de tiempo Los datos se han leído mediante Matlab, y se han creado tablas con datos para que la lectura sea sencilla mediante R. Se ha obtenido la media de todos los valores mayores que cero, dentro de los siguientes periodos: a) Cada mes de 2015 b) 3 días de junio. c) Una hora de junio. d) Media de los valores de un mismo periodo de 15 min de todos los días de junio. El objetivo de obtener datos de periodos distintos es comparar las diferencias entre los modelos de los variogramas generados y las discrepancias que se producen en las predicciones. Se ha generado una nube de puntos, mediante la comparación del valor de la radiación de una pareja de puntos en cada instante, y se ha ajustado la recta de mayor proximidad a los puntos: 15 Figura 9. Nube de puntos. Los puntos van en azul y en rojo se representa la recta con mejor ajuste a todos los puntos. Los puntos 3 y 11 están muy cerca, mientras que el 16 y el 2 están más alejados, resulta obvio que cuanto mayor sea la distancia entre los nodos, mayor será la dispersión de los valores en un mismo instante. El análisis de los datos es indispensable antes de trabajar con ellos con el fin de comprender los patrones que siguen y ayudar a analizarlos mediante la obtención de sus características. Se dispone de librerías para realizar el análisis de forma sencilla mediante R, por lo que se van a utilizar. Se guardan los ficheros .txt, generados mediante Matlab con los datos en forma de tabla, con los datos en tres columnas separados mediante comas. Es mejor colocarlos dentro del directorio, no es necesario, pero facilita la ruta a los ficheros, que se leen mediante la función read.table de R, de la siguiente forma: datos <- read.table("Radiacion_media_mes.txt", header=T,sep = ",") Recordamos que en el anexo se puede leer el código completo. Si usamos la función summary(),se puede obtener información sobre los datos como la mediana, la media o los valores de los cuartiles. Por ejemplo, con summary(datos$junio) Min. 1st Qu. Median Mean 3rd Qu. Max. 131.4 132.9 133.5 133.5 134.3 135.5 Mediante la función Desc(),obtenemos una descripción más extensa de la variable, con valores para la desviación estándar, asimetría entre otros, además realiza una representación gráfica de los datos, que mostramos a continuación: Origen y Análisis de los datos. Simulación de una red de sensores 16 16 Figura 10. Gráfica generada mediante Desc del mes de junio. El método de kriging genera las mejores estimaciones posibles si nuestros datos siguen una distribución normal, aunque no es un requisito fundamental, nos interesa obtener el menor error posible en las estimaciones, por lo que cuanto más cerca nos encontremos de una distribución normal, mejores resultados. En (17), se realiza un estudio con los valores normalizados y no normalizar para aplicar kriging ordinario y no se encuentran diferencias substanciales. En el caso del mes de junio (que tenemos en la figura 10), los datos presentan una distribución normal, pero otros meses no siguen esta distribución. Para comprobarlo, aplicamos el test de Shapiro-Wilk, si el resultado de p es menor que 0.05, la distribución de los datos no es normal, si es mayor es posible que sea normal. Estos tipos de tests se ven afectados por el tamaño de la muestra, por lo que se suelen complementar con un gráfico q-q, que es una comprobación subjetiva, de manera gráfica, de la distribución que siguen los datos, si los puntos en la gráfica siguen una distribución normal, deberían formar una línea recta diagonal (18). Shapiro-Wilk normality test. data: datosorig$Junio W = 0.98903, p-value = 0.9257 A continuación, representamos el gráfico q-q para el mes de junio: 17 Figura 11. Gráfica q-q para el mes de junio. En el caso de que no fuera normal, se realiza una normalización de los datos mediante métodos como aplicar el logaritmo, raíces cuadradas o un escalado a los datos. También se ha utilizado una librería, bestNormalize (19), que permite aplicar una función con el mismo nombre, cuyo objetivo es realizar posibles transformaciones de normalización con los datos. Incluso utilizando este recurso, algunos datos no han pasado el test de normalidad y no se aprecian diferencias significativas en las gráficas q-q, por lo que hemos decidido trabajar con datos que no se ajustan totalmente a la distribución normal, vemos un ejemplo con los datos de los días de junio: Figura 12. Distribución de los valores de los datos de los días de junio Origen y Análisis de los datos. Simulación de una red de sensores 18 18 Figura 13. Gráfica q-q para los datos de días de junio Se puede observar que una parte de los datos de la izquierda no se encuentran sobre la recta, si observamos el histograma, vemos que existe una asimetría negativa, esto quiere decir que la distribución de los datos se alarga en el lado izquierdo, la causa de esto es posiblemente debido a factores meteorológicos como que hubiera estado más nublado por la zona de esos puntos durante el tiempo en el que se han tomado los datos. Por otro lado, el resto de los puntos se ajustan de manera adecuada a la recta, por lo que no tendremos muchos problemas con el estimador. La distancia entre los puntos se ha vuelto a calcular mediante dos fórmulas, Vincenty y Haversine, para comprobar la variación entre las distancias en los puntos, que resultan ser pequeñas. Tabla 1. Distancia entre los puntos según la fórmula. Fórmula Distancia máxima [m] Distancia media [m] Haversine 8780 3426 Vincenty 8772 3424 A continuación, representamos un histograma de los datos, en rojo marcamos la media de los datos. 19 Figura 14. Histograma del mes de junio Para visualizar gráficamente los datos, realizamos una representación espacial en forma de mapa, con un gradiente de colores para el valor de la radiación. Figura 15. Mapa con los valores de radiación de los nodos. 4.4 Reproyectar La reproyección consiste en cambiar los valores de coordenadas del conjunto de los datos de un sistema de coordenadas a otro. En el momento en el que se calcularon las posiciones de los puntos en el espacio, se asignaron valores de coordenadas según el sistema de coordenadas geográficas de latitud y longitud, pero para calcular las distancias de forma correcta con las funciones que vamos a utilizar para calcular el variograma, hay que reproyectar las coordenadas a otro sistema de coordenadas de referencia (CRS) que permita calcular distancias de forma Origen y Análisis de los datos. Simulación de una red de sensores 20 20 sencilla. Se utiliza el sistema de coordenadas universal transversal de Mercator (UTM por sus siglas en inglés), que se expresa únicamente en metros al nivel del mar. En el sistema UTM, se separa las latitudes de la tierra en 60 zonas, con 6º grados de anchura cada zona, también llamado huso y se proyecta cada zona mediante la proyección Mercator con una baja distorsión. Los valores de los ejes X e Y tienen un centro de coordenadas en cada zona, que se identifica con un número, en España pasan las siguientes zonas por el territorio: Figura 16. Zonas UTM de España (20). A nuestra localización, en Almería, le corresponde la zona 30 del hemisferio norte. Los sistemas de proyección se pueden identificar mediante un código EPSG, en este caso el 32630. Se puede comprobar en (21), usamos este código para transformar de forma sencilla. Para reproyectar los nodos se instala la librería “rgdal”, con funciones que facilitan las operaciones a realizar, que mostramos a continuación: coordinates(datos) <- ~lon+lat # necesita libreria sp #proj4string(datos) = "+proj=longlat +datum=WGS84" proj4string(datos) <- CRS("+init=epsg:4326") # sistema de coordenadas lat+long #+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs epgs:32630 UTM zone 30N datos=spTransform(datos, CRS("+init=epsg:32630")) De esta manera se obtiene el objeto espacial datos con las coordenadas pasadas al sistema UTM, en la zona 30N, cuyas distancias entre nodos se pueden calcular de forma sencilla con nuestras funciones. 21 5 VARIABILIDAD ESPACIAL Y ESTIMACIÓN DE LA RADIACIÓN SOLAR n este capítulo se realiza el análisis estructural de los datos en los cuatro periodos de tiempo elegidos, se observan las variaciones de los datos en las distintas direcciones del espacio y se representan los modelos teóricos de variograma que mejor se ajustan a nuestros datos en los distintos periodos, también se muestran los variogramas de los distintos meses del año. Finalmente, se calculan las predicciones de la radiación solar en nuestra zona, siendo representadas de forma gráfica y se compara los resultados de la predicción obtenidos mediante los distintos modelos generados por los datos de los distintos periodos. 5.1 Variabilidad espacial de la radiación solar El paquete gstat (22) contiene la mayor parte de las funciones de geoestadística utilizadas, un ejemplo sobre la puesta en práctica de estas se puede consultar en este documento (23). Como trabajo previo a la predicción de la radiación solar se ha llevado a cabo el análisis estructural de los datos. En dicho análisis se estudia la variabilidad espacial de la radiación solar en los cuatro periodos de tiempo en los que estamos trabajando. El único modelo que se ha normalizado ha sido el de la media de los instantes del mes de junio y se han comparado los modelos de variograma normalizados y sin normalizar. El variograma se calcula mediante la función variogram, incluida en el paquete gstat, antes de poder usar esta función hay que crear un objeto estadístico con la función gstat, cuando ejecutamos variogram obtenemos una variable que podemos dibujar con el comando plot. Figura 17. Variograma experimental con la muestra de datos del mes de junio. E Variabilidad espacial y estimación de la radiación solar 28 28 Se observa una región con menos radiación en la parte inferior izquierda de la imagen, esto es debido a la presencia de nubes por esa zona, el resto del mapa presente una radiación uniforme sobre los 255 Wh/m2, según la predicción realizada con el modelo del mes de junio, en la figura 29, podemos ver las diferencias entre las predicciones según el modelo de variograma utilizado con kriging. También es posible utilizar un parámetro obtenido de la función, una variable que nos indica la varianza de cada región calculada, si hacemos la raíz cuadrada para calcular la desviación típica, podemos mapear el error que pueden tener las predicciones de nuestro modelo de variograma. Figura 26. Desviación típica de la predicción de horas con el modelo de variograma de junio. Las zonas azules corresponden a las zonas más cercanas a los puntos de los que disponemos el dato, por eso el error cometido debe ser menor, conforme nos vamos alejando de nuestros datos, aumenta el error que cometemos al estimar. Para calcular el error cometido en las estimaciones de krige con nuestros propios datos, podemos utilizar la función krige.cv, que utiliza validación cruzada, de forma que dejando fuera del cálculo de krige un punto observado, predice el valor mediante krige ordinario en la localización del punto eliminado y lo compara con el observado, procediendo de esta forma con todos los puntos que tenemos como datos, y representamos el resultado: 29 Figura 27. Resultado de la validación cruzada con el modelo de junio del variograma. Los puntos verdes son los valores predichos frente al dato que tenemos, para comprobar la correlación se calcula el coeficiente R2. La función nos devuelve una variable que contiene el residual que es la diferencia entre observado y estimación para cada dato observado de la muestra, y los valores estimados, procedemos a calcular las siguientes medidas: 𝑀𝐴𝐸=1 𝑁∑| 𝑁 𝑖=1 𝑍󰆹(𝑠𝑖) − 𝑍(𝑠𝑖)| 𝑅𝑀𝑆𝐸= √∑[𝑍󰆹(𝑠𝑖) −𝑍(𝑠𝑖)]2 𝑁 𝑖=1 𝑁 Donde: N= Número de datos. 𝑍󰆹=Estimaciones. Z= Valores observados. Con los modelos que se han calculado para los periodos del mes de junio, los días, horas e instantes, es una buena idea utilizar la función krige.cv para estimar con los propios datos cada periodo (datos del mes y modelo del mes, datos de días y modelo de días, …) y se procede a comparar los resultados de aplicar estas fórmulas. Variabilidad espacial y estimación de la radiación solar 30 30 Tabla 2. Resultados de la validación cruzada de los datos. Periodo de datos R2 RMSE MAE Mes 0,884 0,318 0,255 Días 0,844 0,515 0,356 Hora 0,661 7,565 4,666 Media de instantes 0,349 0,758 0,577 Cuanto más próximo sea el valor de R2 a 1, más exacto resultan las estimaciones generadas mediante el modelo y a menor RMSE o MAE menores son los errores cometidos. Se observa que cuanto menor sea el periodo temporal en el que se realiza la media para obtener los datos (por ejemplo, la media en el tiempo del mes tiene más que la del día) para calcular el modelo del variograma, peor parece ser la predicción de nuestros datos, probablemente debido a que las nubes afectan a la radiación solar y cuantos menos datos se tengan en cuenta mayor sea el efecto perturbador de las nubes sobre el cálculo de los valores, que recordemos que se ha realizado la media de los valores mayores que cero dentro del periodo, por lo que si tenemos muchos datos las nubes que puedan haber afectado algunos días no tendrán mucho efecto sobre el total. Por esta razón, el valor de RMSE y MAE es mayor en la hora de junio, donde podemos encontrar nubes por una parte de nuestra zona como veremos mejor después. La media de instantes de junio se ha calculado con datos normalizados por lo que el RMSE sale pequeño debido a que son valores normalizados, pero tiene el menor R2 de todos. Los cálculos de kriging se realizan con datos sin normalizar debido a que si se transforman los datos se pierde los valores reales de las variables y es más complicado realizar la transformación a los valores originales. Se procede a comparar la diferencia entre usar datos normalizados y no normalizados en la predicción representando el mapa de la predicción en la figura 25. Figura 28. Comparación de las estimaciones con datos normalizados y sin normalizar. 31 Los datos normalizados están a la izquierda, aunque los resultados puedan ser mejores por ajuste del R2 si calculamos el krige.cv (explicado a continuación), perdemos el significado físico de los datos ya que en vez de conocer la radiación conocemos un valor normalizado. 5.2.2 Comparación de predicciones con distintos modelos de variograma. Primero vamos a representar las predicciones que realizan los distintos modelos, para comparar las diferencias existentes. Figura 29. Comparación de las predicciones de la hora de los distintos modelos de variograma. Se observan claras diferencias entre las predicciones realizadas por los distintos modelos de los variogramas, pero no se puede estar seguro sobre que predicción se ajusta mejor a los datos reales de radiación, para resolver esto, se comprueban las predicciones con valores reales, tenemos 10 coordenadas con datos que no se han usado para predecir, ahora los utilizamos para comprobar si las predicciones se ajustan a la realidad. Se va a representar los puntos que vamos a añadir para comprobar los resultados de la predicción para poder localizarlos de forma gráfica en el mapa. Variabilidad espacial y estimación de la radiación solar 32 32 Figura 30. Puntos de validación sobre la predicción. A continuación, se guardan las predicciones en un archivo tipo .csv, para ser utilizados para calcular el error en Matlab, en ese archivo se tienen los valores de las coordenadas y de la predicción para cada punto. En Matlab se utilizan las coordenadas en el sistema latitud y longitud por lo que antes hemos de reproyectar a ese sistema de coordenadas (WSG84). Entonces, en Matlab se obtiene el punto de la predicción más cercano a los 10 que se usan para comprobar con datos reales (recordar que tenemos una cuadrícula de puntos con predicciones) y se calcula el error absoluto de la predicción, restando al valor real el estimado en valor absoluto y el error relativo se calcula igual, pero dividiendo entre el valor real del dato y se expresa en porcentaje multiplicando por 100. Se ha calculado la media de los errores de los 10 puntos de validación y se muestra en la siguiente tabla para la estimación de las horas. Tabla 3. Media de los errores de la predicción de la hora (c) comparando con datos reales. Modelo Error absoluto medio Error relativo (%) Mes 3,36 1,51 Días 6,50 2,83 Horas 4,20 1,83 Instantes 7,82 3,49 33 Se muestran los errores de las predicciones de los distintos modelos de variograma calculados realizadas con los propios datos con los que se han calculado los modelos, por ejemplo, los datos del mes que se han utilizado para generar el modelo del mes hacemos kriging usando esos datos con los modelos b), c) y d) de variograma calculados. Esto se ha realizado con todos los datos de los modelos de variograma del mes de junio, b), c) y d) y los mostramos en las siguientes tablas con los errores cometidos. Tabla 4. Media de los errores de la predicción de los datos del mes (a) comparando con datos reales. Modelo Error absoluto medio Error relativo (%) Mes 0,287 0,21 Días 0,404 0,30 Horas 0,424 0,32 Instantes 0,445 0,33 Tabla 5. Media de los errores de la predicción de los datos de días (b) comparando con datos reales. Modelo Error absoluto medio Error relativo (%) Mes 0,344 0,26 Días 0,465 0,35 Horas 0,369 0,28 Instantes 0,844 0,64 Tabla 6. Media de los errores de la predicción de los datos de instantes (d) comparando con datos reales. Modelo Error absoluto medio Error relativo (%) Mes 1,09 0,58 Días 1,8 0,95 Horas 1,98 1,05 Instantes 1,01 0,54 Se puede comprobar que, aunque el modelo coincida exactamente con los datos, el modelo del mes obtiene las mejores predicciones excepto en la media de los instantes del mes, aunque la diferencia no sea grande entre los errores cometidos. Variabilidad espacial y estimación de la radiación solar 34 34 También se observa que cuanto menor sea la cantidad de datos observados utilizados para generar la media que se usa como valor en los modelos, mayores son los errores cometidos, esto se debe a la influencia que tienen las nubes que existan en un determinado momento sobre cada punto, cuanto menor sea el número de datos, mayor efecto tienen la caída de la radiación provocadas por un cielo nublado sobre su valor. En algunos estudios, observan que los modelos de las estaciones del año para la radiación solar, funcionan igual de bien que los modelos de días y horas, sería posible reducir la necesidad de cálculos al ser posible usar el modelo de variograma del mes para estimar los valores, por lo que no sería necesario calcular un modelo teórico de variograma para los propios datos antes de estimar, ya que el modelo del mes produce buenas estimaciones para datos horarios (8). 35 6 CONCLUSIONES l objetivo de este trabajo era lograr obtener estimaciones de los valores de la radiación para nuestra planta solar hipotética en el desierto de Almería sin sensores colocados en el terreno, se ha investigado para obtener los datos que se ajusten a una red de sensores distribuidos irregularmente en nuestra planta y se han utilizado datos obtenidos mediante satélites de la Unión Europea, disponibles de forma gratuita en la red. Entonces se ha utilizado la geoestadística, haciendo un análisis estructural de los datos para obtener la variabilidad espacial de la radiación solar por distintos periodos de tiempo y los tratamientos necesarios para, posteriormente, aplicar las técnicas de interpolación de kriging para obtener las predicciones de los valores de la radiación solar, por lo que se podría aplicar este proceso para mejorar la eficiencia de nuestros colectores solares ya que se conocen las zonas con la mayor radiación solar. Al comparar las predicciones con otros valores reales vemos que los errores no son demasiado grandes, por lo que podemos decir que kriging es un método válido para predecir los valores de radiación solar en una zona geográfica de tamaño extenso. Nuestros datos no se han obtenido de forma experimental mediante sensores sino mediante datos de satélite por lo que no es exactamente una aplicación a una situación real, pero puede servir de base a desarrollar para posteriores trabajos con datos reales. Los datos obtenidos mediante satélite no son exactos, pero suponen un ahorro en los costes de los sensores piranómetros, su instalación y su mantenimiento, por lo que podría considerarse su uso en aplicaciones no críticas. Se ha realizado una búsqueda de las bases de datos disponibles de forma gratuita en Internet, se ha hecho una pequeña recopilación con las más accesibles y se ha explicado cómo funcionan. Líneas de investigación futuras. • Utilizar datos obtenidos de forma real mediante sensores, ya sea por redes de sensores u otros métodos como drones y comparar la diferencia entre los valores obtenidos mediante satélites y si las predicciones varían de forma significativa. • Comparar distintos métodos de kriging y comparar las desviaciones de las predicciones. • Comparar con distintos métodos de interpolación, por ejemplo, con la distancia inversa ponderada (IDW por sus siglas en inglés), red irregular triangulada o splines, y estudiar cómi varían los resultados de la predicción de la radiación solar. • Aplicar las predicciones de los valores al control de los paneles solares para estudiar el posible aumento de la eficacia en la captación de energía solar. • Estudiar el número de nodos de red óptimos dentro de una región y su disposición mediante kriging con el fin de reducir el coste en el número de sensores sin perder información de la radiación. E Anexo a: Código realizado 36 36 7 ANEXO A: CÓDIGO REALIZADO 7.1 Código realizado en R. library(gstat) library(sp) library(mapview) library(DescTools) library(ggplot2) library(rgeos) library(automap) #autofitvariogram library(rgdal) #reproyectar library(latticeExtra) #AÑADIR cosas a spplot library(magrittr, warn.conflicts = FALSE) #pipe operator library(MASS) library(bestNormalize) library(rngtools) rm(list = ls()) #borrar workspace como clear all en matlab # opciones para ggplot theme_set(theme_bw()) theme_update(legend.position='right') #Cargamos todos los datos. #datosmes<- read.table("Radiacion_media_mes.txt", header=T,sep = ",") #datosdias<- read.table("Radiacion_media_dias.txt", header=T,sep = ",") #datoshoras <- read.table("Radiacion_media_hora.txt", header=T,sep = ",") #datosinst<- read.table("Rad_media_instantes_junio.txt", header=T,sep = ",") #Elegimos los datos a cargar. datos_orig <- read.table("Radiacion_media_dias.txt", header=T,sep = ",") 37 datos <- datos_orig coordinates(datos) <- ~lon+lat # necesita libreria sp #proj4string(datos) = "+proj=longlat +datum=WGS84" proj4string(datos) <- CRS("+init=epsg:4326") # sistema de coordenadas normal -> lat+long #necesario reproyectar los datos. #+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs epgs:32630 UTM zone 30N #EPSG:25830 PARA EUROPA. PROYECCION UTM HUSO 30 N datos=spTransform(datos, CRS("+init=epsg:32630")) #dev.off() shapiro.test(datos_orig$mediadia) #>0.05 gaussiana summary(datos_orig) Desc(datos_orig$mediadia) ##Intentos de normalizar ## pru1<- log(datos_orig$mediadia) Desc(pru1) shapiro.test(pru1) pru2<- scale(datos_orig$mediadia) #Desc(pru2) pru4<- datos_orig$mediadia normal=(pru4-min(pru4))/(max(pru4)-min(pru4)) Desc(normal) shapiro.test(normal) Anexo a: Código realizado 44 44 ( p2=spplot(resultado_dias,"var1.des", asp=1, col.regions=cm.colors(80), main='Desviacion', scales=list(draw=TRUE),xlab="X (m)", ylab="Y (m)", sp.layout=comprobacion.layer )) • Análisis exploratorio de los datos. #Estudio estructural de los datos. library(gstat) library(sp) library(mapview) library(DescTools) library(ggplot2) library(rgeos) library(automap) #autofitvariogram library(rgdal) #reproyectar library(geosphere) #distHaversine # opciones para ggplot theme_set(theme_bw()) theme_update(legend.position='right') # datos #list.files() #datos <- read.table("~/R/geostadistica proyecto prueba/proyecto/BroomsBarn.txt", sep="\t") datos <- read.table("Radiacion_media_mes.txt", header=T,sep = ",") datosorig=datos #copia seg coordinates(datos) <- ~lon+lat # necesita libreria sp #proj4string(datos) = "+proj=longlat +datum=WGS84" proj4string(datos) <- CRS("+init=epsg:4326") # sistema de coordenadas normal -> lat+long 45 #exploracion de datos summary(datosorig) Desc(datosorig$Junio) #Test de normalidad shapiro.test(datosorig$Enero) #Enero no tiene distribución normal. shapiro.test(datosorig$Junio) qqnorm(datos$Enero, pch = 1) qqline(datos$Enero, col = "steelblue", lwd = 2) apply(datosorig,2,var) head(datosorig) #quantile(dist(datosorig[,1:2])) Para ver las distancias, pero es para planos. ma <- distm(datosorig[,1:2]) #Calcula la matriz de distancias quantile(distm(datosorig[,1:2],datosorig[,1:2], fun=distVincentyEllipsoid)) #Distancia máxima en metros, media,.. quantile(distm(datosorig[,1:2],datosorig[,1:2], fun=distHaversine)) #Varía un poco según el método usado #pointDistance(datosorig,lonlat=TRUE) ggplot(datosorig,aes(Enero))+ geom_histogram(aes(y=..density..),bins=10,col=1,fill=4,alpha=.5)+ geom_vline(xintercept = mean(datosorig$Enero),col=2)+ geom_density(col=4)+ labs(y='Densidad') # preprocesamiento para analisis espacial #calculamos los grids, la distancia dint= signif(max(c(diff(range(datos$lon))/length(datos$lon),diff(range(datos$lat))/len gth(datos$lat))),3) # generar secuencia de las x e y xint= seq(min(datos$lon),max(datos$lon),dint) #coord. x del grid para la interpolacion yint= seq(min(datos$lat),max(datos$lat),dint) datosint= expand.grid(x=xint,y=yint) Anexo a: Código realizado 46 46 gridded(datosint)=~x+y # grid como objeto espacial mapview(datos, burst= T, hide=T) mapview(datos, zcol='Junio', legend=T) # contorno que contiene los datos, edge de la figura #No funciona porque necesica planar coordinates #q=min(c(diff(range(datos$lon)),diff(range(datos$lat)))) #outline= gBuffer(datos, byid=FALSE, id=NULL, width=dint*q, # joinStyle = "ROUND", quadsegs=10) #outline= gBuffer(outline, byid=FALSE, id=NULL, width=-dint*q, # joinStyle = "ROUND", quadsegs=10) #plot(outline) 47 7.2 Código en Matlab Generación de los puntos aleatoriamente: generar_puntos.m % Generar los puntos del mapa clear, close NPUNTO=10; %Numero estaciones AREA=50; %Se ha hecho con 50 km2 TAM_SUPERF=sqrt(AREA); %AREA CUADRADA, LADO EN KM %N números aleatorios en el intervalo (a,b) con la fórmula r=a+(b-a).*rand(N,1) puntos=TAM_SUPERF.*rand(NPUNTO,2); for i=1:NPUNTO hold on plot(puntos(i,1),puntos(i,2),'ro-') end axis([ 0 TAM_SUPERF 0 TAM_SUPERF]) title('Superficie') xlabel('X') ylabel('Y') % Guardar una generacion • Cálculo de coordenadas, plot de los puntos y cálculo de distancia entre los puntos. %Aqui tenemos los datos de los puntos en un plano, los pasamos a coordenadas y se mapean finalmente calculamos la distancia. close NPUNTO=50; %Numero estaciones AREA=50; %Elegir area TAM_SUPERF=sqrt(AREA); %AREA CUADRADA, LADO EN KM lat1=37.105; lat2=37.041; lon1=-2.488; lon2=-2.5675; Anexo a: Código realizado 48 48 diflat=lat1-lat2; diflon=lon1-lon2; %Puntos elegidos de forma aleatoria y guardados, se cargan. load estaciones1 %declaracion de vector de coord coordenadas_est=ones(NPUNTO,2); %Grafica de las estaciones en un plano figure() for i=2:NPUNTO hold on plot(estaciones(i,1),estaciones(i,2),'ro-') end axis([ 0 TAM_SUPERF 0 TAM_SUPERF]) title('Superficie') xlabel('X') ylabel('Y') %Pasamos a coordenadas con una relación lineal. for i=1:NPUNTO coordenadas_est(i,1)=estaciones(i,1)*diflon/TAM_SUPERF+lon2; coordenadas_est(i,2)=estaciones(i,2)*diflat/TAM_SUPERF+lat2; end filename = 'testdata.xlsx'; xlswrite(filename,coordenadas_est,1,'C2') %El primer punto no lo tengo como estación, por eso lo salto. figure() for i=2:NPUNTO hold on plot(coordenadas_est(i,1),coordenadas_est(i,2),'go-') end %axis([ lon2 lon1 lat2 lat1]) 49 title('Superficie') xlabel('X') ylabel('Y') %calculo de la matriz de distancias. %Las distancias se van a poner en una matriz de 49*49 simétrica, donde el %componente i,j es la distancia de i hasta j. %Para calcular las distancias usar haversine (o distance) coor=[coordenadas_est(:,2),coordenadas_est(:,1)]; %coor (:,1) contiene latitudes, (:,2) longitudes distancia=ones(49); %49*49 for i=2:50 for j=2:50 distancia(i-1,j-1)=haversine(coor(i,1),coor(i,2),coor(j,1),coor(j,2)); end end Anexo a: Código realizado 50 50 • Función Haversine para cálculo de distancias en la Tierra con coordenadas. function dist= haversine(lat1,lon1,lat2,lon2) %Calculo de la distancia usando haversine R= earthRadius('km'); % 6371 (el radio es aprox) dlat=deg2rad(lat2-lat1); dlon=deg2rad(lon2-lon1); lat1=deg2rad(lat1); lat2=deg2rad(lat2); a = (sin(dlat./2)).^2 + cos(lat1).*cos(lat2).*(sin(dlon./2)).^2; c = 2.*asin(sqrt(a)); dist=c*R; end 51 • Lectura de los datos de los ficheros, cálculos y transformaciones de los datos. %Varias versiones de este programa, esta calcula la media de los dias %14 a 17 de junio. %run puntoelegido.m %clear all close all sheet=1; %Elegir el rango de los datos con el que queramos calcular la media. Rango= 'G34411:G34794'; % seleccionar carpeta de forma interactiva (datos buenos o datos comprobar) ubicacion = uigetdir; %Guardamos los nombres como vector D=dir([ubicacion, '\*.csv']); %Extraer los nombres nombre_archivos={D(:).name} datos = cell(length(D),1); for i = length(D):-1:1 % Crear el nombre completo y parcial nombrecompleto = [ubicacion filesep D(i).name]; %filesep en windows es '\' % Leer datos datos{i} = xlsread(nombrecompleto,sheet,Rango); end %Pasar a matriz [L,N]=size(datos); radiacion_globalmediainstantes=[]; M = max(cellfun(@numel, datos)); % maxEl=M ;%max number of elements of a cell element in your cell % C=cell2mat(cellfun(@(x) [cell2mat(x) zeros(1,maxEl-numel(x))], Anexo a: Código realizado 52 52 datos,'UniformOutput',0)) % for j=1:L % datos(L)=[datos(L) %Nota importante: Para hacer la matriz tiene que tener el mismo tamaño for i=1:L dif=M-numel(datos{i}); if dif>0 radiacion_globalmediainstantes(:,i)=([cell2mat(datos(i)) ; zeros(1,dif)]); %todos tengan el mismo tamaño else radiacion_globalmediainstantes(:,i)=([cell2mat(datos(i))]); end end aux2=radiacion_globalmediainstantes'; %Quiero la matriz de esta forma. %Voy a eliminar los elementos cero para sumar la media sin ceros. [tamf,tamc]=size(aux2); cont=0; total=0; media_dia=zeros(tamf,1) % dias_mes=[31,28,31,30,31,30,31,31,30,31,30,31]; % datos_dia=4*24; %4 datos cada hora (cada 15 min). for i=1:tamf for j=1:tamc if aux2(i,j)>0 total=total+aux2(i,j); cont=cont+1; 53 end end media_dia(i)=total/cont; cont=0; total=0; end media_dias_datos=media_dia; datos2=media_dias_datos; %Guardar datos de validación. %save vali_dias datos2 %Ejecutar antes de esta parte puntoelegido.m %run puntoelegido.m %Tabla para usar en R. tabla_media=table(coordenadas_est(2:50,1), coordenadas_est(2:NPUNTO,2), datos2,... 'VariableNames',{'lon','lat','mediadia',}) writetable(tabla_media,'Radiacion_media_dias') %Vamos a pintar todos los puntos para estudiar los datos. % run nube_puntos.m • Nube de puntos: %Nube de puntos. Comparamos los valores en el mismo t entre estaciones. %HAY QUE EJECUTAR Media_enero.m antes! close all Referencias 60 60 21. Consejería de agricultura, ganadería, pesca y desarrollo sostenible. [Online] http://www.juntadeandalucia.es/medioambiente/site/rediam/menuitem.04dc44281e5d53cf8ca78ca73152 5ea0/?vgnextoid=2a412abcb86a2210VgnVCM1000001325e50aRCRD&lr=lang_es. 22. Multivariable geostatistics in S: the gstat package. Pebesma, E.J. 7, s.l. : Computers & Geosciences, 2004, Vol. 30. 23. Pebesma, Edzer. [Online] 16 Mayo 2019. https://cran.rproject.org/web/packages/gstat/vignettes/gstat.pdf. 61