Homec4science
Diffusion solarPV (aw) Browse

solarPV_identification/Available_area/panel_fitting_batch_new.py
20a1108820e7aw

Last Change
Raw File

panel_fitting_batch_new.py
View Options

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon
import fiona
import geo_ops
import os
import glob
import sys
import time
################# USER INPUTS ################
# file information
ROOFS_FILE = sys.argv[1] # NOTE: must be a shapefile
TARGET_DIR = sys.argv[2]
# Current batch information
BATCH_ID = int(sys.argv[3])
BATCH_SIZE = int(sys.argv[4])
# panel dimensions
PANEL_WIDTH = float(sys.argv[5])
PANEL_HEIGHT = float(sys.argv[6])
################ CONSTANTS ###################
# special rules for flat roofs
FLAT_ROOF_THRESH = 5 # check which roofs count as "flat" --> OUTCOME OF VALIDATION: USE 5°
FLAT_PANEL_TILT = 15 # panel tilt of flat roofs --> OUTCOME OF VALIDATION: GOOD FIT
# field names in shapefile
AREA_FIELD = 'FLAECHE'
TILT_FIELD = 'NEIGUNG'
ORIENTATION_FIELD = 'aspect'
BUFF_DIST = 0.0
TOLERANCE = 0.1
LOG_COUNT = 500 # at which iteration frequency should the time be displayed in loops
SAVE_PANELS = True
try:
os.mkdir(TARGET_DIR)
except:
'Directory already exists - no new directory created (results will be overwritten)'
#============================================#
#============================================#
############### Main function ###############
def main():
print_inputs()
timer = time.time()
# load data and perform pre-processing steps (clean & rotate)
roofs, rotated_roofs = load_from_file(ROOFS_FILE, batch_id = BATCH_ID, batch_size = BATCH_SIZE, tolerance = TOLERANCE)
print("Loaded data in %.2f seconds\n" %(time.time() - timer))
timer = time.time()
# create a buffer to obtain the shapes that the panels will be fitted on
buffered_roofs = make_buffer(rotated_roofs, BUFF_DIST)
buffered_roofs.to_file( os.path.join(TARGET_DIR, "TMP_%d_roofs_buf.shp" %BATCH_ID) )
print("Created buffer in %.2f seconds\n" %(time.time() - timer))
timer = time.time()
# obtain information necessary for fitting panels (panel dimensions etc.)
panelling_info = get_panel_fitting_info(buffered_roofs)
panelling_info.to_csv( os.path.join(TARGET_DIR, "TMP_%d_panel_info.csv" %BATCH_ID), index = False)
print("Computed info for panel fitting in %.2f seconds\n" %(time.time() - timer))
timer = time.time()
# fit panels on rooftops and evaluate counts
panel_count, panels, matched_roofs = fit_panels(rotated_roofs, panelling_info, buff_dist = BUFF_DIST, tolerance = TOLERANCE)
panels.crs = roofs.crs
if SAVE_PANELS: panels.to_file( os.path.join(TARGET_DIR, "PANELS_%d.shp" %BATCH_ID) )
print("Computed panel counts for rooftops and saved panels in %.2f seconds\n" %(time.time() - timer))
timer = time.time()
# compute statistics for panel counts (panelled_area_ratio)
panel_stats = compute_panel_stats(panel_count, roofs, matched_roofs)
panel_stats.to_csv( os.path.join(TARGET_DIR, 'panelled_area_stats_%d.csv' %BATCH_ID), index = False )
print("Saved %s in %.2f seconds\n" %(os.path.join(TARGET_DIR, 'panelled_area_stats_%d.csv' %BATCH_ID), time.time() - timer))
print(panel_stats.head())
# remove temporary files
[os.remove(x) for x in glob.glob( os.path.join(TARGET_DIR, "TMP_%d_*" %BATCH_ID) )]
print('Returning from main')
################ SUBROUTINES ################
def load_from_file(filepath, batch_id, batch_size, tolerance):
# load a given batch of lines from file and pre-process the data
# returns geopandas dataframe with polygons and database information
# open file with fiona
# file = fiona.open(filepath)
# load selected rows from file
start = batch_id * batch_size
end = (batch_id + 1) * batch_size
# roofs = gpd.GeoDataFrame.from_features(file[start : end]) --> DEPRECATED
# roofs.crs = file.crs
roofs = gpd.read_file(filepath, rows = slice(start, end))
if len(roofs) == 0:
raise ValueError('Rooftop dataframe is empty - no data loaded from file.')
# PRE-PROCESSING:
# part 1: check that roof orientation is correct (flat roofs should be south-facing (0))
# roofs.loc[roofs[TILT_FIELD] <= FLAT_ROOF_THRESH, ORIENTATION_FIELD] = 0 # --> BETTER FIT IF PLACED ALONG ASPECT
# part 2: simplify roof shapes
roofs['geometry'] = roofs.simplify(tolerance = tolerance)
# part 3: add panel tilt information to roofs (different from tilt for flat roofs)
roofs['panel_tilt'] = roofs[ TILT_FIELD ]
roofs.loc[roofs[ TILT_FIELD ] <= FLAT_ROOF_THRESH, 'panel_tilt'] = FLAT_PANEL_TILT
# part 4: rotate roofs to face south
rotated_roofs = geo_ops.rotate_by_field(roofs, ORIENTATION_FIELD, origin = (0,0))
return roofs, rotated_roofs
def make_buffer(roofs, buff_dist, cap_style = 2, join_style = 2):
# create a buffer of width buff_dist INSIDE the polygon
return geo_ops.create_buffer( roofs,
buff_dist = - buff_dist,
cap_style = cap_style,
join_style = join_style )
def get_panel_fitting_info(roofs):
# create a pandas dataframe (without geometries) that contains all data necessary to fit panels
# required fields: tilt, orientation, n_rows, n_cols, cell_width, cell_height, origin_x_rot, origin_y_rot
# Step 1: rotate shapes around the origin and get the bounding box of each shape (facing south)
roof_bounds = roofs.bounds
# Step 2: create dataframe with rooftop data and get additional information
panel_info_df = roofs.loc[:, ['DF_UID', ORIENTATION_FIELD, TILT_FIELD, 'panel_tilt']]
# calculate cell width & cell height
panel_info_df['cell_width' ] = PANEL_WIDTH
panel_info_df['cell_height'] = PANEL_HEIGHT * np.cos(np.deg2rad( panel_info_df.panel_tilt ))
# calculate number of rows and number of columns
x_extent = roof_bounds.maxx - roof_bounds.minx
y_extent = roof_bounds.maxy - roof_bounds.miny
panel_info_df['n_cols'] = (np.ceil( x_extent / panel_info_df.cell_width )).astype(int)
panel_info_df['n_rows'] = (np.ceil( y_extent / panel_info_df.cell_height)).astype(int)
# add rotated origin (for grid creation, i.e. bottom-left corner of south-facing roof)
panel_info_df['origin_x'] = roof_bounds.minx
panel_info_df['origin_y'] = roof_bounds.miny
return panel_info_df
def fit_panels(roofs, panel_info, buff_dist, tolerance):
# fit panels on rooftops based on information in panel_info
# uses a tolerance around the buff_dist to find the spatial joins
# returns:
# - panel count: Pandas dataframe with columns DF_UID, panel_count
# - panels: GeoDataFrame of panel polygons (in original rotation) with fields geometry, DF_UID
# - matched_roofs: rooftop shapes based on which the panel fitting was performed
# Step 1: create roofs to match on (buffered with buff_dist - tolerance)
shapes_to_match = geo_ops.create_buffer( roofs,
buff_dist = - (buff_dist - tolerance),
columns = ['geometry', 'DF_UID'],
cap_style = 2,
join_style = 2 )
match_roofs = shapes_to_match.merge(panel_info, how = 'inner', on = 'DF_UID')
timer = time.time()
loop_timer = time.time()
# Step 2: loop through all rooftops, generate grid and find panels inside the roof area
panel_count = []
panels = []
for idx, DF_UID in enumerate(match_roofs.DF_UID):
# 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')
panel_count.append([DF_UID, len(inside_roofs)])
panels.append(inside_roofs)
if (idx+1) % LOG_COUNT == 0:
print('Computing panel counts: Iteration %d out of %d - Loop duration: %.2f seconds' %(idx+1, len(match_roofs), time.time()-loop_timer))
curr_data = pd.DataFrame(np.asarray(panel_count), columns = ['DF_UID', 'panel_count'])
curr_data.to_csv(os.path.join(TARGET_DIR, "TMP_%d_panel_count.csv" %BATCH_ID))
# convert panel count to dataframe
panel_count = pd.DataFrame(np.asarray(panel_count), columns = ['DF_UID', 'panel_count'])
print('Computed all panel counts - time since loop start: %.2fs' %(time.time() - timer))
if SAVE_PANELS:
# merge all panels to one geodataframe, rotate them back and remove unncessary data
panels = pd.concat(panels).reset_index().drop('index', axis = 1)
panels_rotated = geo_ops.rotate_by_field(panels, ORIENTATION_FIELD, origin = (0,0), backwards = True)
panels_out = panels_rotated.loc[:, ['geometry', 'DF_UID']]
panels_out.crs = match_roofs.crs
print('Created panel shapes - time since loop start: %.2fs' %(time.time() - timer))
else:
panels_out = []
return panel_count, panels_out, match_roofs
def compute_panel_stats(panel_count, roofs, matched_roofs, add_ncorners = True):
# compute panelled area and panelled area ratio
# returns dataframe with the following fields:
# - DF_UID, TILT_FIELD, AREA_FIELD, panel_tilt : directly from roofs geodataframe
# - panel_count : from fitting panels on rooftops
# - panelled_area, panelled_area_ratio
# - n_corners : number of feature vertices, can be used as additional ML feature
# compute the panelled area for each rooftop
panel_area = matched_roofs.cell_width * matched_roofs.cell_height / np.cos(np.deg2rad(matched_roofs[TILT_FIELD]))
panel_count['panelled_area'] = panel_count.panel_count * panel_area.values
# add panel counts to basic rooftop information and set all NaNs to 0
# (counts for roofs that were not in match_roofs)
all_roofs = roofs.loc[:, ['DF_UID', AREA_FIELD, TILT_FIELD, 'panel_tilt']]
panel_eval = all_roofs.merge( panel_count, on = 'DF_UID', how = 'left')
panel_eval.fillna(0, inplace = True)
# compute panelled area ratio
panel_eval['panelled_area_ratio'] = panel_eval.panelled_area / panel_eval[ AREA_FIELD ]
if add_ncorners:
try:
panel_eval['n_corners'] = roofs.geometry.exterior.apply(lambda x: len(x.coords)-1)
except:
print('Adding n_corners failed: Likely due to failure because of multipart polygons')
return panel_eval
############ AUXILARY FUNCTIONS ##############
def print_inputs():
print('\nINPUT DATA:')
print("Input file: %s" %ROOFS_FILE)
print("Output directory: %s" %TARGET_DIR)
print("Batch ID: %d" %BATCH_ID)
print("Batch size: %d" %BATCH_SIZE)
print("Panel width: %f" %PANEL_WIDTH)
print("Panel height: %f" %PANEL_HEIGHT)
print("")
def rotate(angle,x,y):
# rotate a point (or pandas Series) by.a given angle around the origin
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
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
# - TILT_FIELD: tilt of rooftop (different rule of placing panels for flat and for tilted roofs)
# 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) + i * float(rooftop.cell_height) for i in range( split_y + 1 )]
x_coords = [ float(rooftop.origin_x) + i * float(rooftop.cell_width) for i in range( split_x + 1 )]
# if roof is flat, skip every other row: --> DONT DO THIS ANY MORE (EW SCENARIO)
#if rooftop[TILT_FIELD].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
Xleft = x_coords[ left_idx ]
Xright = x_coords[ right_idx ]
Ytop = y_coords[ top_idx ]
Ybottom = y_coords[ bottom_idx]
polygons.append(Polygon([(Xleft, Ytop), (Xright, Ytop), (Xright, Ybottom), (Xleft, Ybottom)]))
grid = gpd.GeoDataFrame({ 'geometry': polygons })
grid.crs = rooftop.crs
return grid
##############################################
if __name__ == '__main__':
main()
c4science · Help