create_BLD_ID_batch.py
cdf83e4b6b33master

create_BLD_ID_batch.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; 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))

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

			if (len(TLM) == 0) & (len(ROOFS) == 0):
			print('No buildings found')
			exit()

			# # Dissolve and split polygons
			TLM_BLD = dissolve_overlapping_polygons(TLM[['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 by itself 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.overlay( TLM, TLM_BLD, how = 'intersection' )
			print('Overlaid dissolved TLM with original TLM in %.2fs' %(time.time()-tt))

			TLM_remaining = TLM[np.logical_not(TLM.UUID.isin(TLM_w_ID.UUID))]

			if len(TLM_remaining) > 0:
			tt = time.time()
			TLM_w_ID_rem = gpd.sjoin( TLM_remaining, TLM_BLD, op = 'intersects' ).drop(columns = ['index_right'])
			print('Spatially joined remaining TLM in %.2fs' %(time.time()-tt))

			TLM_w_ID = pd.concat([TLM_w_ID, TLM_w_ID_rem])

			# Clean duplicates
			TLM_w_ID['area'] = TLM_w_ID.geometry.area
			TLM_w_ID = TLM_w_ID.sort_values('area', ascending = False)
			TLM_w_ID = TLM_w_ID.drop_duplicates(subset = 'UUID', keep = 'first').drop(columns = ['area'])

			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_BLD = ROOFS.merge( TLM_w_ID[['UUID', 'BLD_ID']].rename({'UUID' : 'SB_UUID'}, axis = 1).dropna(),
			on = 'SB_UUID', how = 'inner')
			print('Merged rooftop ID to BLD in %.2fs' %(time.time()-tt))

			print('Lengths of ROOF subsets: ROOFS_w_BLD, ROOFS, ROOFS_w_BLD.SB_UUID.unique(), ROOFS.SB_UUID.unique()')
			print(len(ROOFS_w_BLD))
			print(len(ROOFS))

			print(len(ROOFS_w_BLD.SB_UUID.unique()))
			print(len(ROOFS.SB_UUID.unique()))


			# ### Check roofs that are not found

			remaining_roofs = ROOFS[np.logical_not(ROOFS.SB_UUID.isin( ROOFS_w_BLD.SB_UUID.values ))]

			if len(remaining_roofs) > 0:
			invalid_geometry = remaining_roofs[remaining_roofs.geometry.is_valid == False ]
			remaining_roofs.loc[invalid_geometry.index, 'geometry'] = invalid_geometry.geometry.buffer(0)

			# ### For roofs that are not found, merge and intersect with TLM_BLD

			tt = time.time()
			SB_UUIDs = remaining_roofs[['geometry', 'SB_UUID']].dissolve(by = 'SB_UUID').reset_index()
			print('Dissolved remaining roofs in %.2fs' %(time.time()-tt))

			tt = time.time()
			SB_UUID_w_ID = gpd.overlay( SB_UUIDs, TLM_BLD, how = 'intersection' )
			print('Overlaid SB_UUID with TLM_BLD in %.2fs' %(time.time()-tt))

			SB_UUID_rem = SB_UUIDs[np.logical_not(SB_UUIDs.SB_UUID.isin(SB_UUID_w_ID.SB_UUID))]

			if len(SB_UUID_rem) > 0:
			tt = time.time()
			SB_UUID_w_ID_rem = gpd.sjoin( SB_UUID_rem, TLM_BLD, op = 'intersects' ).drop(columns = ['index_right'])
			print('Spatially joined remaining SB_UUIDs in %.2fs' %(time.time()-tt))

			SB_UUID_w_ID = pd.concat([SB_UUID_w_ID, SB_UUID_w_ID_rem])

			# Clean duplicates
			SB_UUID_w_ID['area'] = SB_UUID_w_ID.geometry.area
			SB_UUID_w_ID = SB_UUID_w_ID.sort_values('area', ascending = False)
			SB_UUID_w_ID = SB_UUID_w_ID.drop_duplicates(subset = 'SB_UUID', keep = 'first').drop(columns = ['area'])

			# ### Merge information with mapped roofs

			ROOFS_w_BLD = pd.concat([ ROOFS_w_BLD, SB_UUID_w_ID[['SB_UUID', 'BLD_ID']].merge(ROOFS, on = 'SB_UUID') ],
			sort = True)

			print('Lengths of ROOF subsets after cleaning: ROOFS_w_BLD, ROOFS, ROOFS_w_BLD.SB_UUID.unique(), ROOFS.SB_UUID.unique()')
			print(len(ROOFS_w_BLD))
			print(len(ROOFS))

			print(len(ROOFS_w_BLD.SB_UUID.unique()))
			print(len(ROOFS.SB_UUID.unique()))


			# ### Get yet remaining roofs

			remaining_SB = SB_UUIDs[np.logical_not(SB_UUIDs.SB_UUID.isin(ROOFS_w_BLD.SB_UUID.values))]
			# remaining_SB['geometry'] = remaining_SB.geometry.centroid

			# ### Assign new BLD_ID to yet remaning SB_UUIDs

			remaining_SB['BLD_ID'] = np.arange( TLM_BLD.BLD_ID.max() + 1, TLM_BLD.BLD_ID.max() + len(remaining_SB) + 1 )

			print('Added BLD_ID for non-found geometries')
			# print(remaining_SB)

			# ### Merge with roof info

			ROOFS_w_BLD = pd.concat([ ROOFS_w_BLD, ( remaining_SB.drop(columns = ['geometry']).merge(ROOFS, on = 'SB_UUID') )],
			sort = True )

			print('Lengths of remaining ROOF subsets: ROOFS_w_BLD, ROOFS_w_BLD.dropna(), ROOFS, ROOFS_w_BLD.SB_UUID.unique(), ROOFS.SB_UUID.unique()')
			print(len(ROOFS_w_BLD))
			print(len(ROOFS_w_BLD.dropna()))
			print(len(ROOFS))

			print(len(ROOFS_w_BLD.SB_UUID.unique()))
			print(len(ROOFS.SB_UUID.unique()))

			# # Add additional shapes to TLM_BLD

			TLM_BLD = pd.concat([TLM_BLD, remaining_SB[TLM_BLD.columns]])

			else:
			print('No remaining roofs')


			# ## Rooftops
			ROOFS_w_info = ROOFS_w_BLD[['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) )
			print('Saved %s' %(os.path.join(OUT_FOLDER, 'SOLKAT_DACH_BLD_ID_batch%d.csv' %BATCH_ID)))

			# ## 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_BLD[['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) )
			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))

solarPV_identification/PV_postprocessing/Create_BLD/create_BLD_ID_batch.pycdf83e4b6b33master

create_BLD_ID_batch.pyView Options

solarPV_identification/PV_postprocessing/Create_BLD/create_BLD_ID_batch.py
cdf83e4b6b33master

create_BLD_ID_batch.py
View Options