Generación de Paneles Solares en techos Bogotá

Document created by xander_bakker on May 31, 2017Last modified by jesarmientoesri-co-esridist on Jun 20, 2017
Version 5Show Document
  • View in full screen mode

Durante el evento Transformación Digital e Inteligencia Operacional en las empresas de Servicios Públicos en la presentación de Jennyfer Sarmiento y Xander Bakker (ver presentación: "Evento Esri Colombia - Soluciones sobre ArcGIS "),  se mostró un ejemplo de generar paneles solares sobre techos en Bogotá.

 

 

En este documento explicaremos el proceso que aplicamos para obtener los resultados.

 

 

Datos de entrada

Partimos de los siguientes datos:

  • Imágenes capturados por un dron del sector de Chapinero de Bogotá
  • Polígonos de edificios del mismo sector

 

Primer paso: Procesar las imágenes con Drone2Map

La aplicación Drone2Map | ArcGIS, nos permite planear y procesar imágenes capturados por un dron. Algunos de productos que se pueden generar son:

  • Ortofoto mosaico
  • Modelo Digital de Superficie
  • Modelo Digital de Terreno (a partir de la versión 1.1)
  • Nubes de puntos 3D (formato LAS, PLY, XYZ)
  • Curvas de nivel
  • Indices de vegetación en caso de contar con imágenes multiespectrales
  • 3D textured mesh

 

 

Para este ejercicio solamente se usaron el ortofoto mosaico y los modelos digitales de superficie y terreno. Los rasters fueron generados con una resolución de 10cm.

 

Segundo paso: Calcular altura de objetos

Para generar un raster con las alturas de los edificios usamos la Calculadora ráster—Ayuda | ArcGIS Desktop disponible en la extensión Spatial Analyst. 

Básicamente restamos el terreno del superficie para obtener las alturas:

 

 

Para corregir celdas aisladas donde el terreno tiene un valor mayor a la superficie podemos aplicar la siguiente formula:

Con((Superficie - Terreno) < 0, 0, Superficie - Terreno)

 

Tercer paso: Asignar alturas a los edificios

Al contar con los polígonos de los edificios, podemos usar la herramienta Estadísticas zonales como tabla—Ayuda | ArcGIS Desktop la cual determina estadísticas de valores en un raster para zonas. En este caso nuestras zonas son los edificios y las estadísticas son valores mínimo, máximo y promedio (mean). Hacemos este proceso tanto para el terreno como para el superficie. Luego se puede usar la herramienta Campo de unión—Caja de herramientas Administración de datos | ArcGIS Desktop para unir los datos de las tablas de estadísticas a los edificios. Se puede usar el Calcular campo—Caja de herramientas Administración de datos | ArcGIS Desktop para calcular la altura al restar el promedio de terreno del promedio de superficie.

 

Se puede usar la herramienta De entidad a 3D por atributo—Ayuda | ArcGIS Desktop para convertir los edificios a entidades 3D. Esta herramienta está disponible con la extensión 3D Analyst.

 

 

Cuarto paso: Análisis de radiación solar

En la caja de herramientas Spatial Analyst se encuentre la herramienta Radiación solar de áreas—Ayuda | ArcGIS Desktop que nos permite generar un raster donde cada celda represente la radiación solar potencial. En este caso usamos el modelo digital de superficie (DSM) como insumo, puesto que queremos incluir el efecto de sombras de edificios altos sobre su alrededor. Además, queremos hacer el cálculo para todo un año para modelar la posición del sol durante el año y obviamente incluyendo su posición a diferentes horas del día. Abajo se muestre como se puede configurar este análisis: 

 

 

En este caso como resultado opcional también se va generar un raster con la cantidad de horas sol al año.

Horas sol por año. La zona en el círculo muestra el efecto de la sombra de los edificios altos sobre zonas cercanas. Mientras más rojo, más horas sol, las zonas azules representan zonas que reciben menos sol al año.

 

Raster de radiación solar expresado en vatios hora por metro cuadrado (WH/m²)

 

Usando la herramienta Estadísticas zonales como tabla—Ayuda | ArcGIS Desktop podemos determinar la radiación solar potencial por cada techo y usando  la herramienta Campo de unión—Caja de herramientas Administración de datos | ArcGIS Desktop podemos obtener estos datos en nuestra capa de edificios. 

 

 

Quinto paso, seleccionar techos y generar paneles solares

Para este ejercicio, el objetivo era escoger los techos que tienen mayor potencial de generar energía solar (radiación solar mayor a 2.000.000 WH/m² al año) y que tienen una superficie mínima de 50 m² . 

 

