La semana pasada recibí la entrega del proyecto de Fin de Master GIS con Python de dos alumnos de la Agencia de Medio Ambiente y Agua de la Junta de Andalucía. Han montado un proyecto muy interesante que cubre una de sus necesidades del día a día. La gestión de riesgo del incendios forestales es parte del núcleo de su actividad principal y la lectura e interpretación de datos meteorológicos es, como se puede entender, una rama fundamental para su desempeño diario. Pues bien, nos han planteado una herramienta que lee de forma automática los ficheros de datos meteorológicos, los une con los puntos de estaciones de medida y hace una análisis de interpolación con uso incluso de técnicas de Machine Learning, como al Cluster and Outliers Analysis.
La herramienta se implementa sobre un script de Python en ArcPy. Y para usarla se puede:
– Instalar un Addin
– Lanzarla directamente desde la ToolBox desarrollada a tal efecto
La Toolbox y los archivos que utiliza se pueden descargar en esta url y más abajo os dejamos el código Python completo.
Para utilizar la herramienta debemos lanzar la Toolbox o el Addin y nos va a pedir los siguientes parámetros.
– Máscara: sería el raster que sirve para definir la máscara, tamaño de pixel del raster de salida y snap a ajustas. Adjuntamos un ráster de referencia llamado malla100.
– Workspace: para definir el entrono de trabajo
– Estaciones: Fichero .txt con el identificador y coordenadas de la red de estaciones
– Valores DC: Fichero .txt con los valores DC de un día determinado. El DC es un índice de sequía muy empleado para definir el riesgo de incendios forestales. Su cálculo lo tenemos implementado en un desarrollo Visual Basic, donde introduciendo datos de la red de estaciones se calcula el índice para un día determinado, dando como resultado el archivo .txt a introducir en esta casilla.
– Leyenda: opción para aplicar al mapa de salida una leyenda definida. Adjuntamos una leyenda de referencia llamada Legend.lyr
– Borrar shapes intermedio: opción para mantener o borrar los .shp intermedios que genera la herramienta.
Desde el punto de vista de la programación es interesante, por ejemplo:
-La lectura del fichero de texto
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
# crea el fichero de puntos y sus campos outputFC = arcpy.CreateFeatureclass_management( workspace, nombre_fich, "POINT", plantilla) lstESTs = f1.readlines() lstDCs = f2.readlines() cur = arcpy.InsertCursor(outputFC) cntr = 1 for estDC in lstDCs: if 'N' in estDC: continue estacion = lstESTs[cntr] vals = estacion.split(";") X = float(vals[1]) Y = float(vals[2]) vals = estDC.split(";") N = int( vals[0]) iDC = float(vals[1]) pnt = arcpy.Point(X, Y) feat = cur.newRow() feat.shape=pnt feat.setValue("N", N) feat.setValue(campoDC, iDC) cur.insertRow(feat) arcpy.AddMessage("Record number " + str(cntr) + " written to feature class") cntr = cntr + 1 |
– La generación de Outlayers con el módulo de Análisis espacial
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
try: # Check out the ArcGIS Spatial Analyst extension license arcpy.CheckOutExtension("Spatial") # Detección de valores atípicos en el shape de entrada mediante los siguientes métodos concatenados arcpy.ClustersOutliers_stats(OF1, campoDC, OF2,"INVERSE_DISTANCE","EUCLIDEAN_DISTANCE", "NONE") arcpy.Select_analysis(OF2, OF3, "COType <> 'HL'") arcpy.HotSpots_stats(OF3, campoDC, OF4, "INVERSE_DISTANCE", "EUCLIDEAN_DISTANCE", "NONE") arcpy.Select_analysis(OF4, OF5, "Gi_Bin <> 3") # Execute IDW power = "3" Search_radius = "VARIABLE 12" arcpy.gp.Idw_sa(OF5, "DC", nombre_raster, CellSZ, power, Search_radius, "") # Carga ráster resultado al TC mxd = arcpy.mapping.MapDocument("current") df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0] miLayer = arcpy.mapping.Layer(nombre_raster) |
Os animamos a utilizar esta interesante herramienta y a compartir con nosotros vuestras impresiones.
Además, el script completo queda como sigue.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 |
#Nombre: cencina y jvenegas #Fecha:08/04/2016 #Objetivo: Importar datos DC y estaciones #Python # -*- coding: cp1252 -*- # Import system modules import os import arcpy from arcpy import env from arcpy.sa import * import pythonaddins campoDC = "DC" temporal1 = arcpy.env.overwriteOutput arcpy.env.overwriteOutput = True try: ambito = arcpy.GetParameterAsText(0) workspace = arcpy.GetParameterAsText(1) entorno = workspace + '/' + ambito m_Raster=arcpy.Raster(ambito) CellSZ=m_Raster.meanCellHeight arcpy.env.snapRaster = ambito arcpy.env.mask = ambito plantilla = workspace + '/SalidaDC_Plantilla' # abre el fichero de estaciones f1 = open(arcpy.GetParameterAsText(2)) # abre el fichero de DC f2 = open(arcpy.GetParameterAsText(3)) # saca la fecha DC del archivo aux1 = f2.name.split(".") aux2 = aux1[0].split("_") aux3 = aux2[2].split("-") sufijo = aux3[2] + aux3[1] + aux3[0] # el nombre de la clase de elementos de salida salidaFC = "DC" # el nombre de la leyenda del ráster de salida leyenda = arcpy.GetParameterAsText(4) # borrar los ficheros intermedios borrar = arcpy.GetParameterAsText(5) # el nombre del fichero shape de salida nombre_fich = salidaFC + sufijo # el nombre del fichero ráster de salida nombre_raster = workspace + '/r' + salidaFC + sufijo # crea el fichero de puntos y sus campos outputFC = arcpy.CreateFeatureclass_management( workspace, nombre_fich, "POINT", plantilla) lstESTs = f1.readlines() lstDCs = f2.readlines() cur = arcpy.InsertCursor(outputFC) cntr = 1 for estDC in lstDCs: if 'N' in estDC: continue estacion = lstESTs[cntr] vals = estacion.split(";") X = float(vals[1]) Y = float(vals[2]) vals = estDC.split(";") N = int( vals[0]) iDC = float(vals[1]) pnt = arcpy.Point(X, Y) feat = cur.newRow() feat.shape=pnt feat.setValue("N", N) feat.setValue(campoDC, iDC) cur.insertRow(feat) arcpy.AddMessage("Record number " + str(cntr) + " written to feature class") cntr = cntr + 1 except Exception as e: pythonaddins.MessageBox(e,"Errores generados") finally: del cur f1.close() f2.close() #pythonaddins.MessageBox(entorno,"Prueba") OF1 = workspace + '/' + nombre_fich OF2 = workspace + '/' + nombre_fich + 'cluster' OF3 = workspace + '/' + nombre_fich + 'select1' OF4 = workspace + '/' + nombre_fich + 'hotspot' OF5 = workspace + '/' + nombre_fich + 'select2' try: # Check out the ArcGIS Spatial Analyst extension license arcpy.CheckOutExtension("Spatial") # Detección de valores atípicos en el shape de entrada mediante los siguientes métodos concatenados arcpy.ClustersOutliers_stats(OF1, campoDC, OF2,"INVERSE_DISTANCE","EUCLIDEAN_DISTANCE", "NONE") arcpy.Select_analysis(OF2, OF3, "COType <> 'HL'") arcpy.HotSpots_stats(OF3, campoDC, OF4, "INVERSE_DISTANCE", "EUCLIDEAN_DISTANCE", "NONE") arcpy.Select_analysis(OF4, OF5, "Gi_Bin <> 3") # Execute IDW power = "3" Search_radius = "VARIABLE 12" arcpy.gp.Idw_sa(OF5, "DC", nombre_raster, CellSZ, power, Search_radius, "") # Carga ráster resultado al TC mxd = arcpy.mapping.MapDocument("current") df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0] miLayer = arcpy.mapping.Layer(nombre_raster) # Aplica la Leyenda if leyenda <> "": arcpy.ApplySymbologyFromLayer_management(miLayer,leyenda) # Refresca la vista arcpy.mapping.AddLayer(df, miLayer) arcpy.RefreshActiveView() arcpy.RefreshTOC() del mxd, miLayer # Shapes intermedios generados: if borrar == 'true': arcpy.Delete_management(OF1) arcpy.Delete_management(OF2) arcpy.Delete_management(OF3) arcpy.Delete_management(OF4) arcpy.Delete_management(OF5) arcpy.RefreshCatalog(workspace) except Exception as e: pythonaddins.MessageBox(e,"Errores generados") finally: arcpy.env.overwriteOutput = temporal1 |
Leave A Comment