Generación de Paneles Solares en techos Bogotá
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 SPX 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.
- target="_self">Datos de entrada
- Primer paso: Procesar las imágenes con Drone2Map
- Segundo paso: Calcular altura de objetos
- Tercer paso: Asignar alturas a los edificios
- Cuarto paso: Análisis de radiación solar
- Quinto paso, seleccionar techos y generar paneles solares
- Script de Python
- Carencias del script de Python
- Crearlabels 3D con Arcade
- Ejemplo de visor web 3D para el público
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 labelResultado:
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:
