Merge pull request #192 from bacetiner/master

FootprintHandler can now identify small towns and create GeoJSON files for boundary polygons

brails/utils/geoTools.py

@@ -37,11 +37,13 @@
 # Barbaros Cetiner
 #
 # Last updated:
-# 03-13-2024
+# 03-15-2024
 
 from math import radians, sin, cos, atan2, sqrt
+from shapely import to_geojson
 from shapely.geometry import Polygon
 import matplotlib.pyplot as plt
+import json
 
 def haversine_dist(p1,p2):
     """
@@ -132,4 +134,32 @@ def plot_polygon_cells(bpoly, rectangles, fout=False):
             pass        
     if fout:
         plt.savefig(fout, dpi=600, bbox_inches="tight")
-    plt.show()
+    plt.show()
+
+def write_polygon2geojson(poly,outfile):
+    """
+    Function that writes a single Shapely polygon into a GeoJSON file
+    
+    Input: A Shapely polygon or multi polygon
+    """
+
+    if 'geojson' not in outfile.lower():
+        outfile = outfile.replace(outfile.split('.')[-1],'geojson')
+    geojson = {'type':'FeatureCollection', 
+               "crs": {"type": "name", 
+                       "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}},
+               'features':[]}
+    if poly.geom_type == 'MultiPolygon':
+        polytype = 'MultiPolygon'
+    elif poly.geom_type == 'Polygon':
+        polytype = 'Polygon'
+
+    feature = {'type':'Feature',
+               'properties':{},
+               'geometry':{'type': polytype,
+                           'coordinates':[]}}
+    feature['geometry']['coordinates'] = json.loads(
+        to_geojson(poly).split('"coordinates":')[-1][:-1])
+    geojson['features'].append(feature)
+    with open(outfile, 'w') as outputFile:
+        json.dump(geojson, outputFile, indent=2)    

brails/workflow/FootprintHandler.py

@@ -37,7 +37,7 @@
 # Barbaros Cetiner
 #
 # Last updated:
-# 03-13-2024  
+# 03-15-2024  
 
 import math
 import json
@@ -50,7 +50,7 @@
 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
+from brails.utils.geoTools import *
 import concurrent.futures
 from requests.adapters import HTTPAdapter, Retry
 
@@ -62,7 +62,7 @@ def __init__(self):
         self.footprints = []
         self.attributes = {}
 
-    def __fetch_roi(self,queryarea):
+    def __fetch_roi(self,queryarea,outfile=False):
         # Search for the query area using Nominatim API:
         print(f"\nSearching for {queryarea}...")
         queryarea = queryarea.replace(" ", "+").replace(',','+')
@@ -83,13 +83,7 @@ def __fetch_roi(self,queryarea):
         for data in datalist:
             queryarea_osmid = data['osm_id']
             queryarea_name = data['display_name']
-            if(data['osm_type']=='relation' and 
-                'university' in queryarea.lower() and
-                data['type']=='university'):
-                areafound = True
-                break
-            elif (data['osm_type']=='relation' and 
-                  data['type']=='administrative'): 
+            if data['osm_type']=='relation':
                 areafound = True
                 break
 
@@ -136,9 +130,11 @@ def __fetch_roi(self,queryarea):
             sys.exit(f"Could not retrieve the boundary for {queryarea}. " + 
                      'Please check your location query to make sure ' +
                      'it was entered correctly.')    
+        if outfile:
+            write_polygon2geojson(bpoly,outfile)   
         return bpoly, queryarea_printname, queryarea_osmid
 
-    def __bbox2poly(self,queryarea):
+    def __bbox2poly(self,queryarea,outfile=False):
         # Parse the entered bounding box into a polygon:
         if len(queryarea)%2==0 and len(queryarea)!=0:                        
             if len(queryarea)==4:
@@ -158,6 +154,8 @@ def __bbox2poly(self,queryarea):
         else:
                 raise ValueError('Incorrect number of elements detected in the tuple for the bounding box. ' 
                                  'Please check to see if you are missing a longitude or latitude value.')  
+        if outfile:
+            write_polygon2geojson(bpoly,outfile)  
         return bpoly, queryarea_printname
 
     def __write_fp2geojson(self,footprints,attributes,outputFilename):