xander_bakker

Generar un corredor de ruta óptima según el método EPRI

Blog Post created by xander_bakker on Sep 28, 2015

Uno de los ejemplos que se mostró durante la presentación de Python en ArcGIS Pro (CCU2015) fue la generación de un corredor según el método de EPRI.

En este documento detallamos cual fue el proceso.

 

El primer paso es generar un raster donde cada celda indique cual es el costo que genere para pasar la ruta por esta celda. Esto se hizo por 3 componentes, donde cada componente, cada capa en cada componente y cada elemento de las capas fueron asignado una ponderación y un porcentaje de influencia (esto fue demostrado por LIA PATERNINA de EPM en su ponencia "Modelo Predictivo Espacial de Ruta Óptima en Proyectos de Líneas de Transmisión Eléctrica"):

Con base es este raster se elaboren los siguientes pasos:

En primera instancia se definieron los accesos a los datos de entrada (el raster con los costos y los dos featureclasses con el punto de inicio y el punto final:

import os ras_cost = r"C:\CCU2015\2 Analisis Raster\2 Analisis Raster.gdb\cmp_total" fc_inicio = r"C:\CCU2015\2 Analisis Raster\2 Analisis Raster.gdb\punto_inicio" fc_final = r"C:\CCU2015\2 Analisis Raster\2 Analisis Raster.gdb\punto_fin"

 

Luego se definió el porcentaje de área a usar y el featureclass de salida:

porcentaje = 0.10 fc_salida = r"C:\CCU2015\2 Analisis Raster\2 Analisis Raster.gdb\MacroCorredor_v01"

 

Con base en el nombre del featureclass de salida, se determinó el espacio de trabajo (workspace) y se configuró el ambiente de geoprocesamiento:

ws, fc_nombre = os.path.split(fc_salida) arcpy.env.workspace = ws arcpy.env.overwriteOutput = True

 

Luego se utilizó el Cost Distance para calcular los costos acumulados desde el punto de inicio y también desde el punto final y se sumaron estos rasters:

costdist1 = arcpy.sa.CostDistance(fc_inicio, ras_cost) costdist2 = arcpy.sa.CostDistance(fc_final, ras_cost) costdist_suma = costdist1 + costdist2

 

Para poder sacar un porcentaje (10%) de las celdas, es necesario poder tener acceso a la tabla de atributos del raster con la suma de los costos acumulados. Dado que este raster es de tipo punto flotante, no tiene una tabla de atributos. Es necesario convertirlo a un raster entero, pero antes de esto redistribuir los valores a un rango fijo (en este caso 0 - 1000). Esto podemos hacer utilizando las propiedades mínimo y máximo del objeto raster:

cd_rango1000fl = (costdist_suma - costdist_suma.minimum) / (costdist_suma.maximum - costdist_suma.minimum) * 1000 cd_rango1000 = arcpy.sa.Int(cd_rango1000fl)

 

Es necesario hacer sumar el conteo de las celdas, empezando desde el valor de pixel inferior (más apto) hasta llegar a la cantidad de celda que corresponde al porcentaje configurado (en este caso 10%). Para este propósito generamos un diccionario con el valor de celda como clave y el conteo de celdas como valor a asignar (hay que utilizar el "viejo SearchCursor", ya que el SearchCursor del módulo da (data access) no puede ser utilizado con las tablas de atributos de los raster) :

fld_val = "Value" fld_cnt = "Count" dct = {r.getValue(fld_val): r.getValue(fld_cnt) for r in arcpy.SearchCursor(cd_rango1000)}

 

Para sacar el valor de pixel acumulado que corresponde al porcentaje que buscamos, se puede ejecutar el siguiente código:

total = sum(dct.values()) frac = total * porcentaje val, suma = 0, 0 for v, cnt in sorted(dct.items()):     suma += cnt     if suma <= frac:         val = v

 

En este caso definimos el total de celdas haciendo la suma de los valores (conteos) en nuestro diccionario. La cantidad de celdas que estamos buscando corresponde al porcentaje (10%) multiplicado al total de celdas.

Luego se asigne valores 0 a las variables val y suma y se comience con el bucle a través del diccionario. En el bucle se suma los conteos de las celdas y mientras esta suma sea inferior o igual a la fracción se asigna el valor de la celda a la variable "val". Al final del bucle, la variable "val" contiene el valor de celda que buscamos.

Este valor se usa en una "superposición" condicional, asignando el valor 1 a las celda que tienen un valor inferior o igual a la variable "val" y el resto de celdas obtendrá el valor NoData.

macro_corr = arcpy.sa.Con(cd_rango1000 <= val, 1)

 

Al final este raster se convierte a polígonos:

arcpy.RasterToPolygon_conversion(macro_corr, fc_salida)

Outcomes