diff --git a/brails/utils/geoTools.py b/brails/utils/geoTools.py index ef2fc07..337600a 100644 --- a/brails/utils/geoTools.py +++ b/brails/utils/geoTools.py @@ -37,10 +37,11 @@ # Barbaros Cetiner # # Last updated: -# 03-12-2024 +# 03-13-2024 from math import radians, sin, cos, atan2, sqrt from shapely.geometry import Polygon +import matplotlib.pyplot as plt def haversine_dist(p1,p2): """ @@ -108,4 +109,27 @@ def mesh_polygon(polygon, rows:int, cols:int): # a valid cell: if poly.area > 0: rectangles.append(poly.envelope) - return rectangles \ No newline at end of file + return rectangles + +def plot_polygon_cells(bpoly, rectangles, fout=False): + """ + Function that plots the mesh 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 + """ + + if bpoly.geom_type == 'MultiPolygon': + for poly in bpoly.geoms: + plt.plot(*poly.exterior.xy,'k') + else: + plt.plot(*bpoly.exterior.xy,'k') + for rect in rectangles: + try: + plt.plot(*rect.exterior.xy) + except: + pass + if fout: + plt.savefig(fout, dpi=600, bbox_inches="tight") + plt.show() \ No newline at end of file diff --git a/brails/workflow/FootprintHandler.py b/brails/workflow/FootprintHandler.py index 380bddd..5f3768d 100644 --- a/brails/workflow/FootprintHandler.py +++ b/brails/workflow/FootprintHandler.py @@ -37,7 +37,7 @@ # Barbaros Cetiner # # Last updated: -# 03-08-2024 +# 03-13-2024 import math import json @@ -50,6 +50,9 @@ from shapely.geometry import Point, Polygon, LineString, MultiPolygon, box from shapely.ops import linemerge, unary_union, polygonize from shapely.strtree import STRtree +from brails.utils.geoTools import haversine_dist, mesh_polygon, plot_polygon_cells +import concurrent.futures +from requests.adapters import HTTPAdapter, Retry class FootprintHandler: def __init__(self): @@ -420,35 +423,258 @@ def download_ms_tiles(quadkeys,bpoly): print(f"\nFound a total of {len(footprints)} building footprints in {queryarea_printname}") return footprints,attributes - def get_usastruct_footprints(queryarea): - if isinstance(queryarea,tuple): - queryarea_printname = (f"the bounding box: [{queryarea[0]}," - f"{queryarea[1]}, {queryarea[2]}, " - f"{queryarea[3]}]") + def get_usastruct_footprints(queryarea): + def get_usastruct_bldg_counts(bpoly): + # Get the coordinates of the bounding box for input polygon bpoly: + bbox = bpoly.bounds + + # Get the number of buildings in the computed bounding box: + query = ('https://services2.arcgis.com/FiaPA4ga0iQKduv3/ArcGIS/' + + 'rest/services/USA_Structures_View/FeatureServer/0/query?' + + f'geometry={bbox[0]},{bbox[1]},{bbox[2]},{bbox[3]}' + + '&geometryType=esriGeometryEnvelope&inSR=4326' + + '&spatialRel=esriSpatialRelIntersects' + + '&returnCountOnly=true&f=json') + + s = requests.Session() + retries = Retry(total=5, + backoff_factor=0.1, + status_forcelist=[500, 502, 503, 504]) + s.mount('https://', HTTPAdapter(max_retries=retries)) - elif isinstance(queryarea,str): - sys.exit("This option is not yet available for FEMA USA Structures footprint data") + r = s.get(query) + totalbldgs = r.json()['count'] + + return totalbldgs + + def get_polygon_cells(bpoly, totalbldgs = None, nfeaturesInCell=4000, + plotfout=False): + if totalbldgs is None: + # Get the number of buildings in the input polygon bpoly: + totalbldgs = get_usastruct_bldg_counts(bpoly) + + if totalbldgs>nfeaturesInCell: + # Calculate the number of cells required to cover the polygon area with + # 20 percent margin of error: + ncellsRequired = round(1.2*totalbldgs/nfeaturesInCell) + + # Get the coordinates of the bounding box for input polygon bpoly: + bbox = bpoly.bounds + + # Calculate the horizontal and vertical dimensions of the bounding box: + xdist = haversine_dist((bbox[0],bbox[1]),(bbox[2],bbox[1])) + ydist = haversine_dist((bbox[0],bbox[1]),(bbox[0],bbox[3])) + + # Determine the bounding box aspect ratio defined (as a number greater + # than 1) and the long direction of the bounding box: + if xdist>ydist: + bboxAspectRatio = math.ceil(xdist/ydist) + longSide = 1 + else: + bboxAspectRatio = math.ceil(ydist/xdist) + longSide = 2 + + # Calculate the cells required on the short side of the bounding box (n) + # using the relationship ncellsRequired = bboxAspectRatio*n^2: + n = math.ceil(math.sqrt(ncellsRequired/bboxAspectRatio)) + + # Based on the calculated n value determined the number of rows and + # columns of cells required: + if longSide==1: + rows = bboxAspectRatio*n + cols = n + else: + rows = n + cols = bboxAspectRatio*n + + # Determine the coordinates of each cell covering bpoly: + rectangles = mesh_polygon(bpoly, rows, cols) + else: + rectangles = [bpoly.envelope] + # Plot the generated mesh: + if plotfout: + plot_polygon_cells(bpoly, rectangles, plotfout) + + return rectangles - print(f"\nFetching FEMA USA Structures footprint data for {queryarea_printname}...") + def refine_polygon_cells(premCells, nfeaturesInCell=4000): + # Download the building count for each cell: + pbar = tqdm(total=len(premCells), desc='Obtaining the number of buildings in each cell') + results = {} + with concurrent.futures.ThreadPoolExecutor() as executor: + future_to_url = { + executor.submit(get_usastruct_bldg_counts, rect): rect + for rect in premCells + } + for future in concurrent.futures.as_completed(future_to_url): + rect = future_to_url[future] + pbar.update(n=1) + try: + results[rect] = future.result() + except Exception as exc: + results[rect] = None + print("%r generated an exception: %s" % (rect, exc)) + + indRemove = [] + cells2split = [] + cellsKeep = premCells.copy() + for ind,rect in enumerate(premCells): + totalbldgs = results[rect] + if totalbldgs is not None: + if totalbldgs==0: + indRemove.append(ind) + elif totalbldgs>nfeaturesInCell: + indRemove.append(ind) + cells2split.append(rect) + + for i in sorted(indRemove, reverse=True): + del cellsKeep[i] + + cellsSplit = [] + for rect in cells2split: + rectangles = get_polygon_cells(rect, totalbldgs=results[rect]) + cellsSplit+=rectangles + + return cellsKeep, cellsSplit - query = ('https://services2.arcgis.com/FiaPA4ga0iQKduv3/ArcGIS/' + - 'rest/services/USA_Structures_View/FeatureServer/0/query?' + - f'geometry={queryarea[0]},{queryarea[1]},{queryarea[2]},{queryarea[3]}' - '&geometryType=esriGeometryEnvelope&inSR=4326&spatialRel=esriSpatialRelIntersects&f=pjson') + def download_ustruct_bldgattr(cell): + rect = cell.bounds + s = requests.Session() + retries = Retry(total=5, + backoff_factor=0.1, + status_forcelist=[500, 502, 503, 504]) + s.mount('https://', HTTPAdapter(max_retries=retries)) + query = ('https://services2.arcgis.com/FiaPA4ga0iQKduv3/ArcGIS/' + + 'rest/services/USA_Structures_View/FeatureServer/0/query?' + + f'geometry={rect[0]},{rect[1]},{rect[2]},{rect[3]}' + + '&outFields=BUILD_ID,HEIGHT'+ + '&geometryType=esriGeometryEnvelope&inSR=4326' + + '&spatialRel=esriSpatialRelIntersects&outSR=4326&f=json') + + r = s.get(query) + datalist = r.json()['features'] + ids = [] + footprints = [] + bldgheight = [] + for data in datalist: + footprint = data['geometry']['rings'][0] + bldgid = data['attributes']['BUILD_ID'] + if bldgid not in ids: + ids.append(bldgid) + footprints.append(footprint) + height = data['attributes']['HEIGHT'] + try: + height = round(float(height*3.28084),1) + except: + height = None + bldgheight.append(height) + + return (ids, footprints, bldgheight) - r = requests.get(query) - if 'error' in r.text: - sys.exit("Data server is currently unresponsive." + - " Please try again later.") + def download_ustruct_bldgattr4region(cellsFinal,bpoly): + # Download building attribute data for each cell: + pbar = tqdm(total=len(cellsFinal), desc='Obtaining the building attributes for each cell') + results = {} + with concurrent.futures.ThreadPoolExecutor() as executor: + future_to_url = { + executor.submit(download_ustruct_bldgattr, cell): cell + for cell in cellsFinal + } + for future in concurrent.futures.as_completed(future_to_url): + cell = future_to_url[future] + pbar.update(n=1) + try: + results[cell] = future.result() + except Exception as exc: + results[cell] = None + print("%r generated an exception: %s" % (cell, exc)) + pbar.close() + + # Parse the API results into building id, footprints and height information: + ids = [] + footprints = [] + bldgheight = [] + for cell in tqdm(cellsFinal): + res = results[cell] + ids+=res[0] + footprints+=res[1] + bldgheight+=res[2] - datalist = r.json()['features'] - footprints = [] - for data in datalist: - footprint = data['geometry']['rings'][0] - footprints.append(footprint) + # Remove the duplicate footprint data by recording the API outputs to a + # dictionary: + data = {} + for ind,bldgid in enumerate(ids): + data[bldgid] = [footprints[ind],bldgheight[ind]] + + # Calculate building centroids and save the API outputs into their + # corresponding variables: + footprints = [] + attributes = {'buildingheight':[]} + centroids = [] + for value in data.values(): + fp = value[0] + centroids.append(Polygon(fp).centroid) + footprints.append(fp) + attributes['buildingheight'].append(value[1]) + + # Identify building centroids and that fall outside of bpoly: + results = {} + with concurrent.futures.ThreadPoolExecutor() as executor: + future_to_url = { + executor.submit(bpoly.contains, cent): cent + for cent in centroids + } + for future in concurrent.futures.as_completed(future_to_url): + cent = future_to_url[future] + try: + results[cent] = future.result() + except Exception as exc: + results[cell] = None + print("%r generated an exception: %s" % (cent, exc)) + indRemove = [] + for ind, cent in enumerate(centroids): + if not results[cent]: + indRemove.append(ind) + + # Remove data corresponding to centroids that fall outside bpoly: + for i in sorted(indRemove, reverse=True): + del footprints[i] + del attributes['buildingheight'][i] + + return footprints, attributes + + if isinstance(queryarea,tuple): + bpoly, queryarea_printname = self.__bbox2poly(queryarea) + elif isinstance(queryarea,str): + bpoly, queryarea_printname, _ = self.__fetch_roi(queryarea) + + plotCells = True + + if plotCells: + meshInitialfout = queryarea_printname.replace(' ','_') + '_Mesh_Initial.png' + meshFinalfout = queryarea_printname.replace(' ','_') + '_Mesh_Final.png' - print(f"\nFound a total of {len(footprints)} building footprints in {queryarea_printname}") - return footprints + print('\nMeshing the defined area...') + cellsPrem = get_polygon_cells(bpoly, plotfout=meshInitialfout) + + if len(cellsPrem)>1: + cellsFinal = [] + cellsSplit = cellsPrem.copy() + while len(cellsSplit)!=0: + cellsKeep, cellsSplit = refine_polygon_cells(cellsSplit) + cellsFinal+=cellsKeep + print(f'\nMeshing complete. Split {queryarea_printname} into {len(cellsFinal)} cells') + else: + cellsFinal = cellsPrem.copy() + print(f'\nMeshing complete. Covered {queryarea_printname} with a rectangular cell') + + if plotCells: + plot_polygon_cells(bpoly, cellsFinal, meshFinalfout) + + footprints, attributes = download_ustruct_bldgattr4region(cellsFinal,bpoly) + print(f"\nFound a total of {len(footprints)} building footprints in {queryarea_printname}") + + return footprints, attributes def load_footprint_data(fpfile,fpSource,attrmapfile): """ @@ -718,10 +944,12 @@ def fp_source_selector(self): if self.fpSource=='osm': footprints, attributes = get_osm_footprints(self.queryarea) elif self.fpSource=='ms': - footprints,attributes = get_ms_footprints(self.queryarea) + footprints, attributes = get_ms_footprints(self.queryarea) elif self.fpSource=='usastr': - footprints = get_usastruct_footprints(self.queryarea) - attributes = {} + footprints, attributes = get_usastruct_footprints(self.queryarea) + else: + print('Unimplemented footprint source. Setting footprint source to OSM') + footprints, attributes = get_osm_footprints(self.queryarea) return footprints, attributes self.queryarea = queryarea @@ -764,18 +992,19 @@ def fp_source_selector(self): # Calculate centroids of the footprints and remove footprint data that # does not form a polygon: + self.footprints = [] self.centroids = [] - ind_remove = [] + indRemove = [] for (ind,footprint) in enumerate(footprints): try: self.centroids.append(Polygon(footprint).centroid) self.footprints.append(footprint) except: - ind_remove.append(ind) + indRemove.append(ind) pass # Remove attribute corresponding to the removed footprints: - for i in sorted(ind_remove, reverse=True): + for i in sorted(indRemove, reverse=True): for key in self.attributes.keys(): del self.attributes[key][i] diff --git a/docs/conf.py b/docs/conf.py index b8ea780..de27b11 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -21,7 +21,7 @@ .. |s3harkName| replace:: SURF .. |surfName| replace:: SURF .. |brailsName| replace:: BRAILS -.. |full tool name| replace:: Building Recognition using AI at Large-Scale +.. |full tool name| replace:: Building and Infrastructure Recognition using AI at Large-Scale .. _MessageBoard: https://simcenter-messageboard.designsafe-ci.org/smf/index.php?board=10.0 .. |messageBoard| replace:: `MessageBoard`_ .. |short tool name| replace:: BRAILS