panel_fitting.py
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, May 14, 03:03

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))

panel_fitting.pyNo OneTemporaryActions

File Metadata

panel_fitting.pyView Options

Event Timeline

panel_fitting.py
No OneTemporary
Actions

panel_fitting.py
View Options