Coercing UTM conversions to same zone

cityseer/util/graphs.py

@@ -98,13 +98,17 @@ def nX_simple_geoms(networkX_graph: nx.Graph) -> nx.Graph:
     return g_copy
 
 
-def nX_wgs_to_utm(networkX_graph: nx.Graph) -> nx.Graph:
+def nX_wgs_to_utm(networkX_graph: nx.Graph, force_zone_number=None) -> nx.Graph:
     if not isinstance(networkX_graph, nx.Graph):
         raise TypeError('This method requires an undirected networkX graph.')
 
     logger.info('Converting networkX graph from WGS to UTM.')
     g_copy = networkX_graph.copy()
 
+    zone_number = None
+    if force_zone_number is not None:
+        zone_number = force_zone_number
+
     logger.info('Processing node x, y coordinates.')
     for n, d in tqdm(g_copy.nodes(data=True), disable=tqdm_suppress):
         # x coordinate
@@ -118,9 +122,12 @@ def nX_wgs_to_utm(networkX_graph: nx.Graph) -> nx.Graph:
         # check for unintentional use of conversion
         if abs(lng) > 180 or abs(lat) > 90:
             raise AttributeError('x, y coordinates exceed WGS bounds. Please check your coordinate system.')
+        # to avoid issues across UTM boundaries, use the first point to set (and subsequently force) the UTM zone
+        if zone_number is None:
+            zone_number = utm.from_latlon(lat, lng)[2]  # zone number is position 2
         # be cognisant of parameter and return order
         # returns in easting, northing order
-        easting, northing = utm.from_latlon(lat, lng)[:2]
+        easting, northing = utm.from_latlon(lat, lng, force_zone_number=zone_number)[:2]
         # write back to graph
         g_copy.nodes[n]['x'] = easting
         g_copy.nodes[n]['y'] = northing
@@ -135,7 +142,7 @@ def nX_wgs_to_utm(networkX_graph: nx.Graph) -> nx.Graph:
                 raise TypeError(f'Expecting LineString geometry but found {line_geom.type} geometry.')
             # be cognisant of parameter and return order
             # returns in easting, northing order
-            utm_coords = [utm.from_latlon(lat, lng)[:2] for lng, lat in line_geom.coords]
+            utm_coords = [utm.from_latlon(lat, lng, force_zone_number=zone_number)[:2] for lng, lat in line_geom.coords]
             # write back to edge
             g_copy[s][e]['geom'] = geometry.LineString(utm_coords)
 

docs/util/graphs.md

@@ -34,9 +34,9 @@ Returns a `networkX` graph with `shapely` [`Linestring`](https://shapely.readthe
 nX\_wgs\_to\_utm
 ----------------
 
-<FuncSignature>nX_wgs_to_utm(networkX_graph)</FuncSignature>
+<FuncSignature>nX_wgs_to_utm(networkX_graph, force_zone_number=None)</FuncSignature>
 
-Converts `x` and `y` node attributes from [WGS84](https://epsg.io/4326) `lng`, `lat` geographic coordinates to the local UTM projected coordinate system. If edge `geom` attributes are found, the associated `LineString` geometries will also be converted.
+Converts `x` and `y` node attributes from [WGS84](https://epsg.io/4326) `lng`, `lat` geographic coordinates to the local UTM projected coordinate system. If edge `geom` attributes are found, the associated `LineString` geometries will also be converted. The UTM zone derived from the first processed node will be used for the conversion of all other nodes and geometries contained in the graph. This ensures consistent behaviour when a graph spans a UTM boundary.
 
 <FuncHeading>Parameters</FuncHeading>
 
@@ -46,6 +46,12 @@ A `networkX` graph with `x` and `y` node attributes in the WGS84 coordinate syst
 
 </FuncElement>
 
+<FuncElement name="force_zone_number" type="int">
+
+An optional UTM zone number for coercing all conversions to an explicit UTM zone. Use with caution, because mismatched UTM zones may introduce substantial distortions in the results. 
+
+</FuncElement>
+
 <FuncHeading>Returns</FuncHeading>
 <FuncElement name="graph" type="nx.Graph">
 

package.json

@@ -6,7 +6,7 @@
   "dependencies": {},
   "devDependencies": {
     "standard": "12.0.1",
-    "vuepress": "^1.0.0-alpha.44",
+    "vuepress": "^1.0.0-alpha.47",
     "vuepress-plugin-mathjax": "1.2.4",
     "vuepress-theme-cityseer": "0.3.0"
   }

setup.py

@@ -27,7 +27,7 @@
 
 setup (
     name = 'cityseer',
-    version='0.8.0',
+    version='0.8.1',
     packages=['cityseer', 'cityseer.algos', 'cityseer.metrics', 'cityseer.util'],
     description = 'Computational tools for urban analysis',
     url='https://github.com/cityseer/cityseer-api',

tests/util/test_graphs.py

@@ -75,6 +75,42 @@ def test_nX_wgs_to_utm():
     with pytest.raises(AttributeError):
         graphs.nX_wgs_to_utm(G_utm)
 
+    # check that non-matching UTM zones are coerced to the same zone
+    # this scenario spans two UTM zones
+    G_wgs_b = nx.Graph()
+    nodes = [
+        (1, {'x': -0.0005, 'y': 51.572}),
+        (2, {'x': -0.0005, 'y': 51.571}),
+        (3, {'x': 0.0005, 'y': 51.570}),
+        (4, {'x': -0.0005, 'y': 51.569}),
+        (5, {'x': -0.0015, 'y': 51.570})
+    ]
+    G_wgs_b.add_nodes_from(nodes)
+    edges = [
+        (1, 2),
+        (2, 3),
+        (3, 4),
+        (4, 5),
+        (5, 2)
+    ]
+    G_wgs_b.add_edges_from(edges)
+    G_utm_30 = graphs.nX_wgs_to_utm(G_wgs_b)
+    G_utm_30 = graphs.nX_simple_geoms(G_utm_30)
+
+    # if not consistently coerced to UTM zone, the distances from 2-3 and 3-4 will be over 400km
+    for s, e, d in G_utm_30.edges(data=True):
+        assert d['geom'].length < 200
+
+    # check that explicit zones are respectively coerced
+    G_utm_31 = graphs.nX_wgs_to_utm(G_wgs_b, force_zone_number=31)
+    G_utm_31 = graphs.nX_simple_geoms(G_utm_31)
+    for n, d in G_utm_31.nodes(data=True):
+        assert d['x'] != G_utm_30.nodes[n]['x']
+
+    # from cityseer.util import plot
+    # plot.plot_nX(G_wgs_b, labels=True)
+    # plot.plot_nX(G_utm_b, labels=True)
+
 
 def make_messy_graph(G):