Merge pull request #29 from LocationMind/dev

Update OSM process and documentation for beta test

Data/Results/Automatic_result.md

@@ -1,6 +1,6 @@
 # Quality criterias result between OSM and OMF datasets
 
-The test were run on 14/08/2024 14:22, using the 2024-06-13-beta.1 release of OvertureMap data and the OSM data until 2024/06/07.
+The test were run on 21/08/2024 23:51, using the 2024-06-13-beta.1 release of OvertureMap data and the OSM data until 2024/06/07.
 
 ## General results
 

Documentation/user-doc.md

@@ -133,7 +133,7 @@ python.exe -m pip install --upgrade pip
 **Install dependencies**
 
 ```cmd
-pip install pip-tools && pip-compile Requirements\requirements.txt && pip install -r Requirements\requirements.txt
+pip install pip-tools && pip-compile Requirements\requirements.in && pip install -r Requirements\requirements.txt
 ```
 
 Of course, you have to adapt the path and names if you have changed anything.

Python/osm.py

@@ -101,13 +101,15 @@ def addMissingColumns(connection:psycopg2.extensions.connection,
 def createTableToAggregateEdges(connection:psycopg2.extensions.connection,
                                 edgeTable:str,
                                 area:str,
+                                utmProj:int,
                                 schema:str = 'public'):
     """Create a table to join parallel edges.
 
     Args:
         connection (psycopg2.extensions.connection): Database connection token.
         edgeTable (str): Name of the edge table.
         area (str): Name of the area.
+        utmProj (int): UTM projection for the area.
         schema (str, optional): Name of the schema. Defaults to 'public'.
     """
     sql = f"""
@@ -167,8 +169,8 @@ def createTableToAggregateEdges(connection:psycopg2.extensions.connection,
     FROM {schema}.{edgeTable} AS e1
     LEFT JOIN {schema}.{edgeTable} AS e2 ON e1.u = e2.v AND e1.v = e2.u
     AND e1.id != e2.id AND e1.highway = e2.highway
-    AND ST_Contains(ST_Buffer(ST_Transform(e1.geom, 6691), 0.5), ST_Transform(e2.geom, 6691))
-    AND ST_Contains(ST_Buffer(ST_Transform(e2.geom, 6691), 0.5), ST_Transform(e1.geom, 6691))
+    AND ST_Contains(ST_Buffer(ST_Transform(e1.geom, {utmProj}), 0.5), ST_Transform(e2.geom, {utmProj}))
+    AND ST_Contains(ST_Buffer(ST_Transform(e2.geom, {utmProj}), 0.5), ST_Transform(e1.geom, {utmProj}))
     ORDER BY e1.id;
     
     -- Add cost and reverse cost columns
@@ -186,8 +188,8 @@ def createTableToAggregateEdges(connection:psycopg2.extensions.connection,
     UPDATE {schema}.{area}
     SET reverse_cost = ST_Length(geom1::geography)
     WHERE u1 = v2 AND v1 = u2 AND highway1 = highway2 AND id1 != id2
-    AND ST_Contains(ST_Buffer(ST_Transform(geom1, 6691), 0.5), ST_Transform(geom2, 6691))
-    AND ST_Contains(ST_Buffer(ST_Transform(geom2, 6691), 0.5), ST_Transform(geom1, 6691));
+    AND ST_Contains(ST_Buffer(ST_Transform(geom1, {utmProj}), 0.5), ST_Transform(geom2, {utmProj}))
+    AND ST_Contains(ST_Buffer(ST_Transform(geom2, {utmProj}), 0.5), ST_Transform(geom1, {utmProj}));
     
     CREATE INDEX {area}_geom1_idx
     ON {schema}.{area} USING gist (geom1);
@@ -203,14 +205,16 @@ def createTableToAggregateEdges(connection:psycopg2.extensions.connection,
 
 
 def getBidirectionalRoads(engine:sqlalchemy.engine.base.Engine,
-                           area:str,
-                           schema:str = 'public',
-                           geomColumn:str = 'geom1') -> gpd.GeoDataFrame:
+                          area:str,
+                          utmProj:int,
+                          schema:str = 'public',
+                          geomColumn:str = 'geom1') -> gpd.GeoDataFrame:
     """Get only bidirectional roads from the parallel edge table.
 
     Args:
         engine (sqlalchemy.engine.base.Engine): Engine with the database connection.
         area (str): Name of the area.
+        utmProj (int): UTM projection for the area.
         schema (str, optional): Name of the schema. Defaults to 'public'.
         geomColumn (str, optional): Name of the geomColumn to use. Defaults to 'geom1'.
     
@@ -221,23 +225,25 @@ def getBidirectionalRoads(engine:sqlalchemy.engine.base.Engine,
     sql_bi_roads = f"""
     SELECT * FROM {schema}.{area}
     WHERE u1 = v2 AND v1 = u2 AND highway1 = highway2
-    AND ST_Contains(ST_Buffer(ST_Transform(geom1, 6691), 0.5), ST_Transform(geom2, 6691))
-    AND ST_Contains(ST_Buffer(ST_Transform(geom2, 6691), 0.5), ST_Transform(geom1, 6691));"""
+    AND ST_Contains(ST_Buffer(ST_Transform(geom1, {utmProj}), 0.5), ST_Transform(geom2, {utmProj}))
+    AND ST_Contains(ST_Buffer(ST_Transform(geom2, {utmProj}), 0.5), ST_Transform(geom1, {utmProj}));"""
 
     bi = gpd.read_postgis(sql_bi_roads, engine, geom_col=geomColumn)
 
     return bi
 
 
 def getUnidirectionalRoads(engine:sqlalchemy.engine.base.Engine,
-                            area:str,
-                            schema:str = 'public',
-                            geomColumn:str = 'geom1') -> gpd.GeoDataFrame:
+                           area:str,
+                           utmProj:int,
+                           schema:str = 'public',
+                           geomColumn:str = 'geom1') -> gpd.GeoDataFrame:
     """Get only unidirectional roads from the parallel edge table.
 
     Args:
         engine (sqlalchemy.engine.base.Engine): Engine with the database connection.
         area (str): Name of the area.
+        utmProj (int): UTM projection for the area.
         schema (str, optional): Name of the schema. Defaults to 'public'.
         geomColumn (str, optional): Name of the geomColumn to use. Defaults to 'geom1'.
     
@@ -250,8 +256,8 @@ def getUnidirectionalRoads(engine:sqlalchemy.engine.base.Engine,
     WHERE id1 not in (
         SELECT id1 FROM {schema}.{area}
         WHERE u1 = v2 AND v1 = u2 AND highway1 = highway2
-        AND ST_Contains(ST_Buffer(ST_Transform(geom1, 6691), 0.5), ST_Transform(geom2, 6691))
-        AND ST_Contains(ST_Buffer(ST_Transform(geom2, 6691), 0.5), ST_Transform(geom1, 6691)));"""
+        AND ST_Contains(ST_Buffer(ST_Transform(geom1, {utmProj}), 0.5), ST_Transform(geom2, {utmProj}))
+        AND ST_Contains(ST_Buffer(ST_Transform(geom2, {utmProj}), 0.5), ST_Transform(geom1, {utmProj})));"""
 
     uni = gpd.read_postgis(sql_uni_road, engine, geom_col=geomColumn)
 
@@ -429,6 +435,13 @@ def createGraph(connection:psycopg2.extensions.connection,
     end = time.time()
     log(f"Save edge to postgis: {end - start} seconds")
 
+    # Get utm proj for the area
+    utmProj = utils.getUTMProjFromArea(connection, area)
+
+    end = time.time()
+    log(f"UTM Proj is {utmProj}")
+    log(f"Utm proj: {end - start} seconds")
+
     # Add missing columns if needed to the edge table
     addMissingColumns(connection, edgeTable, schema=schema)
 

Python/utils.py

@@ -1,6 +1,7 @@
 import psycopg2
 import sqlalchemy
 import os
+import utm
 
 def bboxCSVToBboxWKT(bboxCSV:str) -> str:
     """Transform a bounding box in CSV format to its equivalent in OGC WKT format.
@@ -383,7 +384,8 @@ def insertBoundingBox(connection:psycopg2.extensions.connection,
         connection (psycopg2.extensions.connection): Database connection token.
         wktGeom (str): Geometry in OGC WKT format.
         aeraName (str): Name of the area.
-        tableName (str, optional): Name of the table to create. Defaults to 'bounding_box'.
+        tableName (str, optional): Name of the bounding box table.
+        Defaults to 'bounding_box'.
     
     Return:
         int: Id of the inserted row
@@ -495,4 +497,44 @@ def downloadOMFTypeBbox(bbox: str,
 
     print(f"{dataType} data has been downloaded")
 
-    return path
+    return path
+
+
+def getUTMProjFromArea(connection:psycopg2.extensions.connection,
+                       aeraName:str,
+                       tableName:str = 'bounding_box') -> int:
+    """Get UTM projection for a specific point.
+
+    Args:
+        connection (psycopg2.extensions.connection): Database connection token.
+        aeraName (str): Name of the area.
+        tableName (str, optional): Name of the bounding box table.
+        Defaults to 'bounding_box'.
+
+    Returns:
+        int: UTM projection crs id
+    """
+    # Get centroid of the area
+    sql = f"""
+    WITH area_centroid AS (
+        SELECT public.ST_Centroid(geom) as center FROM public.{tableName}
+        WHERE name = '{aeraName.capitalize()}'
+        LIMIT 1
+    )
+    SELECT public.ST_X(ac.center) AS lon, public.ST_Y(ac.center) AS lat
+    FROM area_centroid AS ac;
+    """
+
+    # Execute query
+    cursor = executeSelectQuery(connection, sql)
+
+    # Get lat and lon from row
+    (lat, lon) = cursor.fetchone()
+
+    # Get zone number
+    zone = utm.from_latlon(latitude = lat, longitude = lon)[2]
+
+    # Return formatted string
+    epsgCode = f"326{zone:02d}" if lat >= 0 else f"327{zone:02d}"
+
+    return int(epsgCode)