Merge pull request #191 from bacetiner/master

USA Structures downloads are now auto-grided to stay under ESRI API download limits

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
+    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()

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]
 

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