Con el fin de generar los paneles solares en los techos, se tuvo que definir una metodología para crear los rectángulos (2 x 1m) que representen los paneles. Para iniciar la generación se buscó para cada polígono (techo) el lado más largo. Este se obtiene después de aplicar una generalización para eliminar vertices redundantes y evaluar la distancia entre cada pareja de vertices. Una vez determinado el lado más largo, aplicando un espacio entre los paneles se generaron files de paneles y para cada panel se evalúe si queda completamente en el techo o no.

 

El resultado se escribe a un featureclass nuevo. 

 

Script de Python

El script de Python que usamos para crear los paneles solares en los techos es: 

#-------------------------------------------------------------------------------
# Name:        crear_paneles_solares_3d.py
# Purpose:
#
# Author:      xbakker
#
# Created:     20/05/2017
#-------------------------------------------------------------------------------

def main():
    import arcpy
    import os
    arcpy.env.overwriteOutput = True

    fc = r'C:\Drone2Map\Chapinero\gdb\Chapinero.gdb\Construcciones3D'
    fc_out = r'C:\Drone2Map\Chapinero\gdb\Chapinero.gdb\tmp12'
    fc_out_pol = r'C:\Drone2Map\Chapinero\gdb\Chapinero.gdb\Panels_v08' # Panels_v01
    fld_rad = 'RADmean'
    min_area = 50
    rad_min_value = 2000000

    # information to transfer
    fld_dtm = 'DSMmean'
    fld_dsm = 'DTMmean'
    fld_height = 'Height'
    fld_height1 = 'HeightPlus1'

    # panel settings
    pan_width = 2.0
    pan_height = 1.0
    pan_space_w = 0.2
    pan_space_h = 0.5

    fld_area = arcpy.Describe(fc).AreaFieldName
    where = '{0} >= {1} AND {2} >= {3}'.format(fld_rad, rad_min_value, fld_area, min_area)
    flds = ('OID@', 'SHAPE@', fld_dtm, fld_dsm, fld_height, fld_rad)
    dct_roofs = {r[0]: [r[1], r[2], r[3], r[4], r[5]] for r in arcpy.da.SearchCursor(fc, flds, where)}

    print len(dct_roofs.keys())

    feats = []
    feats_pol = []
    for oid, lst in dct_roofs.items():
        polygon = lst[0]
        dtm = lst[1]
        dsm = lst[2]
        height = lst[3]
        height1 = height + 1
        radiacion = lst[4]
        atts = [dtm, dsm, height, height1, radiacion]
        side, max_length = GetLargestLine(polygon, pan_space_w)
        print oid, max_length, side
        feats.append(side)
        angle = 90
        new_line = GetOffsetLine(side, angle, pan_space_h)
        inside = True

        if new_line.within(polygon, "BOUNDARY"):
            angle = 90
        else:
            angle = -90

        i = 0
        pols_panels = []
        while inside:
            start_dist = i * (pan_space_h + pan_height)  + pan_space_h
            end_dist = start_dist + pan_height
            new_line = GetOffsetLine(side, angle, start_dist)
            new_line2 = GetOffsetLine(side, angle, end_dist)
            panels = CreatePolygonsFromLines(new_line, new_line2, pan_width, pan_space_w, polygon)
            #polygon_out = CreatePolygonFromLines(new_line, new_line2)
            # if polygon_out.overlaps(polygon) == False:
            if len(panels) == 0:
                inside = False
            else:
                # feats_pol.append(polygon_out)
                pols_panels.extend(panels)
                feats.extend([new_line, new_line2])
            i += 1
        feats_pol.append([pols_panels, atts])

    # arcpy.CopyFeatures_management(feats, fc_out)

    # create output
    sr =  arcpy.Describe(fc).spatialReference
    ws, fc_name = os.path.split(fc_out_pol)
    arcpy.CreateFeatureclass_management(ws, fc_name, "POLYGON", None, None, None, sr)
    arcpy.AddField_management(fc_out_pol, fld_dtm, "DOUBLE")
    arcpy.AddField_management(fc_out_pol, fld_dsm, "DOUBLE")
    arcpy.AddField_management(fc_out_pol, fld_height, "DOUBLE")
    arcpy.AddField_management(fc_out_pol, fld_height1, "DOUBLE")
    arcpy.AddField_management(fc_out_pol, fld_rad, "DOUBLE")

    cnt = 0
    flds = ('SHAPE@', fld_dtm, fld_dsm, fld_height, fld_height1, fld_rad)
    with arcpy.da.InsertCursor(fc_out_pol, flds) as curs:
        for data in feats_pol:
            panels = data[0]
            atts = data[1]
            for polygon in panels:
                cnt += 1
                if cnt % 100 == 0:
                    print "Processing:", cnt
                row = (polygon, atts[0], atts[1], atts[2], atts[3], atts[4], )
                curs.insertRow(row)

    # arcpy.CopyFeatures_management(feats_pol, fc_out_pol)


