diff --git a/data/cad/cad.xlsx b/data/cad/cad.xlsx new file mode 100644 index 0000000..dc2788f Binary files /dev/null and b/data/cad/cad.xlsx differ diff --git a/examples/map_cad2geojson.py b/examples/map_cad2geojson.py new file mode 100644 index 0000000..8336fa2 --- /dev/null +++ b/examples/map_cad2geojson.py @@ -0,0 +1,187 @@ +import logging +import math +import os +from collections import Counter +from typing import Callable, Dict, List, Set, Optional + +import numpy as np +import pandas as pd +from scipy.spatial import KDTree + +from mosstool.map.builder import Builder +from mosstool.map.osm import RoadNet +from mosstool.type import Map +from mosstool.util.format_converter import dict2pb + + +# The following will outline the brief steps for importing from CAD to xlsx files, using AutoCAD 2025 - 简体中文 (Simplified Chinese) as an example. +# 1. Navigate to and click on "插入-链接" and "提取-提取数据". In the "界面-“数据提取-开始", click "创建新数据提取" then proceed to next. +# 2. After clicking "next", a prompt appears to save the data extraction as a "*.dxe" file. Name it and click "save". +# 3. After saving, the interface "界面-“数据提取-定义数据源" opens. Click on "数据源-图形/图纸集" and check "包括当前图形". Then click "next". +# 4. In the "数据提取-选取对象" interface, uncheck "显示所有对象类型", click "仅显示块", then click "next". +# 5. In the "数据提取-选择特性" interface, keep default selections and click "next". +# 6. In the "数据提取-优化数据" interface, keep default settings and click "next". +# 7. In the "数据提取-选择输出" interface, check "将数据提取处理表插入图形" and "将数据输入至外部文件". Click "...", choose ".xlsx" type in "另存为" dialog, input save path and filename, click "save", then click "next" in "选择输出". +# 8. In the "数据提取-表格样式" interface, keep default settings and click "next". +# 9. In the "数据提取-完成" interface, click finish to complete the attribute extraction and get the output Excel file. +def cad2osm( + cad_path: str, + node_start_id: int = 0, + way_start_id: int = 0, + merge_gate: float = 0.0, + x_transform_func: Optional[Callable[[float], float]] = None, + y_transform_func: Optional[Callable[[float], float]] = None, +) -> List[Dict]: + df = pd.read_excel(cad_path) + df = df.fillna(0) + if x_transform_func is None: + x_transform_func = lambda x: x + if y_transform_func is None: + y_transform_func = lambda y: y + node_id, way_id = node_start_id, way_start_id + osm_data: List[Dict] = [] + ANGLE_STEP = 1 + _edge_node_ids:Set[int] = set() + for i, now_row in df.iterrows(): + if now_row["名称"] == "直线": + _edge_node_ids.add(node_id) + n = { + "type": "node", + "id": node_id, + "x": float(now_row["起点 X"]), + "y": float(now_row["起点 Y"]), + } + osm_data.append(n) + node_id += 1 + + _edge_node_ids.add(node_id) + n = { + "type": "node", + "id": node_id, + "x": float(now_row["端点 X"]), + "y": float(now_row["端点 Y"]), + } + osm_data.append(n) + node_id += 1 + + w = { + "type": "way", + "id": way_id, + "nodes": [node_id - 2, node_id - 1], + "tags": {"highway": "tertiary"}, # here we set the same tag for all ways + "others": now_row, + } + way_id += 1 + + osm_data.append(w) + elif now_row["名称"] == "圆弧": + c_x = float(now_row["中心 X"]) + c_y = float(now_row["中心 Y"]) + start_angle = float(now_row["起点角度"]) + r = float(now_row["半径"]) + agl = 0 + nodes_indexes = [] + total_angle = float(now_row["总角度"]) + for agl in np.linspace(0, total_angle, max(int(total_angle / ANGLE_STEP), 3)): + n = { + "type": "node", + "id": node_id, + "x": c_x + r * math.cos(math.pi * (agl + start_angle) / 180), + "y": c_y + r * math.sin(math.pi * (agl + start_angle) / 180), + } + osm_data.append(n) + nodes_indexes.append(n["id"]) + node_id += 1 + _edge_node_ids.add(nodes_indexes[0]) + _edge_node_ids.add(nodes_indexes[-1]) + w = { + "type": "way", + "id": way_id, + "nodes": nodes_indexes, + "tags": {"highway": "tertiary"}, # here we set the same tag for all ways + "others": now_row, + } + osm_data.append(w) + way_id += 1 + for i in osm_data: + for _key, _transform_func in zip( + ["x", "y"], [x_transform_func, y_transform_func] + ): + if _key in i: + i[_key] = _transform_func(i[_key]) + to_merge_xys_dict = {i["id"]:(i["x"], i["y"]) for i in osm_data if "x" in i and "y" in i} + tree_id_to_node_id = {idx:i for idx,i in enumerate(to_merge_xys_dict.keys())} + tree = KDTree([v for v in to_merge_xys_dict.values()]) # type:ignore + merged = set() + father_id_dict = { + i: i for i in range(node_id) + } # Initialize the parent node id of each node + for _node_idx, xy in to_merge_xys_dict.items(): + if _node_idx in merged: + continue + a = [tree_id_to_node_id[i] for i in tree.query_ball_point(xy, merge_gate)] # type:ignore + if len(a) == 1: + unique_node_idx = a.pop() + father_id_dict[unique_node_idx] = _node_idx + merged.add(unique_node_idx) + else: + visited_nids = {_node_idx} + while len(a) > 0: + b = [] + for i in a: + if i in visited_nids: + continue + _xy = to_merge_xys_dict[i] + b.extend( + [tree_id_to_node_id[j] for j in tree.query_ball_point(_xy, merge_gate)] # type:ignore + ) + visited_nids.add(i) + father_id_dict[i] = _node_idx + a, b = b, [] + merged |= visited_nids + for i in range(node_id): + while father_id_dict[i] != father_id_dict[father_id_dict[i]]: + father_id_dict[i] = father_id_dict[father_id_dict[i]] + to_delete_ids = set() + for i in osm_data: + if i["type"] == "way": + i["nodes"] = [father_id_dict[n] for n in i["nodes"]] + if not len(set(i["nodes"]))>=2: + to_delete_ids.add(("way",i["id"])) + elif i["type"] == "node" and father_id_dict[i["id"]] != i["id"]: + to_delete_ids.add(("node",i["id"])) + return [i for i in osm_data if (i["type"],i["id"]) not in to_delete_ids] + + +logging.basicConfig( + level=logging.INFO, + format="%(asctime)s %(levelname)s %(message)s", +) +PROJ_STR = "+proj=tmerc +lat_0=33.9 +lon_0=116.4" +os.makedirs("cache", exist_ok=True) +rn = RoadNet( + proj_str=PROJ_STR, +) +CAD_PATH = "./data/cad/cad.xlsx" +osm_data = cad2osm( + CAD_PATH, + x_transform_func=lambda x: x - 229036.4782002, + y_transform_func=lambda y: y - 214014.32078879, +) +import pickle + +pickle.dump(osm_data, open("./cache/osm_data.pkl", "wb")) +print(Counter(i["type"] for i in osm_data)) +path = "cache/topo_from_cad.geojson" +net = rn.create_road_net(path, osm_data_cache=osm_data) + +builder = Builder( + net=net, + gen_sidewalk_speed_limit=50 / 3.6, + road_expand_mode="M", + proj_str=PROJ_STR, +) +m = builder.build("test") +pb = dict2pb(m, Map()) +with open("data/temp/cad_map.pb", "wb") as f: + f.write(pb.SerializeToString()) diff --git a/mosstool/map/_map_util/format_checker.py b/mosstool/map/_map_util/format_checker.py index 5439792..65e16d3 100644 --- a/mosstool/map/_map_util/format_checker.py +++ b/mosstool/map/_map_util/format_checker.py @@ -1,15 +1,15 @@ """Check the field type and value of map format""" import logging -import sys from collections import defaultdict +from typing import Dict, List, Optional from geojson import FeatureCollection from ...type import LaneType from .._map_util.const import * -__all__ = ["geojson_format_check", "output_format_check"] +__all__ = ["geojson_format_check", "output_format_check", "osm_format_checker"] class _FormatCheckHandler(logging.Handler): @@ -302,7 +302,7 @@ def geojson_format_check( handler.trigger_warnings() -def output_format_check(output_map: dict, output_lane_length_check: bool): +def output_format_check(output_map: dict, output_lane_length_check: bool) -> None: logging.basicConfig(level=logging.INFO) handler = _FormatCheckHandler() logger = logging.getLogger() @@ -411,3 +411,37 @@ def output_format_check(output_map: dict, output_lane_length_check: bool): f"Junction {lane_type} lane {lane_id} is too long ({length} m), please check input GeoJSON file!" ) handler.trigger_warnings() + + +def osm_format_checker( + osm_cache_check: bool, + osm_data: Optional[List[Dict]] = None, + required_keys_dict: Optional[Dict[str, List[str]]] = None, +) -> None: + if not osm_cache_check: + return + if osm_data is None: + return + logging.basicConfig(level=logging.INFO) + handler = _FormatCheckHandler() + logger = logging.getLogger() + logger.addHandler(handler) + if not hasattr(osm_data, "__iter__"): + logging.warning(f"input OSM data is not iterable!") + else: + if not all(isinstance(d, dict) for d in osm_data): + logging.warning(f"Not all items in OSM data is a `Dict`!") + else: + for _key in ["type", "id"]: + if not all(_key in d for d in osm_data): + logging.warning(f"insufficient key {_key} provided in OSM data!") + if required_keys_dict is None: + required_keys_dict = {} + for d in osm_data: + _keys = required_keys_dict.get(d["type"], []) + if any(_key not in d for _key in _keys): + logging.warning( + f"insufficient keys {_keys} provided in OSM data {d}!" + ) + + handler.trigger_warnings() diff --git a/mosstool/map/builder/builder.py b/mosstool/map/builder/builder.py index f3f0ab9..0695d9e 100644 --- a/mosstool/map/builder/builder.py +++ b/mosstool/map/builder/builder.py @@ -185,11 +185,16 @@ def __init__( ) # (N, 2) xyz_coords = np.column_stack([xy_coords, z_coords]) # (N, 3) else: + z_coords = ( + coords[:, 2] + if coords.shape[1] > 2 + else np.zeros((coords.shape[0], 1), dtype=np.float64) + ) xy_coords = np.array( [c[:2] for c in feature["geometry"]["coordinates"]], dtype=np.float64, ) - xyz_coords = coords + xyz_coords = np.column_stack([xy_coords, z_coords]) # (N, 3) feature["geometry"]["coordinates_xy"] = xy_coords feature["geometry"]["coordinates_xyz"] = xyz_coords if feature["geometry"]["type"] == "LineString": diff --git a/mosstool/map/osm/building.py b/mosstool/map/osm/building.py index 0875026..ee1f131 100644 --- a/mosstool/map/osm/building.py +++ b/mosstool/map/osm/building.py @@ -1,5 +1,5 @@ import logging -from typing import Dict, Optional +from typing import Dict, List, Optional, Tuple import pyproj import requests @@ -8,6 +8,8 @@ from shapely.geometry import MultiPolygon as sMultiPolygon from shapely.geometry import Polygon as sPolygon +from .._map_util.format_checker import osm_format_checker + __all__ = ["Building"] @@ -18,23 +20,23 @@ class Building: def __init__( self, - proj_str: str, - max_longitude: float, - min_longitude: float, - max_latitude: float, - min_latitude: float, + proj_str: Optional[str] = None, + max_longitude: Optional[float] = None, + min_longitude: Optional[float] = None, + max_latitude: Optional[float] = None, + min_latitude: Optional[float] = None, wikipedia_name: Optional[str] = None, proxies: Optional[Dict[str, str]] = None, ): """ Args: - - proj_str (str): projection string, e.g. 'epsg:3857' - - max_longitude (float): max longitude - - min_longitude (float): min longitude - - max_latitude (float): max latitude - - min_latitude (float): min latitude - - wikipedia_name (str): wikipedia name of the area in OSM. - - proxies (dict): proxies for requests, e.g. {'http': 'http://localhost:1080', 'https': 'http://localhost:1080'} + - proj_str (Optional[str]): projection string, e.g. 'epsg:3857' + - max_longitude (Optional[float]): max longitude + - min_longitude (Optional[float]): min longitude + - max_latitude (Optional[float]): max latitude + - min_latitude (Optional[float]): min latitude + - wikipedia_name (Optional[str]): wikipedia name of the area in OSM. + - proxies (Optional[Dict[str, str]]): proxies for requests, e.g. {'http': 'http://localhost:1080', 'https': 'http://localhost:1080'} """ self.bbox = ( min_latitude, @@ -42,7 +44,12 @@ def __init__( max_latitude, max_longitude, ) - self.projector = pyproj.Proj(proj_str) + self.proj_parameter = proj_str + self.projector = ( + pyproj.Proj(self.proj_parameter) if self.proj_parameter else None + ) + if self.projector is None and any(i is None for i in self.bbox): + logging.warning("Using osm cache data for input and xy coords for output!") self.wikipedia_name = wikipedia_name self.proxies = proxies # OSM raw data @@ -55,62 +62,76 @@ def __init__( # generate AOIs self.aois: list = [] - def _query_raw_data(self): + def _query_raw_data(self, osm_data_cache: Optional[List[Dict]] = None): """ Get raw data from OSM API OSM query language: https://wiki.openstreetmap.org/wiki/Overpass_API/Language_Guide Can be run and visualized in real time at https://overpass-turbo.eu/ """ - logging.info("Querying osm raw data") - bbox_str = ",".join(str(i) for i in self.bbox) - query_header = f"[out:json][timeout:120][bbox:{bbox_str}];" - area_wikipedia_name = self.wikipedia_name - if area_wikipedia_name is not None: - query_header += f'area[wikipedia="{area_wikipedia_name}"]->.searchArea;' - osm_data = None - for _ in range(3): # retry 3 times - try: - query_body_raw = [ - ( - "way", - '[!highway][!tunnel][!boundary][!railway][!natural][!barrier][!junction][!waterway][!public_transport][landuse!~"grass|forest"][place!~"suburb|neighbourhood"][amenity!=fountain][historic!=city_gate][artwork_type!=sculpture][man_made!~"bridge|water_well"][building!~"wall|train_station|roof"]', - ), - ( - "rel", - '[landuse][landuse!~"grass|forest"][type=multipolygon]', - ), - ( - "rel", - "[amenity][amenity!=fountain][type=multipolygon]", - ), - ( - "rel", - '[building][building!~"wall|train_station|roof"][type=multipolygon]', - ), - ( - "rel", - "[leisure][leisure!=nature_reserve][type=multipolygon]", - ), - ] - query_body = "" - for obj, args in query_body_raw: - area = ( - "(area.searchArea)" if area_wikipedia_name is not None else "" - ) - query_body += obj + area + args + ";" - query_body = "(" + query_body + ");" - query = query_header + query_body + "(._;>;);" + "out body;" - logging.info(f"{query}") - osm_data = requests.get( - "http://overpass-api.de/api/interpreter?data=" + query, - proxies=self.proxies, - ).json()["elements"] - break - except Exception as e: - logging.warning(f"Exception when querying OSM data {e}") - logging.warning("No response from OSM, Please try again later!") - if osm_data is None: - raise Exception("No AOI response from OSM!") + if osm_data_cache is None: + logging.info("Querying osm raw data") + assert all( + i is not None for i in self.bbox + ), f"longitude and latitude are required without cache file!" + (min_lat, min_lon, max_lat, max_lon) = self.bbox + if self.proj_parameter is None and self.projector is None: + proj_str = f"+proj=tmerc +lat_0={(max_lat+min_lat)/2} +lon_0={(max_lon+min_lon)/2}" # type: ignore + self.proj_parameter = proj_str + self.projector = pyproj.Proj(self.proj_parameter) + bbox_str = ",".join(str(i) for i in self.bbox) + query_header = f"[out:json][timeout:120][bbox:{bbox_str}];" + area_wikipedia_name = self.wikipedia_name + if area_wikipedia_name is not None: + query_header += f'area[wikipedia="{area_wikipedia_name}"]->.searchArea;' + osm_data = None + for _ in range(3): # retry 3 times + try: + query_body_raw = [ + ( + "way", + '[!highway][!tunnel][!boundary][!railway][!natural][!barrier][!junction][!waterway][!public_transport][landuse!~"grass|forest"][place!~"suburb|neighbourhood"][amenity!=fountain][historic!=city_gate][artwork_type!=sculpture][man_made!~"bridge|water_well"][building!~"wall|train_station|roof"]', + ), + ( + "rel", + '[landuse][landuse!~"grass|forest"][type=multipolygon]', + ), + ( + "rel", + "[amenity][amenity!=fountain][type=multipolygon]", + ), + ( + "rel", + '[building][building!~"wall|train_station|roof"][type=multipolygon]', + ), + ( + "rel", + "[leisure][leisure!=nature_reserve][type=multipolygon]", + ), + ] + query_body = "" + for obj, args in query_body_raw: + area = ( + "(area.searchArea)" + if area_wikipedia_name is not None + else "" + ) + query_body += obj + area + args + ";" + query_body = "(" + query_body + ");" + query = query_header + query_body + "(._;>;);" + "out body;" + logging.info(f"{query}") + osm_data = requests.get( + "http://overpass-api.de/api/interpreter?data=" + query, + proxies=self.proxies, + ).json()["elements"] + break + except Exception as e: + logging.warning(f"Exception when querying OSM data {e}") + logging.warning("No response from OSM, Please try again later!") + if osm_data is None: + raise Exception("No AOI response from OSM!") + else: + logging.info("Loading osm raw data from cache") + osm_data = osm_data_cache nodes = [d for d in osm_data if d["type"] == "node"] ways = [d for d in osm_data if d["type"] == "way"] rels = [ @@ -145,10 +166,19 @@ def _make_raw_aoi(self): ways_aoi = self._ways_aoi ways_rel = self._ways_rel logging.info("Making raw aoi") - nodes_dict = { - n["id"]: [a for a in self.projector(n["lon"], n["lat"], inverse=False)] - for n in nodes - } + nodes_dict = {} + for n in nodes: + if "x" in n and "y" in n: + nodes_dict[n["id"]] = [a for a in (n["x"], n["y"])] + elif "lon" in n and "lat" in n: + assert ( + self.projector is not None + ), f"proj_str are required when downloading from OSM!" + nodes_dict[n["id"]] = [ + a for a in self.projector(n["lon"], n["lat"], inverse=False) + ] + else: + raise ValueError(f"Neither lon and lat or x and y in node {n}!") ways_aoi_dict = {w["id"]: w for w in ways_aoi} # nodes, tags ways_rel_dict = {w["id"]: w for w in ways_rel} rels_dict = {r["id"]: r for r in rels} # members: [(ref, role)], tags @@ -311,17 +341,31 @@ def _make_raw_aoi(self): logging.info(f"invalid_relation_cnt: {invalid_relation_cnt}") logging.info(f"raw aoi: {len(self.aois)}") - def create_building(self, output_path: Optional[str] = None): + def _transform_coordinate(self, c: Tuple[float, float]) -> List[float]: + if self.projector is None: + return [c[0], c[1]] + else: + return list(self.projector(c[0], c[1], inverse=True)) + + def create_building( + self, + output_path: Optional[str] = None, + osm_data_cache: Optional[List[Dict]] = None, + osm_cache_check: bool = False, + ): """ Create AOIs from OpenStreetMap. Args: + - osm_data_cache (Optional[List[Dict]]): OSM data cache. - output_path (str): GeoJSON file output path. + - osm_cache_check (bool): check the format of input OSM data cache. Returns: - AOIs in GeoJSON format. """ - self._query_raw_data() + osm_format_checker(osm_cache_check, osm_data_cache) + self._query_raw_data(osm_data_cache) self._make_raw_aoi() geos = [] for aoi in self.aois: @@ -330,7 +374,7 @@ def create_building(self, output_path: Optional[str] = None): geometry=Polygon( [ [ - list(self.projector(c[0], c[1], inverse=True)) + self._transform_coordinate(c) for c in aoi["geo"].exterior.coords ] ] diff --git a/mosstool/map/osm/point_of_interest.py b/mosstool/map/osm/point_of_interest.py index c7aa108..221e15e 100644 --- a/mosstool/map/osm/point_of_interest.py +++ b/mosstool/map/osm/point_of_interest.py @@ -1,9 +1,11 @@ import logging -from typing import Dict, Optional +from typing import Dict, List, Optional import requests from geojson import Feature, FeatureCollection, Point, dump +from .._map_util.format_checker import osm_format_checker + __all__ = ["PointOfInterest"] @@ -14,21 +16,21 @@ class PointOfInterest: def __init__( self, - max_longitude: float, - min_longitude: float, - max_latitude: float, - min_latitude: float, + max_longitude: Optional[float] = None, + min_longitude: Optional[float] = None, + max_latitude: Optional[float] = None, + min_latitude: Optional[float] = None, wikipedia_name: Optional[str] = None, proxies: Optional[Dict[str, str]] = None, ): """ Args: - - max_longitude (float): max longitude - - min_longitude (float): min longitude - - max_latitude (float): max latitude - - min_latitude (float): min latitude - - wikipedia_name (str): wikipedia name of the area in OSM. - - proxies (dict): proxies for requests, e.g. {'http': 'http://localhost:1080', 'https': 'http://localhost:1080'} + - max_longitude (Optional[float]): max longitude + - min_longitude (Optional[float]): min longitude + - max_latitude (Optional[float]): max latitude + - min_latitude (Optional[float]): min latitude + - wikipedia_name (Optional[str]): wikipedia name of the area in OSM. + - proxies (Optional[Dict[str, str]]): proxies for requests, e.g. {'http': 'http://localhost:1080', 'https': 'http://localhost:1080'} """ self.bbox = ( min_latitude, @@ -44,57 +46,66 @@ def __init__( # generate POIs self.pois: list = [] - def _query_raw_data(self): + def _query_raw_data(self, osm_data_cache: Optional[List[Dict]] = None): """ Get raw data from OSM API OSM query language: https://wiki.openstreetmap.org/wiki/Overpass_API/Language_Guide Can be run and visualized in real time at https://overpass-turbo.eu/ """ - logging.info("Querying osm raw data") - bbox_str = ",".join(str(i) for i in self.bbox) - query_header = f"[out:json][timeout:120][bbox:{bbox_str}];" - area_wikipedia_name = self.wikipedia_name - if area_wikipedia_name is not None: - query_header += f'area[wikipedia="{area_wikipedia_name}"]->.searchArea;' - osm_data = None - for _ in range(3): # retry 3 times - try: - query_body_raw = [ - ( - "node", - "", - ), - ] - query_body = "" - for obj, args in query_body_raw: - area = ( - "(area.searchArea)" if area_wikipedia_name is not None else "" - ) - query_body += obj + area + args + ";" - query_body = "(" + query_body + ");" - query = query_header + query_body + "(._;>;);" + "out body;" - logging.info(f"{query}") - osm_data = requests.get( - "http://overpass-api.de/api/interpreter?data=" + query, - proxies=self.proxies, - ).json()["elements"] - break - except Exception as e: - logging.warning(f"Exception when querying OSM data {e}") - logging.warning("No response from OSM, Please try again later!") - if osm_data is None: - raise Exception("No POI response from OSM!") + if osm_data_cache is None: + logging.info("Querying osm raw data") + assert all( + i is not None for i in self.bbox + ), f"longitude and latitude are required without cache file!" + bbox_str = ",".join(str(i) for i in self.bbox) + query_header = f"[out:json][timeout:120][bbox:{bbox_str}];" + area_wikipedia_name = self.wikipedia_name + if area_wikipedia_name is not None: + query_header += f'area[wikipedia="{area_wikipedia_name}"]->.searchArea;' + osm_data = None + for _ in range(3): # retry 3 times + try: + query_body_raw = [ + ( + "node", + "", + ), + ] + query_body = "" + for obj, args in query_body_raw: + area = ( + "(area.searchArea)" + if area_wikipedia_name is not None + else "" + ) + query_body += obj + area + args + ";" + query_body = "(" + query_body + ");" + query = query_header + query_body + "(._;>;);" + "out body;" + logging.info(f"{query}") + osm_data = requests.get( + "http://overpass-api.de/api/interpreter?data=" + query, + proxies=self.proxies, + ).json()["elements"] + break + except Exception as e: + logging.warning(f"Exception when querying OSM data {e}") + logging.warning("No response from OSM, Please try again later!") + if osm_data is None: + raise Exception("No POI response from OSM!") + else: + osm_data = osm_data_cache nodes = [d for d in osm_data if d["type"] == "node"] logging.info(f"node: {len(nodes)}") self._nodes = nodes + def _make_raw_poi(self): """ Construct POI from original OSM data. """ - _raw_pois = [] + _raw_pois = [] for d in self._nodes: - d_tags = d.get("tags",{}) - p_name = d_tags.get("name","") + d_tags = d.get("tags", {}) + p_name = d_tags.get("name", "") # name d["name"] = p_name # catg @@ -129,30 +140,37 @@ def _make_raw_poi(self): continue logging.info(f"raw poi: {len(_raw_pois)}") self.pois = _raw_pois - def create_pois(self, output_path: Optional[str] = None): + + def create_pois( + self, + output_path: Optional[str] = None, + osm_data_cache: Optional[List[Dict]] = None, + osm_cache_check: bool = False, + ): """ Create POIs from OpenStreetMap. Args: + - osm_data_cache (Optional[List[Dict]]): OSM data cache. - output_path (str): GeoJSON file output path. + - osm_cache_check (bool): check the format of input OSM data cache. Returns: - POIs in GeoJSON format. """ - self._query_raw_data() + osm_format_checker(osm_cache_check, osm_data_cache, {"node": ["lon", "lat"]}) + self._query_raw_data(osm_data_cache) self._make_raw_poi() geos = [] - for poi_id,poi in enumerate(self.pois): + for poi_id, poi in enumerate(self.pois): geos.append( Feature( - geometry=Point( - [poi["lon"],poi["lat"]] - ), + geometry=Point([poi["lon"], poi["lat"]]), properties={ "id": poi_id, "osm_tags": poi["tags"], - "name":poi["name"], - "catg":poi["catg"], + "name": poi["name"], + "catg": poi["catg"], }, ) ) diff --git a/mosstool/map/osm/roadnet.py b/mosstool/map/osm/roadnet.py index 5d13eb9..7929b5e 100644 --- a/mosstool/map/osm/roadnet.py +++ b/mosstool/map/osm/roadnet.py @@ -10,8 +10,8 @@ import requests from geojson import Feature, FeatureCollection, LineString, MultiPoint, dump from shapely.geometry import LineString as sLineString -from tqdm import tqdm +from .._map_util.format_checker import osm_format_checker from .._map_util.osm_const import * from ._motif import close_nodes, motif_H, suc_is_close_by_other_way from ._wayutil import merge_way_nodes, parse_osm_way_tags @@ -26,33 +26,37 @@ class RoadNet: def __init__( self, - proj_str: str, - max_longitude: float, - min_longitude: float, - max_latitude: float, - min_latitude: float, + proj_str: Optional[str] = None, + max_longitude: Optional[float] = None, + min_longitude: Optional[float] = None, + max_latitude: Optional[float] = None, + min_latitude: Optional[float] = None, wikipedia_name: Optional[str] = None, proxies: Optional[Dict[str, str]] = None, ): """ Args: - - proj_str (str): projection string, e.g. 'epsg:3857' - - max_longitude (float): max longitude - - min_longitude (float): min longitude - - max_latitude (float): max latitude - - min_latitude (float): min latitude - - wikipedia_name (str): wikipedia name of the area in OSM. - - proxies (dict): proxies for requests, e.g. {'http': 'http://localhost:1080', 'https': 'http://localhost:1080'} + - proj_str (Optional[str]): projection string, e.g. 'epsg:3857' + - max_longitude (Optional[float]): max longitude + - min_longitude (Optional[float]): min longitude + - max_latitude (Optional[float]): max latitude + - min_latitude (Optional[float]): min latitude + - wikipedia_name (Optional[str]): wikipedia name of the area in OSM. + - proxies (Optional[Dict[str, str]]): proxies for requests, e.g. {'http': 'http://localhost:1080', 'https': 'http://localhost:1080'} """ # configs self.proj_parameter = proj_str - self.projector = pyproj.Proj(self.proj_parameter) + self.projector = ( + pyproj.Proj(self.proj_parameter) if proj_str is not None else None + ) self.bbox = ( min_latitude, min_longitude, max_latitude, max_longitude, ) + if self.projector is None and any(i is None for i in self.bbox): + logging.warning("Using osm cache data for input and xy coords for output!") self.proxies = proxies self.wikipedia_name = wikipedia_name self.way_filter = WAY_FILTER @@ -76,6 +80,14 @@ def default_way_settings(self): def _download_osm(self): """Fetch raw data from OpenStreetMap""" + assert all( + i is not None for i in self.bbox + ), f"longitude and latitude are required without cache file!" + (min_lat, min_lon, max_lat, max_lon) = self.bbox + if self.proj_parameter is None and self.projector is None: + proj_str = f"+proj=tmerc +lat_0={(max_lat+min_lat)/2} +lon_0={(max_lon+min_lon)/2}" # type: ignore + self.proj_parameter = proj_str + self.projector = pyproj.Proj(self.proj_parameter) bbox_str = ",".join(str(i) for i in self.bbox) query = f"[out:json][timeout:180][bbox:{bbox_str}];" wikipedia_name = self.wikipedia_name @@ -252,21 +264,36 @@ def to_topo(self): return FeatureCollection(geos) - def _get_osm(self): - osm_data = self._download_osm() + def _get_osm(self, osm_data_cache: Optional[List[Dict]] = None): + if osm_data_cache is not None: + osm_data = osm_data_cache + else: + osm_data = self._download_osm() # find every valuable nodes self.nodes = {} - for doc in tqdm(osm_data): + for doc in osm_data: if doc["type"] != "node": continue - if "lon" in doc and "lat" in doc: + if "x" in doc and "y" in doc: + # doc["x"], doc["y"] = round(doc["x"], 6), round(doc["y"], 6) + if self.projector is not None: + doc["lon"], doc["lat"] = self.projector( + doc["x"], doc["y"], inverse=True + ) + else: + doc["lon"], doc["lat"] = doc["x"], doc["y"] + self.nodes[doc["id"]] = doc + elif "lon" in doc and "lat" in doc: # Latitude and longitude are specified with a precision of six decimal places. doc["lon"], doc["lat"] = round(doc["lon"], 6), round(doc["lat"], 6) + assert ( + self.projector is not None + ), f"proj_str are required when downloading from OSM!" doc["x"], doc["y"] = self.projector(doc["lon"], doc["lat"]) self.nodes[doc["id"]] = doc # all ways self.ways = {} - for way in tqdm(osm_data): + for way in osm_data: # Filter non-way elements if way["type"] != "way": continue @@ -297,7 +324,7 @@ def _get_osm(self): way_id = max(self.ways) # New road collection new_ways = {} - for way in tqdm(self.ways.values()): + for way in self.ways.values(): last_pos = 0 for pos, node in list(enumerate(way["nodes"]))[1:-1]: if node in joints: @@ -442,45 +469,48 @@ def _remove_simple_joints(self): # Two-way processing, find all connected components, confirm that there are 2 points with degree 1, and the rest are points with degree 2 # Select any point with degree 1, calculate the path to another point, and then assemble the new way # Find all Unicom components - for component in list(nx.connected_components(g)): - # Use assert to confirm that there are 2 nodes with degree 1, and the rest are nodes with degree 2 - degree_one_count = 0 - degree_two_count = 0 - for node in component: - if g.degree(node) == 1: # type: ignore - degree_one_count += 1 - elif g.degree(node) == 2: # type: ignore - degree_two_count += 1 - else: - raise AssertionError( - "The degree of the node does not meet the requirements" - ) - assert ( - degree_one_count == 2 - ), "The number of nodes with degree 1 is incorrect" - assert ( - degree_two_count == len(component) - 2 - ), "The number of nodes with degree 2 is incorrect" - # Select a node with degree 1 and calculate the path to another node with degree 1 - degree_one_nodes = [ - node for node in component if g.degree(node) == 1 # type: ignore - ] - start_node = degree_one_nodes[0] - end_node = degree_one_nodes[1] - way_ids = nx.shortest_path(g, start_node, end_node) - # Assemble new ways and delete redundant ways - main_way_id = way_ids[0] - main_way = self.ways[main_way_id] - main_way["remove_simple_joints"] = "main" - for way_id in way_ids[1:]: - way = self.ways[way_id] - try: - main_way["nodes"] = merge_way_nodes(main_way["nodes"], way["nodes"]) - except ValueError as e: - # with open("cache/graph.pkl", "wb") as f: - # pickle.dump([g, component], f) - raise e - removed_ids.add(way_id) + try: + for component in list(nx.connected_components(g)): + # Use assert to confirm that there are 2 nodes with degree 1, and the rest are nodes with degree 2 + degree_one_count = 0 + degree_two_count = 0 + for node in component: + if g.degree(node) == 1: # type: ignore + degree_one_count += 1 + elif g.degree(node) == 2: # type: ignore + degree_two_count += 1 + else: + raise AssertionError( + "The degree of the node does not meet the requirements" + ) + assert ( + degree_one_count == 2 + ), "The number of nodes with degree 1 is incorrect" + assert ( + degree_two_count == len(component) - 2 + ), "The number of nodes with degree 2 is incorrect" + # Select a node with degree 1 and calculate the path to another node with degree 1 + degree_one_nodes = [ + node for node in component if g.degree(node) == 1 # type: ignore + ] + start_node = degree_one_nodes[0] + end_node = degree_one_nodes[1] + way_ids = nx.shortest_path(g, start_node, end_node) + # Assemble new ways and delete redundant ways + main_way_id = way_ids[0] + main_way = self.ways[main_way_id] + main_way["remove_simple_joints"] = "main" + for way_id in way_ids[1:]: + way = self.ways[way_id] + try: + main_way["nodes"] = merge_way_nodes(main_way["nodes"], way["nodes"]) + except ValueError as e: + # with open("cache/graph.pkl", "wb") as f: + # pickle.dump([g, component], f) + raise e + removed_ids.add(way_id) + except Exception as e: + logging.warning(f"Error when trying to simplify roadnet {e}!") for wid in removed_ids: del self.ways[wid] self._assert_ways() @@ -668,18 +698,26 @@ def _clean_topo(self): del self.node2junction[node] del self.junction2nodes[c] - def create_road_net(self, output_path: Optional[str] = None): + def create_road_net( + self, + output_path: Optional[str] = None, + osm_data_cache: Optional[List[Dict]] = None, + osm_cache_check: bool = False, + ): """ Create Road net from OpenStreetMap. Args: - - output_path (str): GeoJSON file output path. + - osm_data_cache (Optional[List[Dict]]): OSM data cache. + - output_path (Optional[str]): GeoJSON file output path. + - osm_cache_check (bool): check the format of input OSM data cache. Returns: - roads and junctions in GeoJSON format. """ + osm_format_checker(osm_cache_check, osm_data_cache) # Get OSM data - self._get_osm() + self._get_osm(osm_data_cache) # self.dump_as_geojson("cache/1.geojson") # Process diff --git a/pyproject.toml b/pyproject.toml index fb6b325..e09d32e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "mosstool" -version = "1.0.12" +version = "1.0.13" description = "MObility Simulation System toolbox " authors = ["Jun Zhang "] license = "MIT"