From 1cce31389a487c0efe10faad63f3a2488fa43411 Mon Sep 17 00:00:00 2001 From: bacetiner Date: Tue, 26 Mar 2024 18:32:46 -0700 Subject: [PATCH 1/3] Added STRTree search function for matching points to polygons --- brails/utils/geoTools.py | 67 ++++++++++++++++++++++++++++++---------- 1 file changed, 51 insertions(+), 16 deletions(-) diff --git a/brails/utils/geoTools.py b/brails/utils/geoTools.py index 1dfd9d1..e648c79 100644 --- a/brails/utils/geoTools.py +++ b/brails/utils/geoTools.py @@ -37,20 +37,21 @@ # Barbaros Cetiner # # Last updated: -# 03-15-2024 +# 03-26-2024 from math import radians, sin, cos, atan2, sqrt -from shapely import to_geojson from shapely.geometry import Polygon +from shapely import to_geojson +from shapely.strtree import STRtree import matplotlib.pyplot as plt import json -def haversine_dist(p1,p2): +def haversine_dist(p1:list,p2:list)->float: """ Function that calculates the Haversine distance between two points. - Input: Two points with coordinates defined as latitude and longitude - with each point defined as a list of two floating-point values + Input: Two points with coordinates defined in latitude and longitude + order with each point defined as a list of two floating-point values Output: Distance between the input points in feet """ @@ -74,14 +75,14 @@ def haversine_dist(p1,p2): distance = R*c*3280.84 return distance -def mesh_polygon(polygon, rows:int, cols:int): +def mesh_polygon(polygon,rows:int,cols:int)->list: """ Function that splits a polygon into individual rectangular polygons. Inputs: A Shapely polygon and number of rows and columns of rectangular cells requested to mesh the area of the polygon. - with each point defined as a list of two floating-point values - Output: Individual rectangular cells stored as Shapely polygons + Output: A list containing the individual rectangular cells stored as + Shapely polygons """ # Get bounds of the polygon: @@ -113,13 +114,13 @@ def mesh_polygon(polygon, rows:int, cols:int): rectangles.append(poly.envelope) return rectangles -def plot_polygon_cells(bpoly, rectangles, fout=False): +def plot_polygon_cells(bpoly,rectangles:list,fout=False)->None: """ - Function that plots the mesh for a polygon and saves the plot into a PNG image + Function that plots meshes for a polygon and saves the plot into a PNG image Inputs: A Shapely polygon and a list of rectangular mesh cells saved as - Shapely polygons. Optional input is a string containing the name of the - output file + Shapely polygons. Optional input is a string containing the name of + the output file """ if bpoly.geom_type == 'MultiPolygon': @@ -136,11 +137,12 @@ def plot_polygon_cells(bpoly, rectangles, fout=False): plt.savefig(fout, dpi=600, bbox_inches="tight") plt.show() -def write_polygon2geojson(poly,outfile): +def write_polygon2geojson(poly,outfile:str)->None: """ - Function that writes a single Shapely polygon into a GeoJSON file + Function that writes a single Shapely polygon into a GeoJSON file. - Input: A Shapely polygon or multi polygon + Inputs: A Shapely polygon or multi polygon and a string containing the name + of the GeoJSON file """ if 'geojson' not in outfile.lower(): @@ -162,4 +164,37 @@ def write_polygon2geojson(poly,outfile): to_geojson(poly).split('"coordinates":')[-1][:-1]) geojson['features'].append(feature) with open(outfile, 'w') as outputFile: - json.dump(geojson, outputFile, indent=2) \ No newline at end of file + json.dump(geojson, outputFile, indent=2) + +def match_points2polygons(points:list,polygons:list)->list: + """ + Function that finds the set of points that match a set of polygons + + Inputs: A list of Shapely points and a list of footprint data defined as a + list of lists of coordinates in EPSG 4326, i.e., [[vert1],.... + [vertlast]]. Vertices are defined in [longitude,latitude] fashion. + Outputs: A list of Shapely points and a dictionary that maps each footprint + list of coordinates (converted to string) to the first matched + Shapely point + """ + + # Create an STR tree for the input points: + pttree = STRtree(points) + + # Find the data points that are enclosed in each polygon: + ptkeepind = [] + fp2ptmap = {} + for poly in polygons: + res = pttree.query(Polygon(poly)) + if res.size!=0: + ptkeepind.extend(res) + fp2ptmap[str(poly)] = points[res[0]] + ptkeepind = set(ptkeepind) + + # Create a list of points that include just the points that have a polygon + # match: + ptskeep = [] + for ind in ptkeepind: + ptskeep.append(points[ind]) + + return ptskeep,fp2ptmap \ No newline at end of file From 8de73f9fefa5492ade06f689f654ba120347048e Mon Sep 17 00:00:00 2001 From: bacetiner Date: Tue, 26 Mar 2024 18:34:13 -0700 Subject: [PATCH 2/3] Streamlined NSI attribute extraction and added length unit option --- brails/workflow/NSIParser.py | 230 +++++++++++++++++++---------------- 1 file changed, 124 insertions(+), 106 deletions(-) diff --git a/brails/workflow/NSIParser.py b/brails/workflow/NSIParser.py index 04ae1b6..75e9342 100644 --- a/brails/workflow/NSIParser.py +++ b/brails/workflow/NSIParser.py @@ -37,26 +37,37 @@ # Barbaros Cetiner # # Last updated: -# 03-25-2024 +# 03-26-2024 import requests -import numpy as np import pandas as pd import os from requests.adapters import HTTPAdapter, Retry -from shapely.strtree import STRtree -from shapely import Polygon, Point import sys import json +from shapely.geometry import Point +from brails.utils.geoTools import match_points2polygons +from brails.workflow.FootprintHandler import FootprintHandler class NSIParser: def __init__(self): - self.inventory = [] - self.attributes = ['erabuilt','numstories','occupancy','constype'] + self.attributes = {} + self.brails2r2dmap = {'lon':'Longitude','lat':'Latitude', + 'fparea':'PlanArea', + 'numstories':'NumberOfStories', + 'erabuilt':'YearBuilt', + 'repaircost':'ReplacementCost', + 'constype':'StructureType', + 'occupancy':'OccupancyClass','fp':'Footprint'} + self.footprints = {} + self.nsi2brailsmap = {'x':'lon','y':'lat','sqft':'fparea', + 'num_story':'numstories','med_yr_blt':'erabuilt', + 'val_struct':'repaircost','bldgtype':'constype', + 'occtype':'occupancy'} def __get_bbox(self,footprints:list)->tuple: """ - Function that determines the extent of the area covered by the + Method that determines the extent of the area covered by the footprints as a tight-fit rectangular bounding box Input: List of footprint data defined as a list of lists of @@ -84,7 +95,8 @@ def __get_bbox(self,footprints:list)->tuple: def __get_nbi_data(self,bbox:tuple)->dict: """ - Function that gets the NBI data for a bounding box entry + Method that gets the NBI data for a bounding box entry + Input: Tuple containing the minimum and maximum latitude and longitude values Output: Dictionary containing extracted NBI data keyed using the @@ -126,7 +138,7 @@ def __get_nbi_data(self,bbox:tuple)->dict: def GetRawDataROI(self,roi,outfile:str)->None: """ - The method that reads NSI buildings data finds the points within a + Method that reads NSI buildings data finds the points within a bounding polygon or matches a set of footprints, then writes the extracted data into a GeoJSON file. @@ -163,29 +175,23 @@ def GetRawDataROI(self,roi,outfile:str)->None: # Get NBI data for the computed bounding box: datadict = self.__get_nbi_data(bbox) - pointsKeep = set() points = list(datadict.keys()) - if bpoly is None: - # Create an STR tree for the building points obtained from NBI: - - pttree = STRtree(points) - - # Find the data points that are enclosed in each footprint: - ptind = [] - for fp in footprints: - res = pttree.query(Polygon(fp)) - ptind.extend(res) - ptind = set(ptind) - - for ind in ptind: - pointsKeep.add(points[ind]) - else: + if bpoly is None: + # Find the NBI points that are enclosed in each footprint: + pointsKeep, _ = match_points2polygons(points,footprints) + else: + # Find the NBI points contained in bpoly: + pointsKeep = set() for pt in points: if bpoly.contains(pt): pointsKeep.add(pt) + pointsKeep = list(pointsKeep) - pointsKeep = list(pointsKeep) + # Display the number of NSI points that are within roi: + print(f'Found a total of {len(pointsKeep)} building points in' + ' NSI that are within the entered region of interest') + # Write the extracted NSI data in a GeoJSON file: attrkeys = list(datadict[pointsKeep[0]].keys()) geojson = {'type':'FeatureCollection', "crs": {"type": "name", "properties": @@ -204,93 +210,76 @@ def GetRawDataROI(self,roi,outfile:str)->None: geojson['features'].append(feature) with open(outfile, 'w') as outputFile: - json.dump(geojson, outputFile, indent=2) - + json.dump(geojson, outputFile, indent=2) - def GenerateBldgInventory(self,footprints,bpoly=None,outFile=None): + def GetNSIData(self,footprints:list,lengthUnit:str='ft',outfile=None): """ Method that reads NSI buildings points and matches the data to a set - of footprints. + of footprints and writes the data in a BRAILS-compatible format - Input: List of footprint data defined as a list of lists of + Input: roi: List of footprint data defined as a list of lists of coordinates in EPSG 4326, i.e., [[vert1],....[vertlast]]. Vertices are defined in [longitude,latitude] fashion. - Output: Building inventory for the footprints containing the attributes + + outfile: string containing the name of the output file + Currently, only GeoJSON format is supported + + Output: Attribute dictionary for the footprints containing the attributes: 1)latitude, 2)longitude, 3)building plan area, 4)number of floors, 5)year of construction, 6)replacement cost, 7)structure type, 8)occupancy class, 9)footprint polygon """ - - - - def get_inv_from_datadict(datadict,footprints): - # Create an STR tree for the building points obtained from NBI: + def get_attr_from_datadict(datadict,footprints,nsi2brailsmap): + # Parsers for building and occupancy types: + def bldgtype_parser(bldgtype): + bldgtype = bldgtype + '1' + if bldgtype=='M1': + bldgtype = 'RM1' + return bldgtype + def occ_parser(occtype): + if '-' in occtype: + occtype = occtype.split('-')[0] + return occtype + + # Extract the NSI points from the data dictionary input: points = list(datadict.keys()) - pttree = STRtree(points) - # Find the data points that are enclosed in each footprint: - ress = [] - for fp in footprints: - res = pttree.query(Polygon(fp)) - ress.extend(res) - if res.size!=0: - ress.append(res) - else: - ress.append(np.empty(shape=(0, 0))) + # Match NBI data to footprints: + _, fp2ptmap = match_points2polygons(points,footprints) + + # Write the matched footprints to a list: + footprintsMatched = fp2ptmap.keys() - # Match NBI data to each footprint: - footprints_out_json = [] - footprints_out = [] - lat = [] - lon = [] - planarea = [] - nstories = [] - yearbuilt = [] - repcost = [] - strtype = [] - occ = [] - for ind, fp in enumerate(footprints): - if ress[ind].size!=0: - footprints_out.append(fp) - footprints_out_json.append(('{"type":"Feature","geometry":' + - '{"type":"Polygon","coordinates":[' + - f"""{fp}""" + - ']},"properties":{}}')) - ptind = ress[ind][0] - ptres = datadict[points[ptind]] - lat.append(ptres['y']) - lon.append(ptres['x']) - planarea.append(ptres['sqft']) - nstories.append(ptres['num_story']) - yearbuilt.append(ptres['med_yr_blt']) - repcost.append(ptres['val_struct']) - 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']) - + # Initialize the attributes dictionary: + attributes = {'fp':[], 'fp_json':[]} + nsikeys = list(nsi2brailsmap.keys()) + brailskeys = list(nsi2brailsmap.values()) + for key in brailskeys: + attributes[key] = [] + + # Extract NSI attributes corresponding to each building footprint + # and write them in attributes dictionary: + for fpstr in footprintsMatched: + fp = json.loads(fpstr) + pt = fp2ptmap[fpstr] + ptres = datadict[pt] + attributes['fp'].append(fp) + attributes['fp_json'].append(('{"type":"Feature","geometry":' + + '{"type":"Polygon","coordinates":[' + + f"""{fp}""" + ']},"properties":{}}')) + for key in nsikeys: + if key=='bldgtype': + attributes[nsi2brailsmap[key]].append(bldgtype_parser(ptres[key])) + elif key=='occtype': + attributes[nsi2brailsmap[key]].append(occ_parser(ptres[key])) + else: + attributes[nsi2brailsmap[key]].append(ptres[key]) # Display the number of footprints that can be matched to NSI points: - print(f'Found a total of {len(footprints_out)} building points in' - ' NSI that match the footprint data.') - - # Write the extracted features into a Pandas dataframe: - inventory = pd.DataFrame(pd.Series(lat,name='Latitude')) - inventory['Longitude'] = lon - inventory['PlanArea'] = planarea - inventory['NumberOfStories'] = nstories - inventory['YearBuilt'] = yearbuilt - inventory['ReplacementCost'] = repcost - inventory['StructureType'] = strtype - inventory['OccupancyClass'] = occ - inventory['Footprint'] = footprints_out - return (inventory, footprints_out_json) + print(f'Found a total of {len(footprintsMatched)} building points' + ' in NSI that match the footprint data.') + return attributes # Determine the coordinates of the bounding box including the footprints: bbox = self.__get_bbox(footprints) @@ -299,13 +288,42 @@ def get_inv_from_datadict(datadict,footprints): datadict = self.__get_nbi_data(bbox) # Create a footprint-merged building inventory from extracted NBI data: - (self.inventory, footprints_out_json) = get_inv_from_datadict(datadict,footprints) + attributes = get_attr_from_datadict(datadict,footprints,self.nsi2brailsmap) + + # Convert to footprint areas to sqm if needed: + if lengthUnit=='m': + fparea = [area/(3.28084**2) for area in attributes['fparea']] + attributes.update({'fparea':fparea}) + + + self.attributes = attributes.copy() + if outfile: + if 'csv' in outfile.lower(): + # Write the extracted features into a Pandas dataframe: + brailskeys = list(self.brails2r2dmap.keys()) + inventory = pd.DataFrame(pd.Series(attributes[brailskeys[0]], + name=self.brails2r2dmap[brailskeys[0]])) + for key in brailskeys[1:]: + inventory[self.brails2r2dmap[key]] = attributes[key] + + # Write the created inventory in R2D-compatible CSV format: + n = inventory.columns[-1] + inventory.drop(n, axis = 1, inplace = True) + inventory[n] = attributes['fp_json'] + inventory.to_csv(outfile, index=True, index_label='id') + else: + # Write the extracted features into a GeoJSON file: + footprints = attributes['fp'].copy() + del attributes['fp_json'] + del attributes['fp'] + fpHandler = FootprintHandler() + fpHandler._FootprintHandler__write_fp2geojson(footprints,attributes,outfile) + outfile = outfile.split('.')[0] + '.geojson' + + print(f'\nFinal inventory data available in {os.getcwd()}/{outfile}') - # 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 + # Reformat the class variables to make them compatible with the + # InventoryGenerator: + self.footprints = self.attributes['fp'].copy() + del self.attributes['fp'] + del self.attributes['fp_json'] \ No newline at end of file From 6616b5eec2a0d1ad5ee68e731e461975181e07a9 Mon Sep 17 00:00:00 2001 From: bacetiner Date: Tue, 26 Mar 2024 18:35:21 -0700 Subject: [PATCH 3/3] Changed attribute keys for consistency with other workflow components --- brails/workflow/FootprintHandler.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/brails/workflow/FootprintHandler.py b/brails/workflow/FootprintHandler.py index 66ba051..dd142e8 100644 --- a/brails/workflow/FootprintHandler.py +++ b/brails/workflow/FootprintHandler.py @@ -37,7 +37,7 @@ # Barbaros Cetiner # # Last updated: -# 03-22-2024 +# 03-26-2024 import math import json @@ -1055,14 +1055,14 @@ def fp_source_selector(self): ' shall be entered for the vertex pairs of any of the two' + ' diagonals of the rectangular bounding box.') - self.attributes['fpAreas'] = [] + self.attributes['fparea'] = [] for fp in footprints: lons = [] lats = [] for pt in fp: lons.append(pt[0]) lats.append(pt[1]) - self.attributes['fpAreas'].append(int(polygon_area(lats,lons, + self.attributes['fparea'].append(int(polygon_area(lats,lons, lengthUnit))) # Calculate centroids of the footprints and remove footprint data that