def CreatePolygonFromLines(line1, line2):
    sr = line1.spatialReference
    pnt1, pnt2, = GetStartEndPoints(line1)
    pnt4, pnt3, = GetStartEndPoints(line2)
    return arcpy.Polygon(arcpy.Array([pnt1, pnt2, pnt3, pnt4, pnt1]), sr)


def CreatePolygonsFromLines(line1, line2, size, space, polygon):
    sr = line1.spatialReference
    length = line1.length
    fit = True
    i = 1
    feats = []
    while fit:
        end_dist = i*size + (i-1)*space
        if end_dist <= length:
            start_dist = end_dist - size
            seg1 = line1.segmentAlongLine(start_dist, end_dist, False)
            seg2 = line2.segmentAlongLine(start_dist, end_dist, False)
            pnt1, pnt2, = GetStartEndPoints(seg1)
            pnt4, pnt3, = GetStartEndPoints(seg2)
            panel = arcpy.Polygon(arcpy.Array([pnt1, pnt2, pnt3, pnt4, pnt1]), sr)
            if panel.within(polygon):
                feats.append(panel)
            i += 1
        else:
            fit = False
    return feats


def GetLargestLine(polygon, offset):
    max_length = 0
    prev_pnt = None
    side = None
    sr = polygon.spatialReference
    polygon2 = polygon.generalize(0.1)
    for part in polygon2:
        for pnt in part:
            # print pnt
            if prev_pnt is None:
                pass
            else:
                if pnt is None:
                    pass
                else:
                    length = GetDist(pnt, prev_pnt)
                    if length > max_length:
                        max_length = length
                        side = arcpy.Polyline(arcpy.Array([prev_pnt, pnt]), sr)
            prev_pnt = pnt

    length = side.length
    side = side.segmentAlongLine(offset, length-offset, False)
    return side, length


def GetDist(pnt1, pnt2):
    import math
    return math.hypot(pnt2.X - pnt1.X, pnt2.Y - pnt1.Y)


def GetAngle(pntg1, pntg2):
    '''determine angle of line based on start and end points geometries'''
    return pntg1.angleAndDistanceTo(pntg2, method='PLANAR')[0]

def GetOffsetLine(line, angle_offset, distance):
    sr = line.spatialReference
    pntg1, pntg2, = GetStartEndPointsGemetries(line)
    angle = GetAngle(pntg1, pntg2)
    angle_new = angle + angle_offset
    pntgla = pntg1.pointFromAngleAndDistance(angle_new, distance, 'PLANAR')
    pntglb = pntg2.pointFromAngleAndDistance(angle_new, distance, 'PLANAR')
    return arcpy.Polyline(arcpy.Array([pntgla.firstPoint, pntglb.firstPoint]), sr)

def GetStartEndPointsGemetries(line):
    sr = line.spatialReference
    pntg1 = arcpy.PointGeometry(line.firstPoint, sr, True, False)
    pntg2 = arcpy.PointGeometry(line.lastPoint, sr, True, False)
    return pntg1, pntg2

def GetStartEndPoints(line):
    return line.firstPoint, line.lastPoint

if __name__ == '__main__':
    main()

 

Carencias del script de Python

Los resultados del script no son perfectos y falta aplicar unos ajustes al script, para optimizar el resultado. Ver el ejemplo abajo:

 

En este caso se determine el lado más largo y se generen files de 7 paneles solares. Sin embargo, hacia el lado izquierdo queda espacio donde se pueden colocar otros paneles, pero este parte con el proceso actual queda sin paneles.

 

 

Crear labels 3D con Arcade

En ArcGIS Pro se pueden generar labels de tipo BillBoard que mantienen una buena legibilidad durante una animación.  El contenido del label se puede formatear usando un script de Arcade (ArcGIS Arcade | ArcGIS for Developers ): 

 

 

Script:

 

var kwh = Round($feature.Shape_Area / 4, 0) * 200;
var carbon = Round(kwh / 14280 * 6056, 0);
var trees = Round(kwh / 14280 * 157, 0);
var carkms = Round(kwh / 14280 * 23365, 0);

var label = 'Energía: ' + kwh + ' kWh\n' + 'Huella Carbon: ' + carbon + ' kg\n' + 'Arboles: ' + trees + '\nRadio carro: ' + carkms + ' km';
return label

 

Resultado:

 

y al rotar la vista:

 

Animación de los resultados:

 

 

Ejemplo de visor web 3D para el público

Esri desarrollo una aplicación web que se implementó en Dubai para compartir con el público que estimativo de rendimiento se puede obtener al implementar unos paneles solares cubriendo todo el techo. Para ver la aplicación se puede entrar a:

 

1 person found this helpful

Attachments

    Outcomes