diff --git a/brails/InventoryGenerator.py b/brails/InventoryGenerator.py index cc07cf2..fc2e337 100644 --- a/brails/InventoryGenerator.py +++ b/brails/InventoryGenerator.py @@ -40,7 +40,7 @@ # Satish Rao # # Last updated: -# 01-29-2024 +# 02-01-2024 import random import sys @@ -58,6 +58,7 @@ YearBuiltClassifier) from .workflow.ImHandler import ImageHandler from .workflow.FootprintHandler import FootprintHandler +from brails.workflow.NSIParser import NSIParser class InventoryGenerator: @@ -136,7 +137,8 @@ def enabled_attributes(self): "InventoryGenerator.generate, simply set attributes='all'") def generate(self,attributes=['numstories','occupancy','roofshape'], - outFile='inventory.csv', lengthUnit='ft'): + useNSIBaseline= True, outFile='inventory.csv', + lengthUnit='ft'): def write_to_dataframe(df,predictions,column,imtype='street_images'): """ @@ -266,7 +268,7 @@ def write_inventory_output(inventorydf,outFile): # Rearrange the column order of dfout such that the Footprint field is # the last: cols = [col for col in dfout.columns if col!='Footprint'] - new_cols = ['Latitude','Longitude'] + cols[:-2] + ['Footprint'] + new_cols = ['Latitude','Longitude'] + cols + ['Footprint'] dfout = dfout[new_cols] # If the inventory is desired in CSV format, write dfout to a CSV: @@ -305,17 +307,25 @@ def write_inventory_output(inventorydf,outFile): json.dump(geojson, output_file, indent=2) print(f'\nFinal inventory data available in {outFile} in {os.getcwd()}') - # Parse/correct the list of user requested building attributes: self.attributes = parse_attribute_input(attributes, self.enabledAttributes) # Create a list of footprints for easier module calls: footprints = self.inventory['Footprint'].values.tolist() + + if useNSIBaseline: + nsi = NSIParser() + nsi.GenerateBldgInventory(footprints) + attributes_process = set(self.attributes) - set(nsi.attributes) + self.inventory = nsi.inventory.copy(deep=True) + footprints = self.inventory['Footprint'].values.tolist() + else: + attributes_process = self.attributes.copy() # Download the images required for the requested attributes: image_handler = ImageHandler(self.apiKey) - if 'roofshape' in self.attributes: #or 'roofcover' in self.attributes: + if 'roofshape' in attributes_process: #or 'roofcover' in attributes_process: image_handler.GetGoogleSatelliteImage(footprints) imsat = [im for im in image_handler.satellite_images if im is not None] self.inventory['satellite_images'] = image_handler.satellite_images @@ -323,12 +333,12 @@ def write_inventory_output(inventorydf,outFile): streetAttributes = self.enabledAttributes[:] streetAttributes.remove('roofshape') #streetAttributes.remove('roofcover') - if set.intersection(set(streetAttributes),set(self.attributes))!=set(): + if set.intersection(set(streetAttributes),set(attributes_process))!=set(): image_handler.GetGoogleStreetImage(footprints) imstreet = [im for im in image_handler.street_images if im is not None] self.inventory['street_images'] = image_handler.street_images - for attribute in self.attributes: + for attribute in attributes_process: if attribute=='chimney': # Initialize the chimney detector object: @@ -449,9 +459,9 @@ def write_inventory_output(inventorydf,outFile): # Bring the attribute values to the desired length unit: if lengthUnit.lower()=='m': self.inventory['PlanArea'] = self.inventory['PlanArea'].apply(lambda x: x*0.0929) - if 'buildingheight' in self.attributes: + if 'buildingheight' in attributes_process: self.inventory['buildingheight'] = self.inventory['buildingheight'].apply(lambda x: x*0.3048) - if 'roofeaveheight' in self.attributes: + if 'roofeaveheight' in attributes_process: self.inventory['roofeaveheight'] = self.inventory['roofeaveheight'].apply(lambda x: x*0.3048) # Write the genereated inventory in outFile: diff --git a/brails/TranspInventoryGenerator.py b/brails/TranspInventoryGenerator.py index ff32d44..f7234c8 100644 --- a/brails/TranspInventoryGenerator.py +++ b/brails/TranspInventoryGenerator.py @@ -157,8 +157,8 @@ def combineAndFormat_HWY(self, minimumHAZUS=True, connectivity=False, maxRoadLen # Dump to json file with open("hwy_inventory.geojson", "w") as f: json.dump(combinedGeoJSON, f, indent = 2) - print('Combined transportation inventory saved in hwy_inventory.geojson' - 'This file is suitable for R2D use and is available in {os.getcwd()}') + print('Combined transportation inventory saved in hwy_inventory.geojson.' + f' This file is suitable for R2D use and is available in {os.getcwd()}') return # Convert common length units diff --git a/brails/workflow/FootprintHandler.py b/brails/workflow/FootprintHandler.py index 044c12f..bf79f79 100644 --- a/brails/workflow/FootprintHandler.py +++ b/brails/workflow/FootprintHandler.py @@ -37,7 +37,7 @@ # Barbaros Cetiner # # Last updated: -# 01-24-2024 +# 02-01-2024 import math import json @@ -58,7 +58,6 @@ def __init__(self): self.availableDataSources = ['osm','ms','usastr'] self.fpSource = 'osm' self.footprints = [] - self.bldgheights = [] def fetch_footprint_data(self,queryarea,fpSource='osm',attrmap=None): """ @@ -494,7 +493,7 @@ def polygon_area(lats, lons): else: #return in ratio of sphere total area return area - def load_footprint_data(fpfile,attrmapfile,fpSource): + def load_footprint_data(fpfile,fpSource,attrmapfile): """ Function that loads footprint data from a GeoJSON file @@ -544,26 +543,70 @@ def fp_download(bbox,fpSource): elif fpSource=='usastr': footprints = get_usastruct_footprints(bbox) return footprints + + def parse_fp_geojson(data, attrmap, attrkeys, fpfile): + # Create the attribute fields that will be extracted from the + # GeoJSON file: + attributes = {} + for attr in attrmap.values(): + attributes[attr] = [] + + footprints_out = [] + discardedfp_count = 0 + correctedfp_count = 0 + for loc in data: + # If the footprint is a polygon: + if loc['geometry']['type']=='Polygon': + # Read footprint coordinates: + temp_fp = loc['geometry']['coordinates'] - with open(attrmapfile) as f: - lines = [line.rstrip() for line in f] - attrmap = {} - for line in lines: - lout = line.split(':') - if lout[1]!='': - attrmap[lout[1]] = lout[0] - - attributes = {} - for attr in attrmap.values(): - attributes[attr] = [] - - with open(fpfile) as f: - data = json.load(f)['features'] - - attrkeys = list(data[0]['properties'].keys()) - print(attrkeys) - ptdata = all(loc['geometry']['type']=='Point' for loc in data) - if ptdata: + # Check down to two levels deep into extracted JSON + # structure to account for inconsistencies in the + # provided footprint data + if len(temp_fp)>1: + fp = temp_fp[:] + elif len(temp_fp[0])>1: + fp = temp_fp[0][:] + elif len(temp_fp[0][0])>1: + fp = temp_fp[0][0][:] + + # If mutliple polygons are detected for a location, + # take the outermost polygon: + if len(fp)==2: + list_len = [len(i) for i in fp] + fp = fp[list_len.index(max(list_len))] + correctedfp_count+=1 + + # Add the footprint and attributes to the output + # variables + footprints_out.append(fp) + if attrkeys: + for key in attrkeys: + try: + attributes[attrmap[key]].append(loc['properties'][key]) + except: + pass + # If the footprint is a multi-polygon, discard the footprint: + elif loc['geometry']['type']=='MultiPolygon': + discardedfp_count+=1 + + # Print the results of the footprint extraction: + if discardedfp_count==0: + print(f"Extracted a total of {len(footprints_out)} building footprints from {fpfile}") + else: + print(f"Corrected {correctedfp_count} building footprint{pluralsuffix(correctedfp_count)} with invalid geometry") + print(f"Discarded {discardedfp_count} building footprint{pluralsuffix(discardedfp_count)} with invalid geometry") + print(f"Extracted a total of {len(footprints_out)} building footprints from {fpfile}") + + return (footprints_out, attributes) + + def parse_pt_geojson(data, attrmap, attrkeys, fpSource): + # Create the attribute fields that will be extracted from the + # GeoJSON file: + attributes = {} + for attr in attrmap.values(): + attributes[attr] = [] + # Write the data in datalist into a dictionary for better data access, # and filtering the duplicate entries: datadict = {} @@ -575,9 +618,21 @@ def fp_download(bbox,fpSource): datadict[pt] = loc['properties'] points = list(datadict.keys()) - # Determine the coordinates of the bounding box including the points: + + # Determine the coordinates of the bounding box containing the + # points: bbox = get_bbox(ptcoords) - footprints = fp_download(bbox,fpSource) + + # Get the footprint data corresponding to the point GeoJSON + # input: + if 'geojson' in fpSource.lower(): + with open(fpSource) as f: + data = json.load(f)['features'] + (footprints,_) = parse_fp_geojson(data, {}, {}, fpSource) + else: + footprints = fp_download(bbox,fpSource) + + # Create an STR tree for efficient parsing of point coordinates: pttree = STRtree(points) # Find the data points that are enclosed in each footprint: @@ -601,44 +656,64 @@ def fp_download(bbox,fpSource): attributes[attrmap[key]].append(ptres[key]) except: pass - + + return (footprints_out, attributes) + + if attrmapfile: + # Create a dictionary for mapping the attributes in the GeoJSON + # file to BRAILS inventory naming conventions: + with open(attrmapfile) as f: + lines = [line.rstrip() for line in f] + attrmap = {} + for line in lines: + lout = line.split(':') + if lout[1]!='': + attrmap[lout[1]] = lout[0] + + # Read the GeoJSON file and check if all the data in the file is + # point data: + with open(fpfile) as f: + data = json.load(f)['features'] + ptdata = all(loc['geometry']['type']=='Point' for loc in data) + + # Identify the attribute keys in the GeoJSON file: + attrkeys0 = list(data[0]['properties'].keys()) + if attrkeys0: + print('Building attributes detected in the input GeoJSON: ' + + ', '.join(attrkeys0)) + + # Check if all of the attribute keys in the GeoJSON have + # correspondence in the map. Ignore the keys that do not have + # correspondence: + attrkeys = set() + for key in attrkeys0: + try: + attrmap[key] + attrkeys.add(key) + except: + pass + ignored_Attr = set(attrkeys0) - attrkeys + if ignored_Attr: + print('Attribute mapping does not cover all attributes detected in' + 'the input GeoJSON. Ignoring detected attributes: ' + + ', '.join(ignored_Attr)) + else: + attrmap = {} + attrkeys = {} + + if ptdata: + (footprints_out, attributes) = parse_pt_geojson(data, + attrmap, + attrkeys, + fpSource) else: - footprints_out = [] - discardedfp_count = 0 - correctedfp_count = 0 - for loc in data: - if loc['geometry']['type']=='Polygon': - temp_fp = loc['geometry']['coordinates'] - if len(temp_fp)>1: - fp = temp_fp[:] - elif len(temp_fp[0])>1: - fp = temp_fp[0][:] - elif len(temp_fp[0][0])>1: - fp = temp_fp[0][0][:] - - if len(fp)==2: - list_len = [len(i) for i in fp] - fp = fp[list_len.index(max(list_len))] - correctedfp_count+=1 - - footprints_out.append(fp) - for key in attrkeys: - try: - attributes[attrmap[key]].append(loc['properties'][key]) - except: - pass - - elif loc['geometry']['type']=='MultiPolygon': - discardedfp_count+=1 - - if discardedfp_count==0: - print(f"Extracted a total of {len(footprints_out)} building footprints from {fpfile}") - else: - print(f"Corrected {correctedfp_count} building footprint{pluralsuffix(correctedfp_count)} with invalid geometry") - print(f"Discarded {discardedfp_count} building footprint{pluralsuffix(discardedfp_count)} with invalid geometry") - print(f"Extracted a total of {len(footprints_out)} building footprints from {fpfile}") + (footprints_out, attributes) = parse_fp_geojson(data, + attrmap, + attrkeys, + fpfile) + fpSource = fpfile - return (footprints_out, attributes) + return (footprints_out, attributes, fpSource) def fp_source_selector(self): if self.fpSource=='osm': @@ -655,20 +730,22 @@ def fp_source_selector(self): self.fpSource = fpSource if isinstance(self.queryarea,str): if 'geojson' in queryarea.lower(): - (self.footprints,self.attributes) = load_footprint_data( + (self.footprints,self.attributes,self.fpSource) = load_footprint_data( self.queryarea, - attrmap, - self.fpSource) + self.fpSource, + attrmap) + bldgheights = [] else: - (self.footprints,self.bldgheights) = fp_source_selector(self) + (self.footprints,bldgheights) = fp_source_selector(self) elif isinstance(queryarea,tuple): - (self.footprints,self.bldgheights) = fp_source_selector(self) + (self.footprints,bldgheights) = fp_source_selector(self) elif isinstance(queryarea,list): self.footprints = [] + bldgheights = [] for query in self.queryarea: (fps, bldghts) = fp_source_selector(query) self.footprints.extend(fps) - self.bldgheights.extend(bldghts) + bldgheights.extend(bldghts) else: sys.exit('Incorrect location entry. The location entry must be defined as' + ' 1) a string or a list of strings containing the name(s) of the query areas,' + @@ -677,8 +754,11 @@ def fp_source_selector(self): ' bounding box of interest in (lon1, lat1, lon2, lat2) format.' + ' For defining a bounding box, longitude and latitude values' + ' shall be entered for the vertex pairs of any of the two' + - ' diagonals of the rectangular bounding box.') - + ' diagonals of the rectangular bounding box.') + + if bldgheights!=[]: + self.attributes['buildingheight'] = bldgheights.copy() + self.fpAreas = [] for fp in self.footprints: lons = [] diff --git a/brails/workflow/NSIParser.py b/brails/workflow/NSIParser.py index 4555a65..efc7956 100644 --- a/brails/workflow/NSIParser.py +++ b/brails/workflow/NSIParser.py @@ -37,7 +37,7 @@ # Barbaros Cetiner # # Last updated: -# 01-24-2024 +# 02-01-2024 import requests import numpy as np @@ -50,8 +50,9 @@ class NSIParser: def __init__(self): self.inventory = [] + self.attributes = ['erabuilt','numstories','occupancy','constype'] - def GenerateBldgInventory(self,footprints): + def GenerateBldgInventory(self,footprints,outFile=None): """ Function that reads NSI buildings points and matches the data to a set of footprints. @@ -122,7 +123,8 @@ def get_nbi_data(box): # Download NBI data using the defined retry strategy, read downloaded # GeoJSON data into a list: - print('\nGetting National Structure Inventory (NSI) building points for the building footprints...') + print('\nGetting National Structure Inventory (NSI) building data' + ' for the building footprints...') response = s.get(url) datalist = response.json()['features'] @@ -149,6 +151,7 @@ def get_inv_from_datadict(datadict,footprints): ress.append(np.empty(shape=(0, 0))) # Match NBI data to each footprint: + footprints_out_json = [] footprints_out = [] lat = [] lon = [] @@ -160,7 +163,8 @@ def get_inv_from_datadict(datadict,footprints): occ = [] for ind, fp in enumerate(footprints): if ress[ind].size!=0: - footprints_out.append(('{"type":"Feature","geometry":' + + footprints_out.append(fp) + footprints_out_json.append(('{"type":"Feature","geometry":' + '{"type":"Polygon","coordinates":[' + f"""{fp}""" + ']},"properties":{}}')) @@ -172,8 +176,16 @@ def get_inv_from_datadict(datadict,footprints): nstories.append(ptres['num_story']) yearbuilt.append(ptres['med_yr_blt']) repcost.append(ptres['val_struct']) - strtype.append(ptres['bldgtype']) - occ.append(ptres['occtype']) + bldgtype = ptres['bldgtype'] + '1' + if bldgtype=='M1': + strtype.append('RM1') + else: + strtype.append(bldgtype) + if '-' in ptres['occtype']: + occ.append(ptres['occtype'].split('-')[0]) + else: + occ.append(ptres['occtype']) + # Display the number of footprints that can be matched to NSI points: print(f'Found a total of {len(footprints_out)} building points in' @@ -189,7 +201,7 @@ def get_inv_from_datadict(datadict,footprints): inventory['StructureType'] = strtype inventory['OccupancyClass'] = occ inventory['Footprint'] = footprints_out - return inventory + return (inventory, footprints_out_json) # Determine the coordinates of the bounding box including the footprints: bbox = get_bbox(footprints) @@ -198,8 +210,13 @@ def get_inv_from_datadict(datadict,footprints): datadict = get_nbi_data(bbox) # Create a footprint-merged building inventory from extracted NBI data: - self.inventory = get_inv_from_datadict(datadict,footprints) + (self.inventory, footprints_out_json) = get_inv_from_datadict(datadict,footprints) - # Write the created inventory in R2D-compatible CSV format: - self.inventory.to_csv('inventory.csv', index=True, index_label='id') - print(f'\nFinal inventory data available in {os.getcwd()}/inventory.csv') \ No newline at end of file + # If requested, write the created inventory in R2D-compatible CSV format: + if outFile: + inventory = self.inventory.copy(deep=True) + n = inventory.columns[1] + inventory.drop(n, axis = 1, inplace = True) + inventory[n] = footprints_out_json + inventory.to_csv(outFile, index=True, index_label='id') + print(f'\nFinal inventory data available in {os.getcwd()}/{outFile}') \ No newline at end of file