Homec4science
Diffusion solarPV (master) Browse

solarPV_identification/PV_postprocessing/spatial_join_pv.py
321a85653075master

Last Change
Raw File

spatial_join_pv.py
View Options

import numpy as np
import pandas as pd
import geopandas as gpd
import os
import sys
import time
import shutil
#### Function definitions
def coords_to_point( df, coord_names = ['x', 'y'] ):
xy = df[ coord_names ].values
return gpd.GeoDataFrame( df, geometry = gpd.points_from_xy(xy[:,0], xy[:,1]))
def load_file(fp, columns = None, rows = None, crs = None, bbox = None, **kwargs):
if columns is not None:
all_cols = gpd.read_file( fp, rows = 0, **kwargs ).columns
ignore_cols = list( all_cols[ np.logical_not( all_cols.isin( columns ))])
else:
ignore_cols = None
return gpd.read_file( fp, ignore_fields = ignore_cols, rows = rows, bbox = bbox, **kwargs)
#### Variable definitions
## INPUTS
PV_FP = sys.argv[1] # '/work/hyenergy/output/CNN_pv_predictions/pv_predictions_geneva.csv'
MASK_FP = sys.argv[2] # geospatial limitation (geodataframe) if perimeter of analysis
OUT_FP = sys.argv[3] # '/scratch/walch/pv_linking/PV_detected.csv'
# Optional inputs
if len(sys.argv) > 4:
BLD_BUFFER = float(sys.argv[4]) # Distance by which to buffer BLD footprint shapes
else:
BLD_BUFFER = 3
if len(sys.argv) > 6:
DEFAULT_TILT = int(sys.argv[5])
DEFAULT_ASPECT = int(sys.argv[6])
else:
DEFAULT_TILT = 0 # Default tilt angle assumed on roofs without rooftop info
DEFAULT_ASPECT = -180 # Default roof aspect assumed on roofs without rooftop info
print('\nDefault tilt for roofs without rooftop info: %d' %DEFAULT_TILT)
print('Default aspect for roofs without rooftop info: %d' %DEFAULT_ASPECT)
print('BLD buffer width: %dm' %BLD_BUFFER)
print('Load PV pixels from %s' %PV_FP)
print('Load mask from %s' %MASK_FP)
print('Save output to %s\n' %OUT_FP)
## Constants
BATCH_SIZE = 100000 # size of batches to process spatial join (executed in for-loop)
PIXEL_RES = 0.25
# CONSTANT FILES (national files)
# TLM_PATH = '/work/hyenergy/raw/Building_data/swissTLM3D/Update swissTLM3D 2017/LV95/swissTLM3D_LV95_LN02.gdb'
# TLM_LAYR = 'TLM_GEBAEUDE_FOOTPRINT'
# TLM_COLS = ['UUID', 'geometry']
BLD_PATH = '/work/hyenergy/output/buildings/buildings_CH/BLD_LV95_CH'
ROOF_PATH = '/work/hyenergy/raw/Sonnendach/Sonnendach_all_roofs/Sonnendach_CH.gdb'
ROOF_LAYR = 'SOKAT_CH_DACH'
# Roof information to inherit ('FLAECHE', 'NEIGUNG', 'AUSRICHTUNG' are compulsory)
ROOF_COLS = ['geometry', 'DF_UID', 'FLAECHE', 'NEIGUNG', 'AUSRICHTUNG']
roof_col_dict = {'FLAECHE' : 'ROOF_AREA', 'NEIGUNG' : 'ROOF_TILT', 'AUSRICHTUNG' : 'ROOF_ASPECT'}
# Files that link to BLD_ID
ROOF_LINK = '/work/hyenergy/output/buildings/buildings_CH/SOLKAT_DACH_CH_BLD_ID.csv'
# TLM_LINK = '/work/hyenergy/output/buildings/buildings_CH/TLM_CH_BLD_ID.csv'
TMP_DIR = os.path.join( os.path.split(OUT_FP)[0], 'tmp' )
try:
os.mkdir( TMP_DIR )
except:
print('%s exists' %TMP_DIR)
#### Load data
## PV pixels
COORDS = ['x', 'y']
PV_data = pd.read_csv( PV_FP, header = None, names = COORDS )
PV_data['pixel_ID'] = list(range(len(PV_data)))
# assign the pixel area
PV_data['pixel_area'] = PIXEL_RES**2
print(PV_data.head())
## Load mask
MASK = gpd.read_file( MASK_FP )
## Load roofs
tt = time.time()
# load roof polygons & link information to BLD
ROOFS = load_file( ROOF_PATH , columns = ROOF_COLS, layer = ROOF_LAYR, mask = MASK )
roof_link = pd.read_csv(ROOF_LINK, index_col = 0).drop(columns = ['batchID'])
print('\nLoaded roof data in %.2fs' %(time.time()-tt))
# Merge BLD_ID to roofs
ROOFS = ROOFS.merge(roof_link, on = 'DF_UID', how = 'left')
# rename columns
# NOTE: roof aspect defined as clockwise degrees from south (flat roofs are assigned a north aspect)
ROOFS = ROOFS.rename( roof_col_dict, axis = 1 )
print(ROOFS.head())
## Buffered BLD data
tt = time.time()
# Load BLD data
BLD = load_file( BLD_PATH, mask = MASK ).drop(columns = ['batchID'])
print('\nLoaded BLD data in %.2fs' %(time.time()-tt))
# Buffer by BLD_BUFFER meters
BLD_BUF = BLD.copy()
BLD_BUF['geometry'] = BLD.buffer( BLD_BUFFER )
BLD_BUF['ROOF_AREA'] = BLD.geometry.area # non-buffered area
BLD_BUF['ROOF_TILT'] = DEFAULT_TILT
BLD_BUF['ROOF_ASPECT'] = DEFAULT_ASPECT
print(BLD_BUF.head())
#### Spatial join pixels to roofs/buildings
N_BATCHES = int(np.ceil( len(PV_data)/ BATCH_SIZE ))
print('\nAssigning pixels to roofs/buildings in %d batches' %N_BATCHES)
all_pixel_roof = []
all_pixel_bld = []
for batch in range( N_BATCHES ):
t0 = time.time()
start = batch * BATCH_SIZE
end = (batch + 1) * BATCH_SIZE
# Convert pixel data to point data
tt = time.time()
PV_pixels = coords_to_point( PV_data[ start:end ].copy(), coord_names = COORDS )
print('Converted coords to point data in %.2fs' %(time.time()-tt))
# Add the coordinate system of PV_pixels
PV_pixels.crs = {'init': 'epsg:21781'}
PV_pixels.crs
if PV_pixels.crs != ROOFS.crs:
tt = time.time()
PV_pixels = PV_pixels.to_crs( ROOFS.crs )
print('Reprojected data in %.2fs' %(time.time()-tt))
# spatial join pixels and roofs (keep only intersecting pixels)
tt = time.time()
pxl_per_roof = gpd.sjoin( PV_pixels, ROOFS, op = 'within' ).drop( columns = ['index_right', 'geometry'] )
print('Spatial joined pixels with roofs in %.2fs' %(time.time()-tt))
if pxl_per_roof.duplicated( subset = COORDS ).any():
print('WARNING: SOME PIXELS ARE DUPLICATED')
all_pixel_roof.append( pxl_per_roof )
### Handle remaining roofs
remaining_pixels = PV_pixels[ np.logical_not(PV_pixels.pixel_ID.isin( pxl_per_roof.pixel_ID.values )) ]
if len(remaining_pixels) == 0:
print('No remaining pixels')
continue
if PV_pixels.crs != BLD_BUF.crs:
print('Reprojecting data')
PV_pixels = PV_pixels.to_crs( BLD_BUF.crs )
# spatial join pixels and buffered buildings (keep only intersecting pixels)
tt = time.time()
pxl_per_bld = gpd.sjoin( remaining_pixels, BLD_BUF, op = 'within' ).drop( columns = ['index_right', 'geometry'] )
print('Spatial joined pixels with buffered buildings in %.2fs' %(time.time()-tt))
# Handle duplicates
if pxl_per_bld.duplicated( subset = COORDS ).any():
# Extract the duplicates
duplicates = pxl_per_bld[pxl_per_bld.duplicated( subset = COORDS, keep = False )]
px_count = {}
# Get the number of pixels on each of the duplicate buildings
for bld in duplicates.BLD_ID.unique():
px_count[bld] = len( pxl_per_bld[pxl_per_bld.BLD_ID == bld] )
# For each duplicate, choose the one with most pixels
duplicates['px_count'] = duplicates.BLD_ID.map( px_count )
bld_w_highest_pxCount = ( duplicates.sort_values('px_count', ascending = False)
.drop_duplicates( subset = COORDS, keep = 'first' ).drop( columns = ['px_count'] ) )
# Add selected values to the pxl_per_bld
pxl_per_bld = pd.concat([ pxl_per_bld.drop_duplicates( subset=COORDS, keep=False ), bld_w_highest_pxCount ])
if pxl_per_bld.duplicated( subset = COORDS ).any():
print('WARNING: SOME PIXELS ARE DUPLICATED')
all_pixel_bld.append( pxl_per_bld )
# Save intermediate data
pxl_per_roof.to_csv( os.path.join( TMP_DIR, 'pxl_per_roof_%d.csv' %batch ), index = False )
pxl_per_bld .to_csv( os.path.join( TMP_DIR, 'pxl_per_bld_%d.csv' %batch ), index = False )
print('Finished iteration %d in %.2fs\n' %(batch, time.time()-t0) )
# Merge all batches
all_pixel_roof = pd.concat(all_pixel_roof)
print(all_pixel_roof.head())
all_pixel_bld = pd.concat(all_pixel_bld )
print(all_pixel_bld.head())
# Merge building and roof data
all_pixels = pd.concat([all_pixel_roof, all_pixel_bld], sort = False)
# Computed the tilted PV area
all_pixels['PV_area'] = all_pixels['pixel_area'] / np.cos( np.deg2rad( all_pixels['ROOF_TILT'] ) )
print(all_pixels.head())
print('\nFiltered pixels: %d of %d\n' %(len(PV_data), len(all_pixels)))
tt = time.time()
all_pixels.to_csv( OUT_FP, index = False )
print('Saved %s - DONE in %.2fs' %(OUT_FP, time.time()-tt))
## Try to remove tree; if failed show an error using try...except on screen
try:
shutil.rmtree(TMP_DIR)
except OSError as e:
print ("Error: %s - %s." % (e.filename, e.strerror))
c4science · Help