Homec4science
Diffusion solar_potential (master) Browse

solarPOT/Panel_fitting/panel_fitting.py
cd293a1ab61cmaster

Last Change
Raw File

panel_fitting.py
View Options

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon
import os
import sys
import time
################# USER INPUTS ################
CORNER_FILE = sys.argv[1] # NOTE: file must be obtained from ArcGIS as described in notes from 17.10.2018
ROOFS_FILE = sys.argv[2] # NOTE: must be an exported shapefile containing at least DF_UID and FLAECHE
TARGET_DIR = sys.argv[3]
# panel dimensions
PANEL_WIDTH = float(sys.argv[4])
PANEL_HEIGHT = float(sys.argv[5])
################ CONSTANTS ###################
# special rules for flat roofs
FLAT_ROOF_THRESH = 0
FLAT_PANEL_TILT = 30
# define database columns to save with roof-fitted panels
DB_COLS = ['geometry', 'DF_UID', 'panel_tilt', 'origin_x', 'origin_y', 'cell_width', 'cell_height']
LOG_COUNT = 50 # at which iteration frequency should the time be displayed
try:
os.mkdir(TARGET_DIR)
except:
'Directory already exists - no new directory created (results will be overwritten)'
#============================================#
#============================================#
########## FUNCTION DEFINITIONS ##############
def rotate(angle,x,y):
# rotate a point (or a pandas dataframe) by a given angle
x_new = np.cos( np.deg2rad(angle) ) * x - np.sin( np.deg2rad(angle) ) * y
y_new = np.cos( np.deg2rad(angle) ) * y + np.sin( np.deg2rad(angle) ) * x
return x_new, y_new
########## Panel fitting algorithm ###########
def generate_panel_grid(rooftop):
# function that generates a grid of rectangles over a given rooftop. The grids are always started from the bottom-left corner of the (south-facing) rooftop
# "rooftop" must have the following attributes:
# - n_rows, n_cols: number of rows and columns of grid to be created, respectively
# - cell_width, cell_height: PROJECTED width and height of cells (original height multiplied by panel_tilt)
# - origin_x_rot, origin_y_rot: bottom left corner of the roof if rotated to face south
# - NEIGING: tilt of rooftop (different rule of placing panels for flat and for tilted roofs)
# - AUSRICHTUNG: orientation of rooftop (for de-rotation of origin)
# get input parameters
split_x = int(rooftop.n_cols) # number of splits in x-direction
split_y = int(rooftop.n_rows) # number of splits in y-direction
# define x and y coordinates of grid corners
y_coords = [ float(rooftop.origin_y_rot) + i * float(rooftop.cell_height) for i in range( split_y + 1 )]
x_coords = [ float(rooftop.origin_x_rot) + i * float(rooftop.cell_width) for i in range( split_x + 1 )]
# if roof is flat, skip every other row:
if rooftop.NEIGUNG.values[0] == 0:
range_y = range(0, split_y, 2)
else:
range_y = range(split_y)
# loop through all polygons and append coordinates of each corner (generate "fishnet")
polygons = []
for bottom_idx in range_y:
top_idx = bottom_idx + 1
for left_idx in range(split_x):
right_idx = left_idx + 1
poly_df = pd.DataFrame([[x_coords[ left_idx ], y_coords[ top_idx ]],
[x_coords[ right_idx ], y_coords[ top_idx ]],
[x_coords[ right_idx ], y_coords[ bottom_idx ]],
[x_coords[ left_idx ], y_coords[ bottom_idx ]]], columns = ['x','y'] )
rot = rotate(- int(rooftop.AUSRICHTUNG), poly_df.x, poly_df.y)
rotated_poly = [(rot[0][0],rot[1][0]),(rot[0][1],rot[1][1]),(rot[0][2],rot[1][2]),(rot[0][3], rot[1][3])]
polygons.append(Polygon(rotated_poly))
grid = gpd.GeoDataFrame({ 'geometry': polygons })
grid.crs = rooftop.crs
return grid
##############################################
# ========================================== #
##### Project corners and get grid origin ####
timer = time.time()
# load data and set flat roofs as south-facing (assume that on all flat roofs solar panels are oriented south)
roof_corners = pd.read_csv(CORNER_FILE)
roof_corners.loc[roof_corners.NEIGUNG == 0,'AUSRICHTUNG'] = 0
# rotate all corners, so that surface are facing south
roof_corners['x_rot'], roof_corners['y_rot'] = rotate(roof_corners.AUSRICHTUNG, roof_corners.POINT_X, roof_corners.POINT_Y)
# group rooftop data per rooftop and compute minimum and maximum values
roof_grouped = roof_corners.groupby('DF_UID')
roof_min = roof_grouped.min()
roof_max = roof_grouped.max()
# get essential information for each roof
rooftops = roof_min.loc[:,['AUSRICHTUNG','NEIGUNG']]
# === ADD ADDITIONAL INFORMATION
# add panel tilt, as this is different for flat roofs
rooftops['panel_tilt'] = rooftops.NEIGUNG
rooftops.loc[rooftops.NEIGUNG <= FLAT_ROOF_THRESH, 'panel_tilt'] = FLAT_PANEL_TILT
# rotate bottom left back (same angle magnitude, but opposite direction)
# origin_x/y is the origin coordinate, axis_x/y is the coordinate for y_axis_coord in create fishnet (only for ArcGIS algorithm!)
rooftops['origin_x'], rooftops['origin_y'] = rotate( -rooftops.AUSRICHTUNG, roof_min.x_rot, roof_min.y_rot )
# rooftops['axis_x' ], rooftops['axis_y' ] = rotate( -rooftops.AUSRICHTUNG, roof_min.x_rot, roof_max.y_rot )
# calculate cell width & cell height
rooftops['cell_width' ] = PANEL_WIDTH
rooftops['cell_height'] = PANEL_HEIGHT * np.cos(np.deg2rad( rooftops.panel_tilt ))
# calculate number of rows and number of columns
roof_extent = roof_max - roof_min
rooftops['n_cols'] = (np.ceil( roof_extent.x_rot / rooftops.cell_width )).astype(int)
rooftops['n_rows'] = (np.ceil( roof_extent.y_rot / rooftops.cell_height)).astype(int)
# add rotated origin (for grid creation, i.e. bottom-left corner of south-facing roof)
rooftops['origin_x_rot'] = roof_min.x_rot
rooftops['origin_y_rot'] = roof_min.y_rot
# add number of corners as indication of roof shape complexity
rooftops['n_corners'] = roof_grouped.count().ORIG_FID - 1 # first point of shape always counted twice!! (as lines touch)
output_file = os.path.join(TARGET_DIR, 'panel_control.csv')
rooftops.to_csv(output_file)
print('Saved panel control file in %s' %output_file)
print('Time for computing rooftop control file: %fs\n' %(time.time()-timer))
# ========================================== #
##### Create panel grids on rooftops #########
# read data and merge with rooftop control information (to find rooftops for which match can be performed)
match_roofs = gpd.read_file(ROOFS_FILE)
match_roofs = match_roofs.loc[:, ['geometry', 'DF_UID','FLAECHE']].merge(rooftops, how = 'inner', on = 'DF_UID')
timer = time.time()
# loop through all rooftops, generate grid and find panels inside the roof area
panels = []
for idx, DF_UID in enumerate(match_roofs.DF_UID):
if idx % LOG_COUNT == 0:
print('Iteration %d out of %d - Loop duration: %f seconds' %(idx, len(match_roofs), time.time()-timer))
# get rooftop surface to match panels on
ROOF = match_roofs[ match_roofs.DF_UID == DF_UID]
# get panel polygons
grid = generate_panel_grid(ROOF)
# find panels inside rooftop and add to all_panels
inside_roofs = gpd.tools.sjoin(grid, ROOF, op = 'within')
panels.append(inside_roofs)
print('Total computational time for computing panels: %fs\n' %(time.time()-timer))
# merge all panels to one geodataframe
panels = pd.concat(panels).reset_index().drop('index', axis = 1)
print('Created GeoDataFrame of panels')
output_file = os.path.join(TARGET_DIR, 'fitted_panels.shp')
panels.loc[:,DB_COLS].to_file(output_file)
print('Saved shapefile with fitted panels in %s\n' %output_file)
# ========================================== #
############ Evaluate results ################
timer = time.time()
# evaluate the panelled area statistics per rooftop (group by rooftops)
panels['panel_area'] = PANEL_WIDTH * PANEL_HEIGHT
panel_groups = panels.groupby('DF_UID')
# get information on count and panelled area
panel_count = panel_groups.count()['geometry']
panel_count.name = 'panel_count'
panelled_area = panel_groups.sum()['panel_area']
panelled_area.name = 'panelled_area'
# get relevant information both from panels and from rooftops
panel_info = pd.DataFrame( [panel_count, panelled_area]).transpose()
roof_info = match_roofs.loc[:,['DF_UID', 'FLAECHE', 'panel_tilt', 'NEIGUNG','n_corners']].set_index('DF_UID')
# merge information and fill missing information with zeors (no panels detected)
panel_eval = panel_info.merge(roof_info, how = 'outer', left_index = True, right_on = 'DF_UID').sort_index()
panel_eval.fillna(0, inplace = True)
# compute panelled area ratio
panel_eval['panelled_area_ratio'] = panel_eval.panelled_area / panel_eval.FLAECHE
output_file = os.path.join(TARGET_DIR, 'panelled_area_stats.csv')
panel_eval.to_csv(output_file)
print('Saved panelled area statistics file in %s' %output_file)
print('Time for computing panelled area statistics: %fs\n' %(time.time()-timer))
c4science · Help