Homec4science
Diffusion solarPV (master) Browse

solarPV_identification/PV_postprocessing/Create_BLD/create_BLD_ID_batch_V2.py
cdf83e4b6b33master

Last Change
Raw File

create_BLD_ID_batch_V2.py
View Options

import numpy as np
import pandas as pd
import xarray as xr
import geopandas as gpd
import os
import time
import sys
from geo_util import *
## ==== USER INPUT ==== ##
OUT_FOLDER = sys.argv[1]
BATCH_PATH = sys.argv[2] # '/scratch/walch/pv_linking/batch25km_PLY'
BATCH_DELTA = int(sys.argv[3]) # 25000
BATCH_ID = int(sys.argv[4])
## ==== CONSTANTS ==== ##
ID_OFFSET = 1e7 # multiplied with BATCH_ID to obtain offset of BLD_ID to avoid overlapping IDs
BUFF_DIST = 1000 # buffer distance for overlapping geometries
TOLERANCE = 0.001 # For geometry simplification;
TOL_OVERLAY = 0.01 # to avoid overlaying issues
TLM_PATH = '/work/hyenergy/raw/Building_data/swissTLM3D/Update swissTLM3D 2017/LV95/swissTLM3D_LV95_LN02.gdb'
TLM_LAYR = 'TLM_GEBAEUDE_FOOTPRINT'
TLM_COLS = ['OBJECTID', 'UUID', 'geometry']
ROOF_PATH = '/work/hyenergy/raw/Sonnendach/Sonnendach_all_roofs/Sonnendach_CH.gdb'
ROOF_LAYR = 'SOKAT_CH_DACH'
ROOF_COLS = ['geometry', 'DF_UID', 'SB_UUID']
## ==== EXECUTION ==== ##
# # Load data
# Load and preprocess batches to obtain mask
batches = gpd.read_file(BATCH_PATH)
batch_core = batches[batches.batchID == BATCH_ID]
mask = batch_core.copy()
mask['geometry'] = batch_core.geometry.buffer(BUFF_DIST, join_style = 2)
print(mask)
# TLM
tt = time.time()
TLM = load_file( TLM_PATH, layer = TLM_LAYR, columns = TLM_COLS, bbox = mask )
print('Loaded TLM in %.2fs' %(time.time()-tt))
invalid_geometry = TLM[TLM.geometry.is_valid == False ]
TLM.loc[invalid_geometry.index, 'geometry'] = invalid_geometry.geometry.buffer(0)
print('Validated TLM geometries')
# Rooftops
tt = time.time()
ROOFS = load_file( ROOF_PATH, layer = ROOF_LAYR, columns = ROOF_COLS, bbox = mask )
print('Loaded roofs in %.2fs' %(time.time()-tt))
invalid_geometry = ROOFS[ROOFS.geometry.is_valid == False ]
ROOFS.loc[invalid_geometry.index, 'geometry'] = invalid_geometry.geometry.buffer(0)
print('Validated roof geometries')
if (len(TLM) == 0) & (len(ROOFS) == 0):
print('No buildings found')
exit()
# # Dissolve and split polygons
TLM_BLD = dissolve_overlapping_polygons(pd.concat([TLM[['geometry']], ROOFS[['geometry']]]))
TLM_BLD = TLM_BLD.drop(columns = 'index').reset_index().rename({'index' : 'BLD_ID'}, axis = 1)
TLM_BLD.crs = TLM.crs
tt = time.time()
TLM_BLD['geometry'] = TLM_BLD['geometry'].simplify(TOLERANCE)
print('Simplified geometries in %.2fs' %(time.time()-tt))
tt = time.time()
TLM_BLD.to_file( os.path.join(OUT_FOLDER, 'TLM_dissolved_batch%d' %BATCH_ID) )
print('Saved %s in %.2fs' %(os.path.join(OUT_FOLDER, 'TLM_dissolved_batch%d' %BATCH_ID), time.time()-tt))
# # Link SB_UUID to BLD_ID
tt = time.time()
TLM_w_ID = gpd.sjoin( TLM, TLM_BLD, op = 'intersects' ).drop(columns = ['index_right'])
print('Spatially joined TLM to BLD in %.2fs' %(time.time()-tt))
# Check for duplicates and remove by overlay
dupl = TLM_w_ID[TLM_w_ID.duplicated(subset = ['UUID'], keep = False)]
dupl_UUIDs = dupl.UUID.unique()
if len(dupl_UUIDs) > 0:
dupl_TLM_w_ID = overlay_and_clean(TLM[TLM.UUID.isin(dupl_UUIDs)], TLM_BLD, tolerance = TOL_OVERLAY)
TLM_w_ID = pd.concat([ TLM_w_ID.drop_duplicates( subset = ['UUID'], keep = False ), dupl_TLM_w_ID ])
else:
print('No duplicates found')
print('Lengths of TLM subsets: TLM_w_ID, TLM, TLM_w_ID.dropna()')
print(len(TLM_w_ID))
print(len(TLM))
print(len(TLM_w_ID.dropna()))
TLM_ID_info = pd.DataFrame( TLM_w_ID.drop(columns = ['geometry']) )
print(TLM_ID_info.head())
TLM_ID_info.to_csv( os.path.join(OUT_FOLDER, 'TLM_BLD_ID_batch%d.csv' %BATCH_ID) , index = False)
print('Saved %s' %os.path.join(OUT_FOLDER, 'TLM_BLD_ID_batch%d.csv' %BATCH_ID))
# # Intersect with roofs
# ### Merge UUIDs
tt = time.time()
ROOFS_w_ID = gpd.sjoin( ROOFS, TLM_BLD, op = 'intersects' ).drop(columns = ['index_right'])
print('Spatially joined ROOFS to BLD in %.2fs' %(time.time()-tt))
# Check for duplicates and remove by overlay
dupl = ROOFS_w_ID[ROOFS_w_ID.duplicated(subset = ['DF_UID'], keep = False)]
dupl_DF_UIDs = dupl.DF_UID.unique()
if len(dupl_DF_UIDs) > 0:
dupl_ROOF_w_ID = overlay_and_clean(ROOFS[ROOFS.DF_UID.isin(dupl_DF_UIDs)], TLM_BLD, tolerance = TOL_OVERLAY)
ROOFS_w_ID = pd.concat([ ROOFS_w_ID.drop_duplicates( subset = ['DF_UID'], keep = False ),
dupl_ROOF_w_ID ]).sort_values(by = 'DF_UID')
else:
print('No duplicate roofs')
print('Lengths of remaining ROOF subsets: ROOFS_w_ID, ROOFS_w_ID.dropna(), ROOFS, ROOFS_w_ID.SB_UUID.unique(), ROOFS.SB_UUID.unique()')
print(len(ROOFS_w_ID))
print(len(ROOFS_w_ID.dropna()))
print(len(ROOFS))
print(len(ROOFS_w_ID.SB_UUID.unique()))
print(len(ROOFS.SB_UUID.unique()))
# ## Rooftops
ROOFS_w_info = ROOFS_w_ID[['DF_UID', 'BLD_ID']]
print(ROOFS_w_info.head())
ROOFS_w_info.to_csv( os.path.join(OUT_FOLDER, 'SOLKAT_DACH_BLD_ID_batch%d.csv' %BATCH_ID), index = False )
print('Saved %s' %(os.path.join(OUT_FOLDER, 'SOLKAT_DACH_BLD_ID_batch%d.csv' %BATCH_ID)))
# ## CHECKSUM for BLD_ID
print('Lengths of BLD data: TLM_BLD.BLD_ID.unique(), TLM_BLD, intersection of BLD_ID of ROOFS_w_ID and TLM_w_ID')
print(len(TLM_BLD.BLD_ID.unique()))
print(len(TLM_BLD))
print(len(list(set(TLM_w_ID.BLD_ID.unique()) | set(ROOFS_w_ID.BLD_ID.unique()))))
# ## Select only the core BLD geometries
tt = time.time()
BLD_centr = TLM_BLD.copy()
BLD_centr['geometry'] = TLM_BLD.centroid
BLD_centr['x'] = BLD_centr.geometry.x
BLD_centr['y'] = BLD_centr.geometry.y
# Match the batchID on the BLD_ID
coord_names = ['x_px', 'y_px']
BLD_centr[coord_names] = map_coords_to_pixel(BLD_centr[['x','y']], BATCH_DELTA/2, BATCH_DELTA )
BLD_centr = BLD_centr.merge( batches.drop(columns = ['geometry']),
on = coord_names, how = 'left').drop(columns = coord_names)[['BLD_ID', 'batchID']]
print('Created and mapped BLD centroids in %.2fs' %(time.time()-tt))
## ============ SAVE OUTPUT
# ## TLM
TLM_ID_info = pd.merge(TLM_w_ID.drop(columns = ['geometry']), BLD_centr, on = 'BLD_ID')
TLM_ID_info['BLD_ID'] = TLM_ID_info['BLD_ID'] + TLM_ID_info['batchID'] * ID_OFFSET
print(TLM_ID_info.head())
TLM_ID_info.to_csv( os.path.join(OUT_FOLDER, 'TLM_BLD_ID_batch%d.csv' %BATCH_ID) , index = False)
print('Saved %s' %os.path.join(OUT_FOLDER, 'TLM_BLD_ID_batch%d.csv' %BATCH_ID))
# ## Rooftops
ROOFS_w_info = pd.merge(ROOFS_w_ID[['DF_UID', 'BLD_ID']], BLD_centr, on = 'BLD_ID')
ROOFS_w_info['BLD_ID'] = ROOFS_w_info['BLD_ID'] + ROOFS_w_info['batchID'] * ID_OFFSET
print(ROOFS_w_info.head())
tt = time.time()
ROOFS_w_info.to_csv( os.path.join(OUT_FOLDER, 'SOLKAT_DACH_BLD_ID_batch%d.csv' %BATCH_ID), index = False )
print('Saved %s in %.2fs' %(os.path.join(OUT_FOLDER, 'SOLKAT_DACH_BLD_ID_batch%d.csv' %BATCH_ID), time.time()-tt))
# ## New building polygons
BLD_out = pd.merge(TLM_BLD, BLD_centr, on = 'BLD_ID')
BLD_out['BLD_ID'] = BLD_out['BLD_ID'] + BLD_out['batchID'] * ID_OFFSET
BLD_out.crs = TLM.crs
print(BLD_out.head())
tt = time.time()
BLD_out.to_file( os.path.join(OUT_FOLDER, 'BLD_LV95_batch%d' %BATCH_ID) )
print('Saved %s in %.2fs' %(os.path.join(OUT_FOLDER, 'BLD_LV95_batch%d' %BATCH_ID), time.time()-tt))
c4science · Help