diff --git a/examples/map_osm2geojson.py b/examples/map_osm2geojson.py index 84886ee..c546fdb 100644 --- a/examples/map_osm2geojson.py +++ b/examples/map_osm2geojson.py @@ -13,8 +13,6 @@ } # Configure log and store it in the file mapbuilder2.log logging.basicConfig( - filename="mapbuilder2.log", - filemode="w", level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s", ) diff --git a/examples/merge_map.py b/examples/merge_map.py index 4885f72..46516ec 100644 --- a/examples/merge_map.py +++ b/examples/merge_map.py @@ -1,19 +1,36 @@ import logging from mosstool.type import Map +from mosstool.map.builder import Builder from mosstool.util.map_merger import merge_map +from mosstool.util.format_converter import dict2pb logging.basicConfig( level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s", ) + +REBUILD_MAP = True maps = [] for i in [0, 1]: with open(f"data/temp/test_{i}.pb", "rb") as f: pb = Map() pb.ParseFromString(f.read()) maps.append(pb) -merge_map( +merged_pb = merge_map( partial_maps=maps, output_path="data/temp/merged_m.pb", ) +if REBUILD_MAP: + builder = Builder( + net=merged_pb, + proj_str=merged_pb.header.projection, + gen_sidewalk_speed_limit=50 / 3.6, + aoi_mode="append", # keep AOIs in merged_pb + road_expand_mode="M", +) + rebuild_m = builder.build(merged_pb.header.name) + rebuild_pb = dict2pb(rebuild_m, Map()) + with open("data/temp/rebuild_merged_m.pb", "wb") as f: + f.write(rebuild_pb.SerializeToString()) + diff --git a/mosstool/map/_map_util/aoi_matcher.py b/mosstool/map/_map_util/aoi_matcher.py index d1808bd..73d6540 100644 --- a/mosstool/map/_map_util/aoi_matcher.py +++ b/mosstool/map/_map_util/aoi_matcher.py @@ -11,25 +11,14 @@ import shapely.ops as ops from scipy.spatial import KDTree from shapely.affinity import scale -from shapely.geometry import ( - LineString, - MultiLineString, - MultiPoint, - MultiPolygon, - Point, - Polygon, -) +from shapely.geometry import (LineString, MultiLineString, MultiPoint, + MultiPolygon, Point, Polygon) from shapely.strtree import STRtree from ...type import AoiType from .._util.angle import abs_delta_angle, delta_angle -from .._util.line import ( - connect_line_string, - get_line_angle, - get_start_vector, - line_extend, - offset_lane, -) +from .._util.line import (connect_line_string, get_line_angle, + get_start_vector, line_extend, offset_lane) from .aoiutils import geo_coords # ATTENTION: In order to achieve longer distance POI merging, the maximum recursion depth needs to be modified. @@ -742,12 +731,10 @@ def _add_poly_aoi_unit(arg): "positions": [{"x": c[0], "y": c[1]} for c in geo_coords(geo)], "area": geo.area, "external": { - "osm_tencent_ids": [ - aoi["id"] - ], # id in aoi_osm_tencent_fudan.aoi_beijing5ring + "osm_input_ids": [aoi["id"]], "ex_poi_ids": aoi["external"].get( "inner_poi", [] - ), # Tencent id of poi contained in aoi + ), # Input id of poi contained in aoi # aoi population "population": aoi["external"].get("population", 0), "land_types": aoi["external"].get("land_types", defaultdict(float)), @@ -849,14 +836,12 @@ def _add_aoi_stop_unit(arg): "positions": [{"x": c[0], "y": c[1]} for c in geo_coords(geo)], "area": geo.area, "external": { - "osm_tencent_ids": [ - aoi["id"] - ], # id in aoi_osm_tencent_fudan.aoi_beijing5ring + "osm_input_ids": [aoi["id"]], "stop_id": aoi["external"]["stop_id"], "station_type": station_type, "ex_poi_ids": aoi["external"].get( "inner_poi", [] - ), # Tencent id of poi contained in aoi + ), # Input id of poi contained in aoi # aoi population "population": aoi["external"].get("population", 0), "land_types": aoi["external"].get("land_types", defaultdict(float)), @@ -1630,7 +1615,7 @@ def _add_pois(aois, pois, projstr): "position": {"x": x, "y": y}, "aoi_id": poi_id_to_aoi_id[poi_id], "external": { - "tencent_poi_id": poi_id, + "input_poi_id": poi_id, }, } poi_uid += 1 @@ -1890,10 +1875,10 @@ def add_aoi_to_map( raw_pois = {doc["id"]: doc for doc in input_pois} aois = _add_aoi(aois=input_aois, stops=input_stops, workers=workers, merge_aoi=True) - added_tencent_poi = [] + added_input_poi = [] for _, aoi in aois.items(): - added_tencent_poi.extend(aoi["external"].get("ex_poi_ids", [])) - logging.info(f"Added Tencent POI num: {len(added_tencent_poi)}") + added_input_poi.extend(aoi["external"].get("ex_poi_ids", [])) + logging.info(f"Added Input POI num: {len(added_input_poi)}") # Post-compute logging.info("Post Compute") aois = _add_aoi_land_use(aois, shp_path, bbox, projstr) diff --git a/mosstool/map/builder/builder.py b/mosstool/map/builder/builder.py index 91a390e..ba54014 100644 --- a/mosstool/map/builder/builder.py +++ b/mosstool/map/builder/builder.py @@ -5,7 +5,8 @@ from copy import deepcopy from math import atan2 from multiprocessing import cpu_count -from typing import Callable, Dict, List, Literal, Optional, Set, Tuple, Union, cast +from typing import (Callable, Dict, List, Literal, Optional, Set, Tuple, Union, + cast) import matplotlib.pyplot as plt import networkx as nx @@ -25,21 +26,15 @@ from .._map_util.aoiutils import generate_aoi_poi from .._map_util.const import * from .._map_util.convert_aoi import convert_aoi, convert_poi -from .._map_util.format_checker import geojson_format_check, output_format_check +from .._map_util.format_checker import (geojson_format_check, + output_format_check) from .._map_util.gen_traffic_light import generate_traffic_light from .._map_util.junctionutils import add_driving_groups, add_overlaps from .._map_util.map_aois_matchers import match_map_aois from .._util.angle import abs_delta_angle, delta_angle -from .._util.line import ( - align_line, - connect_line_string, - get_line_angle, - get_start_vector, - line_extend, - line_max_curvature, - merge_line_start_end, - offset_lane, -) +from .._util.line import (align_line, connect_line_string, get_line_angle, + get_start_vector, line_extend, line_max_curvature, + merge_line_start_end, offset_lane) __all__ = ["Builder"] @@ -1543,8 +1538,8 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): ## junc_lanes = [] # Due to the difference between main roads and auxiliary roads, there are 4 possibilities for connected wids depending on the presence or absence of auxiliary roads. - # 1.Only main road - # 2.in only has the main road and out has the main road and the auxiliary road + # 1. Only main road + # 2. in only has the main road and out has the main road and the auxiliary road # 3. There are main roads and auxiliary roads in and only main roads out. # 4. There are main roads and auxiliary roads in and there are main roads and auxiliary roads out. for in_angle, in_way_ids in in_way_groups: @@ -1571,9 +1566,6 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): ] else: in_auxiliary_turn_config = [] - # if jid in bad_turn_config_jids: - # in_main_turn_config = [] - # in_auxiliary_turn_config = [] # Default turn annotation # Go straight default_in_straight_main_lanes = in_main_lanes @@ -1591,10 +1583,12 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): if len(in_main_lanes) > LARGE_LANE_NUM_THRESHOLD else DEFAULT_TURN_NUM["MAIN_SMALL_RIGHT"] ) + USE_TURN_CONFIG = False # Mark turn according to turn config if in_main_turn_config and len(in_main_turn_config) == len( in_main_lanes ): + USE_TURN_CONFIG = True # Go straight in_straight_main_lanes = [ l @@ -1614,18 +1608,9 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): main_right_lane_num = len( [t for t in in_main_turn_config if "R" in t or "r" in t] ) - - # - - # if ( - # ( - # len(out_around_groups) == 0 and ((main_around_lane_num > 0)) - # ) - # or (len(out_left_groups) == 0 and ((main_left_lane_num > 0))) - # or (len(out_right_groups) == 0 and ((main_right_lane_num > 0))) - # ): # main road if jid in bad_turn_config_jids: + USE_TURN_CONFIG = False # straight in_straight_main_lanes = ( default_in_straight_main_lanes @@ -1675,6 +1660,7 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): if in_auxiliary_turn_config and len( in_auxiliary_turn_config ) == len(in_auxiliary_lanes): + USE_TURN_CONFIG = True # Go straight in_straight_auxiliary_lanes = [ l @@ -1706,23 +1692,9 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): if "R" in t or "r" in t ] ) - # Determine whether turn_config can be used to connect to the lane. If the default configuration cannot be used, - # if ( - # ( - # len(out_around_groups) == 0 - # and ((auxiliary_around_lane_num > 0)) - # ) - # or ( - # len(out_left_groups) == 0 - # and ((auxiliary_left_lane_num > 0)) - # ) - # or ( - # len(out_right_groups) == 0 - # and ((auxiliary_right_lane_num > 0)) - # ) - # ): # auxiliary road if jid in bad_turn_config_jids: + USE_TURN_CONFIG = False # straight in_straight_auxiliary_lanes = ( default_in_straight_auxiliary_lanes @@ -1755,6 +1727,8 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): auxiliary_left_lane_num = default_auxiliary_left_lane_num auxiliary_right_lane_num = default_auxiliary_right_lane_num # All lanes connecting incoming and outgoing main roads + MAIN_HAS_STRAIGHT_LANES = False + AUXILIARY_HAS_STRAIGHT_LANES = False if out_main_group is not None: out_angle, out_way_ids = out_main_group out_main_wid, out_auxiliary_wids = classify_main_auxiliary_wid( @@ -1768,6 +1742,7 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): ] # main road to main road if len(in_straight_main_lanes) > 0 and len(out_main_lanes) > 0: + MAIN_HAS_STRAIGHT_LANES = True junc_lanes += self._connect_lane_group( in_lanes=in_straight_main_lanes, out_lanes=out_main_lanes, @@ -1780,6 +1755,7 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): len(in_straight_auxiliary_lanes) > 0 and len(out_auxiliary_lanes) > 0 ): + AUXILIARY_HAS_STRAIGHT_LANES = True junc_lanes += self._connect_lane_group( in_lanes=in_straight_auxiliary_lanes, out_lanes=out_auxiliary_lanes, @@ -1908,6 +1884,18 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): # Connect left ramp if len(out_left_groups) > 0: + # ATTENTION: add # of left lanes when no turn configs provided + if not USE_TURN_CONFIG: + if not MAIN_HAS_STRAIGHT_LANES: + main_left_lane_num = max( + main_left_lane_num, + len(in_main_lanes) - main_right_lane_num, + ) + if not AUXILIARY_HAS_STRAIGHT_LANES: + auxiliary_left_lane_num = max( + auxiliary_left_lane_num, + len(in_auxiliary_lanes) - auxiliary_right_lane_num, + ) for out_angle, out_way_ids in out_left_groups: ( out_main_wid, @@ -2303,6 +2291,7 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): for in_way_id in in_way_ids: in_way_id2available_conn[in_way_id] = available_groups ## + USE_TURN_CONFIG = False # No side roads if not len(in_auxiliary_lanes) > 0: # total number @@ -2350,6 +2339,7 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): if in_main_turn_config and len(in_main_turn_config) == len( in_main_lanes ): + USE_TURN_CONFIG = True # Go straight main_indexes = [ i @@ -2377,12 +2367,8 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): right_count = len( [t for t in in_main_turn_config if "R" in t or "r" in t] ) - # Determine whether turn_config can be used to connect to lane. If the default configuration cannot be used, - # if ((len(out_around_groups) == 0 and (around_count > 0)) - # or (len(out_left_groups) == 0 and (left_count > 0)) - # or (len(out_right_groups) == 0 and (right_count > 0)) - # ): if jid in bad_turn_config_jids: + USE_TURN_CONFIG = False around_count = ( default_around_count if around_count == 0 @@ -2411,6 +2397,17 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): main_out_start = default_main_out_start main_out_end = default_main_out_end # connect lanes + MAIN_HAS_STRAIGHT_LANES = False + if main_count > 0 and out_main_group: + MAIN_HAS_STRAIGHT_LANES = True + out_road = self.map_roads[out_main_group[1][0]] + junc_lanes += self._connect_lane_group( + in_lanes=in_main_lanes[main_out_start:main_out_end], + out_lanes=out_road["lanes"], + lane_turn=mapv2.LANE_TURN_STRAIGHT, + lane_type=mapv2.LANE_TYPE_DRIVING, + junc_id=j_uid, + ) if around_count > 0: for out_angle, out_way_ids in out_around_groups: for wid in out_way_ids: @@ -2433,6 +2430,10 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): lane_type=mapv2.LANE_TYPE_DRIVING, junc_id=j_uid, ) + if not USE_TURN_CONFIG and not MAIN_HAS_STRAIGHT_LANES: + left_count = max( + left_count, len(in_main_lanes) - right_count + ) if left_count > 0: for out_angle, out_way_ids in out_left_groups: for wid in out_way_ids: @@ -2444,15 +2445,6 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): lane_type=mapv2.LANE_TYPE_DRIVING, junc_id=j_uid, ) - if main_count > 0 and out_main_group: - out_road = self.map_roads[out_main_group[1][0]] - junc_lanes += self._connect_lane_group( - in_lanes=in_main_lanes[main_out_start:main_out_end], - out_lanes=out_road["lanes"], - lane_turn=mapv2.LANE_TURN_STRAIGHT, - lane_type=mapv2.LANE_TYPE_DRIVING, - junc_id=j_uid, - ) else: # There is a auxiliary road # main road # Default turn annotation @@ -2495,10 +2487,12 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): default_main_out_start, default_main_out_end = 0, len( in_main_lanes ) + USE_TURN_CONFIG = False # main road if in_main_turn_config and len(in_main_turn_config) == len( in_main_lanes ): + USE_TURN_CONFIG = True # Go straight main_indexes = [ i @@ -2527,14 +2521,8 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): [t for t in in_main_turn_config if "R" in t or "r" in t] ) - # # Determine whether turn_config can be used to connect to lane. If the default configuration cannot be used, - # if ( - # (len(out_around_groups) == 0 and (around_count > 0)) - # or (len(out_left_groups) == 0 and left_count > 0) - # or (len(out_right_groups) == 0 and right_count > 0) - # ): - if jid in bad_turn_config_jids: + USE_TURN_CONFIG = False around_count = ( default_around_count if around_count == 0 @@ -2626,6 +2614,7 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): if in_auxiliary_turn_config and len( in_auxiliary_turn_config ) == len(in_auxiliary_lanes): + USE_TURN_CONFIG = True # Go straight auxiliary_main_indexes = [ i @@ -2667,23 +2656,8 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): if "R" in t or "r" in t ] ) - # Determine whether turn_config can be used to connect to lane. If the default configuration cannot be used, - # if ( - # ( - # len(out_around_groups) == 0 - # and ((auxiliary_around_count > 0)) - # ) - # or ( - # len(out_left_groups) == 0 - # and ((auxiliary_left_count > 0)) - # ) - # or ( - # len(out_right_groups) == 0 - # and ((auxiliary_right_count > 0)) - # ) - # ): - if jid in bad_turn_config_jids: + USE_TURN_CONFIG = False auxiliary_around_count = ( default_auxiliary_around_count if auxiliary_around_count == 0 @@ -2712,6 +2686,92 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): auxiliary_main_out_start = default_auxiliary_main_out_start auxiliary_main_out_end = default_auxiliary_main_out_end # connect lanes + MAIN_HAS_STRAIGHT_LANES = False + AUXILIARY_HAS_STRAIGHT_LANES = False + if main_count > 0 and out_main_group is not None: + assert ( + out_main_group is not None + ), f"Junction {jid} arranged straight out_ways but no out_ways to connect!" + out_angle, out_way_ids = out_main_group + ( + out_main_wid, + out_auxiliary_wids, + ) = classify_main_auxiliary_wid( + out_way_ids, out_angle, "out_ways" + ) + out_main_lanes = self.map_roads[out_main_wid]["lanes"] + out_auxiliary_lanes = [ + l + for wid in out_auxiliary_wids + for l in self.map_roads[wid]["lanes"] + ] + # main road to main road + if ( + len(in_main_lanes[main_out_start:main_out_end]) + ) > 0 and len(out_main_lanes) > 0: + MAIN_HAS_STRAIGHT_LANES = True + junc_lanes += self._connect_lane_group( + in_lanes=in_main_lanes[main_out_start:main_out_end], + out_lanes=out_main_lanes, + lane_turn=mapv2.LANE_TURN_STRAIGHT, + lane_type=mapv2.LANE_TYPE_DRIVING, + junc_id=j_uid, + ) + # auxiliary road to auxiliary road + if ( + len( + in_auxiliary_lanes[ + auxiliary_main_out_start:auxiliary_main_out_end + ] + ) + > 0 + and len(out_auxiliary_lanes) > 0 + ): + AUXILIARY_HAS_STRAIGHT_LANES = True + junc_lanes += self._connect_lane_group( + in_lanes=in_auxiliary_lanes[ + auxiliary_main_out_start:auxiliary_main_out_end + ], + out_lanes=out_auxiliary_lanes, + lane_turn=mapv2.LANE_TURN_STRAIGHT, + lane_type=mapv2.LANE_TYPE_DRIVING, + junc_id=j_uid, + ) + # When there is no in auxiliary road, connect the in main road to the out auxiliary road + if ( + not len(in_auxiliary_lanes) > 0 + and len(out_auxiliary_lanes) > 0 + and len(in_main_lanes[main_out_start:main_out_end]) > 0 + ): + MAIN_HAS_STRAIGHT_LANES = True + junc_lanes += self._connect_lane_group( + in_lanes=in_main_lanes[main_out_start:main_out_end], + out_lanes=out_auxiliary_lanes, + lane_turn=mapv2.LANE_TURN_STRAIGHT, + lane_type=mapv2.LANE_TYPE_DRIVING, + junc_id=j_uid, + ) + # When there is no out auxiliary road, connect the in auxiliary road to the main road + if ( + len( + in_auxiliary_lanes[ + auxiliary_main_out_start:auxiliary_main_out_end + ] + ) + > 0 + and not len(out_auxiliary_lanes) > 0 + and len(out_main_lanes) > 0 + ): + AUXILIARY_HAS_STRAIGHT_LANES = True + junc_lanes += self._connect_lane_group( + in_lanes=in_auxiliary_lanes[ + auxiliary_main_out_start:auxiliary_main_out_end + ], + out_lanes=out_main_lanes, + lane_turn=mapv2.LANE_TURN_STRAIGHT, + lane_type=mapv2.LANE_TYPE_DRIVING, + junc_id=j_uid, + ) if around_count > 0: if key == (1, 1): for out_angle, out_way_ids in out_around_groups: @@ -2885,6 +2945,17 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): lane_type=mapv2.LANE_TYPE_DRIVING, junc_id=j_uid, ) + # ATTENTION: add # of left lanes when no turn configs provided + if not USE_TURN_CONFIG: + if not MAIN_HAS_STRAIGHT_LANES: + left_count = max( + left_count, len(in_main_lanes) - right_count + ) + if not AUXILIARY_HAS_STRAIGHT_LANES: + auxiliary_left_count = max( + auxiliary_left_count, + len(in_auxiliary_lanes) - auxiliary_right_count, + ) if left_count > 0: for out_angle, out_way_ids in out_left_groups: ( @@ -2943,86 +3014,7 @@ def classify_main_auxiliary_wid(wids, way_angle, group_type): lane_type=mapv2.LANE_TYPE_DRIVING, junc_id=j_uid, ) - if main_count > 0 and out_main_group is not None: - assert ( - out_main_group is not None - ), f"Junction {jid} arranged straight out_ways but no out_ways to connect!" - out_angle, out_way_ids = out_main_group - ( - out_main_wid, - out_auxiliary_wids, - ) = classify_main_auxiliary_wid( - out_way_ids, out_angle, "out_ways" - ) - out_main_lanes = self.map_roads[out_main_wid]["lanes"] - out_auxiliary_lanes = [ - l - for wid in out_auxiliary_wids - for l in self.map_roads[wid]["lanes"] - ] - # main road to main road - if ( - len(in_main_lanes[main_out_start:main_out_end]) - ) > 0 and len(out_main_lanes) > 0: - junc_lanes += self._connect_lane_group( - in_lanes=in_main_lanes[main_out_start:main_out_end], - out_lanes=out_main_lanes, - lane_turn=mapv2.LANE_TURN_STRAIGHT, - lane_type=mapv2.LANE_TYPE_DRIVING, - junc_id=j_uid, - ) - # auxiliary road to auxiliary road - if ( - len( - in_auxiliary_lanes[ - auxiliary_main_out_start:auxiliary_main_out_end - ] - ) - > 0 - and len(out_auxiliary_lanes) > 0 - ): - junc_lanes += self._connect_lane_group( - in_lanes=in_auxiliary_lanes[ - auxiliary_main_out_start:auxiliary_main_out_end - ], - out_lanes=out_auxiliary_lanes, - lane_turn=mapv2.LANE_TURN_STRAIGHT, - lane_type=mapv2.LANE_TYPE_DRIVING, - junc_id=j_uid, - ) - # When there is no in auxiliary road, connect the in main road to the out auxiliary road - if ( - not len(in_auxiliary_lanes) > 0 - and len(out_auxiliary_lanes) > 0 - and len(in_main_lanes[main_out_start:main_out_end]) > 0 - ): - junc_lanes += self._connect_lane_group( - in_lanes=in_main_lanes[main_out_start:main_out_end], - out_lanes=out_auxiliary_lanes, - lane_turn=mapv2.LANE_TURN_STRAIGHT, - lane_type=mapv2.LANE_TYPE_DRIVING, - junc_id=j_uid, - ) - # When there is no out auxiliary road, connect the in auxiliary road to the main road - if ( - len( - in_auxiliary_lanes[ - auxiliary_main_out_start:auxiliary_main_out_end - ] - ) - > 0 - and not len(out_auxiliary_lanes) > 0 - and len(out_main_lanes) > 0 - ): - junc_lanes += self._connect_lane_group( - in_lanes=in_auxiliary_lanes[ - auxiliary_main_out_start:auxiliary_main_out_end - ], - out_lanes=out_main_lanes, - lane_turn=mapv2.LANE_TURN_STRAIGHT, - lane_type=mapv2.LANE_TYPE_DRIVING, - junc_id=j_uid, - ) + # Check whether every road in the junction is connected valid_turn_config = True for _, in_way_ids in in_way_groups: diff --git a/mosstool/trip/generator/_util/utils.py b/mosstool/trip/generator/_util/utils.py index 006def8..fea5c2b 100644 --- a/mosstool/trip/generator/_util/utils.py +++ b/mosstool/trip/generator/_util/utils.py @@ -97,3 +97,136 @@ def recalculate_trip_mode_prob(profile: dict, V: np.ndarray): Filter some invalid trip modes according to the PersonProfile """ return V + + +def gen_departure_times( + rng: np.random.Generator, n_trip: int, mode: str, departure_prob +): + times = [0.0 for _ in range(n_trip)] + if departure_prob is not None: + LEN_TIMES = len(departure_prob) + TIME_INTERVAL = 24 / LEN_TIMES # (hour) + for i in range(n_trip): + time_center = ( + rng.choice(range(LEN_TIMES), p=departure_prob) + ) * TIME_INTERVAL + times[i] = rng.uniform( + time_center - 0.5 * TIME_INTERVAL, time_center + 0.5 * TIME_INTERVAL + ) + else: + # Defaults + ## Go to work + if mode[1] == "W" or mode[1] == "S": + times[0] = rng.normal(8, 2) + else: + times[0] = rng.normal(10.5, 2.5) + ## Go home + end_idx = -1 + for i in range(len(mode)): + if mode[len(mode) - i - 1] == "H": + if mode[len(mode) - i - 2] == "W" or mode[len(mode) - i - 2] == "S": + t = rng.normal(17, 2) + else: + t = rng.normal(14.5, 3) + + times[len(mode) - i - 2] = t + end_idx = len(mode) - i - 2 + break + if times[0] > times[end_idx]: + times[0], times[end_idx] = times[end_idx], times[0] + ## uniform distribute from work and home + for i in range(1, end_idx): + if times[end_idx] - times[0] > 4: + times[i] = rng.uniform(times[0] + 1, times[end_idx] - 1) + elif times[end_idx] - times[0] > 2: + times[i] = rng.uniform(times[0] + 0.5, times[end_idx] - 0.5) + else: + times[i] = rng.uniform(times[0], times[end_idx]) + if mode[-1] != "H": + sleep_time = 22 + rng.exponential(2) + if sleep_time - times[end_idx] > 3: + times[-2:] = rng.uniform(times[end_idx] + 0.5, sleep_time - 0.5, size=2) + elif sleep_time > times[end_idx]: + times[-2:] = rng.uniform(times[end_idx], sleep_time, size=2) + else: + times[-2:] = rng.uniform(sleep_time, times[end_idx], size=2) + times = np.sort(times) + times = times % 24 + return np.array(times) + + +def extract_HWEO_from_od_matrix( + aois: list, + n_region: int, + aoi2region: dict, + aoi_type2ids: dict, + od_prob: np.ndarray, + od_times_length: int, +): + total_popu = np.zeros(n_region) + work_popu = np.zeros(n_region) + educate_popu = np.zeros(n_region) + home_popu = np.zeros(n_region) + LEN_OD_TIMES = od_times_length + for aoi in aois: + # aoi type identify + external = aoi["external"] + aoi_id = aoi["id"] + if external["catg"] in HOME_CATGS: + aoi_type2ids["home"].append(aoi_id) + elif external["catg"] in WORK_CATGS: + aoi_type2ids["work"].append(aoi_id) + elif external["catg"] in EDUCATION_CATGS: + aoi_type2ids["education"].append(aoi_id) + else: + aoi_type2ids["other"].append(aoi_id) + # pop calculation + if aoi_id not in aoi2region: + continue + reg_idx = aoi2region[aoi_id] + total_popu[reg_idx] += external["population"] + work_popu[reg_idx] += ( + external["population"] if external["catg"] in WORK_CATGS else 0 + ) + home_popu[reg_idx] += ( + external["population"] if external["catg"] in HOME_CATGS else 0 + ) + educate_popu[reg_idx] += ( + external["population"] if external["catg"] in EDUCATION_CATGS else 0 + ) + home_dist = home_popu / sum(home_popu) + # initialization + work_od = np.zeros((n_region, n_region, LEN_OD_TIMES)) + educate_od = np.zeros((n_region, n_region, LEN_OD_TIMES)) + other_od = np.zeros((n_region, n_region, LEN_OD_TIMES)) + # calculate frequency + ## work + work_od[:, :, :] = od_prob[:, :, :] * work_popu[:, None] / total_popu[:, None] + work_od[:, total_popu <= 0, :] = 0 + ## study + educate_od[:, :, :] = od_prob[:, :, :] * educate_popu[:, None] / total_popu[:, None] + educate_od[:, total_popu <= 0, :] = 0 + ## other + other_od[:, :, :] = ( + od_prob[:, :, :] + * (total_popu[:, None] - work_popu[:, None] - educate_popu[:, None]) + / total_popu[:, None] + ) + other_od[:, total_popu <= 0, :] = 0 + # to probabilities + ## work + sum_work_od_j = np.sum(work_od, axis=1) + work_od = work_od / sum_work_od_j[:, np.newaxis] + work_od = np.nan_to_num(work_od) + work_od[np.where(sum_work_od_j == 0)] = 1 + ## study + sum_educate_od_j = np.sum(educate_od, axis=1) + educate_od = educate_od / sum_educate_od_j[:, np.newaxis] + educate_od = np.nan_to_num(educate_od) + educate_od[np.where(sum_educate_od_j == 0)] = 1 + ## other + sum_other_od_j = np.sum(other_od, axis=1) + other_od = other_od / sum_other_od_j[:, np.newaxis] + other_od = np.nan_to_num(other_od) + other_od[np.where(sum_other_od_j == 0)] = 1 + return home_dist, work_od, educate_od, other_od diff --git a/mosstool/trip/generator/generate_from_od.py b/mosstool/trip/generator/generate_from_od.py index df1e92b..5e4b77d 100644 --- a/mosstool/trip/generator/generate_from_od.py +++ b/mosstool/trip/generator/generate_from_od.py @@ -1,5 +1,6 @@ import logging from collections import defaultdict +from copy import deepcopy from functools import partial from math import ceil from multiprocessing import Pool, cpu_count @@ -21,10 +22,11 @@ from ...map._map_util.aoiutils import geo_coords from ...map._map_util.const import * -from ...util.format_converter import dict2pb +from ...util.format_converter import dict2pb, pb2dict from ...util.geo_match_pop import geo2pop from ._util.const import * -from ._util.utils import gen_profiles, recalculate_trip_mode_prob +from ._util.utils import (extract_HWEO_from_od_matrix, gen_departure_times, + gen_profiles, recalculate_trip_mode_prob) from .template import DEFAULT_PERSON @@ -130,64 +132,13 @@ def choose_aoi_with_type( home_region = rng.choice(n_region, p=H) else: home_region = a_home_region - mode = rng.choice(len(modes), p=p_mode) - mode = modes[mode] + mode_idx = rng.choice(len(modes), p=p_mode) + mode = modes[mode_idx] aoi_list = [] now_region = home_region # arrange time times[i]: departure time from region i n_trip = len(mode) - 1 if mode[-1] == "H" else len(mode) - times = [0.0 for _ in range(n_trip)] - if departure_prob is not None: - LEN_TIMES = len(departure_prob) - TIME_INTERVAL = 24 / LEN_TIMES # (hour) - for i in range(n_trip): - time_center = ( - rng.choice(range(LEN_TIMES), p=departure_prob) - ) * TIME_INTERVAL - times[i] = rng.uniform( - time_center - 0.5 * TIME_INTERVAL, time_center + 0.5 * TIME_INTERVAL - ) - else: - # Defaults - ## Go to work - if mode[1] == "W" or mode[1] == "S": - times[0] = rng.normal(8, 2) - else: - times[0] = rng.normal(10.5, 2.5) - ## Go home - end_idx = -1 - for i in range(len(mode)): - if mode[len(mode) - i - 1] == "H": - if mode[len(mode) - i - 2] == "W" or mode[len(mode) - i - 2] == "S": - t = rng.normal(17, 2) - else: - t = rng.normal(14.5, 3) - - times[len(mode) - i - 2] = t - end_idx = len(mode) - i - 2 - break - if times[0] > times[end_idx]: - times[0], times[end_idx] = times[end_idx], times[0] - ## uniform distribute from work and home - for i in range(1, end_idx): - if times[end_idx] - times[0] > 4: - times[i] = rng.uniform(times[0] + 1, times[end_idx] - 1) - elif times[end_idx] - times[0] > 2: - times[i] = rng.uniform(times[0] + 0.5, times[end_idx] - 0.5) - else: - times[i] = rng.uniform(times[0], times[end_idx]) - if mode[-1] != "H": - sleep_time = 22 + rng.exponential(2) - if sleep_time - times[end_idx] > 3: - times[-2:] = rng.uniform(times[end_idx] + 0.5, sleep_time - 0.5, size=2) - elif sleep_time > times[end_idx]: - times[-2:] = rng.uniform(times[end_idx], sleep_time, size=2) - else: - times[-2:] = rng.uniform(sleep_time, times[end_idx], size=2) - - times = np.sort(times) - times = times % 24 - + times = gen_departure_times(rng, n_trip, mode, departure_prob) for t in times: assert t >= 0 and t <= 24 @@ -301,6 +252,138 @@ def _process_agent_unit(args, d): ) +def _fill_sch_unit( + O, + p_home, + p_home_region, + p_work, + p_work_region, + p_profile, + modes, + p_mode, + seed, + args, +): + global region2aoi, aoi_map, aoi_type2ids + n_region, projector, departure_prob, LEN_OD_TIMES = args + OD_TIME_INTERVAL = 24 / LEN_OD_TIMES # (hour) + + def choose_aoi_with_type( + region_id, + aoi_type: Union[ + Literal["work"], + Literal["home"], + Literal["education"], + Literal["other"], + ], + ): + region_aois = region2aoi[region_id] + if len(region_aois) > 0: + if len(region2aoi[region_id]) == 1: + return region2aoi[region_id][0] + popu = np.zeros(len(region_aois)) + for i, id in enumerate(region_aois): + popu[i] = aoi_map[id]["external"]["population"] + if sum(popu) == 0: + idx = rng.choice(len(region_aois)) + return region_aois[idx] + p = popu / sum(popu) + idx = rng.choice(len(region_aois), p=p) + return region_aois[idx] + else: + return rng.choice(aoi_type2ids[aoi_type]) + + rng = np.random.default_rng(seed) + mode_idx = rng.choice(len(modes), p=p_mode) + mode = modes[mode_idx] + aoi_list = [] + now_region = p_home_region + # arrange time times[i]: departure time from region i + n_trip = len(mode) - 1 if mode[-1] == "H" else len(mode) + times = gen_departure_times(rng, n_trip, mode, departure_prob) + for t in times: + assert t >= 0 and t <= 24 + + # add hidden H for further process + # HWH+ ->HWH+H + if mode[-1] != "H": + mode = mode + "H" + # home aoi + home_aoi = p_home + home_region = p_home_region + # work aoi + work_aoi = p_work + work_region = p_work_region + # add aois + for idx, mode_i in enumerate(mode): + if mode_i == "H": + aoi_list.append(home_aoi) + now_region = home_region + elif mode_i == "W" or mode_i == "S": + aoi_list.append(work_aoi) + now_region = work_region + elif mode_i == "+" or mode_i == "O": + p = np.zeros((n_region)) + p = O[ + now_region, + :, + min(int(times[idx - 1] / OD_TIME_INTERVAL), LEN_OD_TIMES - 1), + ] # hour to int index + p = p / sum(p) + other_region = rng.choice(n_region, p=p) + other_aoi = choose_aoi_with_type(other_region, "other") + aoi_list.append(other_aoi) + now_region = other_region + trip_modes = [] + for cur_aoi, next_aoi in zip(aoi_list[:-1], aoi_list[1:]): + lon1, lat1 = aoi_map[cur_aoi]["geo"][0][:2] + lon2, lat2 = aoi_map[next_aoi]["geo"][0][:2] + p1 = projector(longitude=lon1, latitude=lat1) + p2 = projector(longitude=lon2, latitude=lat2) + trip_modes.append(_get_mode_with_distribution(p1, p2, p_profile, seed)) + + # determine activity + activities = [] + for next_aoi in aoi_list[1:]: + activities.append(aoi_map[next_aoi]["external"]["catg"]) + assert len(aoi_list) == len(times) + 1 + trip_models = ["" for _ in range(n_trip)] + return ( + aoi_list, + times, + trip_modes, + trip_models, + activities, + ) + + +def _fill_person_schedule_unit(args, d): + global other_od + ( + p_id, + p_home, + p_home_region, + p_work, + p_work_region, + p_profile, + modes, + p_mode, + seed, + ) = d + return _fill_sch_unit( + other_od, + p_home, + p_home_region, + p_work, + p_work_region, + p_profile, + modes, + p_mode, + seed, + args, + ) + + __all__ = ["TripGenerator"] @@ -361,10 +444,6 @@ def __init__( self.template = template self.template.ClearField("schedules") self.template.ClearField("home") - self.area_shapes = [] - self.aoi2region = defaultdict(list) - self.region2aoi = defaultdict(list) - self.aoi_type2ids = defaultdict(list) self.persons = [] # activity proportion if activity_distributions is not None: @@ -469,7 +548,7 @@ def _read_regions(self): self.regions.append(r) def _read_od_matrix(self): - logging.info("reading original ods") + logging.info("Reading original ods") n_region = len(self.regions) assert ( n_region == self.od_matrix.shape[0] and n_region == self.od_matrix.shape[1] @@ -509,7 +588,7 @@ def _match_aoi2region(self): aoi_id, reg_id = r[:2] self.aoi2region[aoi_id] = reg_id self.region2aoi[reg_id].append(aoi_id) - logging.info(f"aoi_matched: {len(self.aoi2region)}") + logging.info(f"AOI matched: {len(self.aoi2region)}") def _generate_mobi( self, @@ -523,77 +602,14 @@ def _generate_mobi( region2aoi = self.region2aoi aoi_map = {d["id"]: d for d in self.aois} n_region = len(self.regions) - total_popu = np.zeros(n_region) - work_popu = np.zeros(n_region) - educate_popu = np.zeros(n_region) - home_popu = np.zeros(n_region) - for aoi in self.aois: - # aoi type identify - external = aoi["external"] - aoi_id = aoi["id"] - if external["catg"] in HOME_CATGS: - self.aoi_type2ids["home"].append(aoi_id) - elif external["catg"] in WORK_CATGS: - self.aoi_type2ids["work"].append(aoi_id) - elif external["catg"] in EDUCATION_CATGS: - self.aoi_type2ids["education"].append(aoi_id) - else: - self.aoi_type2ids["other"].append(aoi_id) - # pop calculation - if aoi_id not in self.aoi2region: - continue - reg_idx = self.aoi2region[aoi_id] - total_popu[reg_idx] += external["population"] - work_popu[reg_idx] += ( - external["population"] if external["catg"] in WORK_CATGS else 0 - ) - home_popu[reg_idx] += ( - external["population"] if external["catg"] in HOME_CATGS else 0 - ) - educate_popu[reg_idx] += ( - external["population"] if external["catg"] in EDUCATION_CATGS else 0 - ) - home_dist = home_popu / sum(home_popu) - # initialization - work_od = np.zeros((n_region, n_region, self.LEN_OD_TIMES)) - educate_od = np.zeros((n_region, n_region, self.LEN_OD_TIMES)) - other_od = np.zeros((n_region, n_region, self.LEN_OD_TIMES)) - # calculate frequency - ## work - work_od[:, :, :] = ( - self.od_prob[:, :, :] * work_popu[:, None] / total_popu[:, None] - ) - work_od[:, total_popu <= 0, :] = 0 - ## study - educate_od[:, :, :] = ( - self.od_prob[:, :, :] * educate_popu[:, None] / total_popu[:, None] - ) - educate_od[:, total_popu <= 0, :] = 0 - ## other - other_od[:, :, :] = ( - self.od_prob[:, :, :] - * (total_popu[:, None] - work_popu[:, None] - educate_popu[:, None]) - / total_popu[:, None] + home_dist, work_od, educate_od, other_od = extract_HWEO_from_od_matrix( + self.aois, + n_region, + self.aoi2region, + self.aoi_type2ids, + self.od_prob, + self.LEN_OD_TIMES, ) - other_od[:, total_popu <= 0, :] = 0 - # to probabilities - ## work - sum_work_od_j = np.sum(work_od, axis=1) - work_od = work_od / sum_work_od_j[:, np.newaxis] - work_od = np.nan_to_num(work_od) - work_od[np.where(sum_work_od_j == 0)] = 1 - ## study - sum_educate_od_j = np.sum(educate_od, axis=1) - educate_od = educate_od / sum_educate_od_j[:, np.newaxis] - educate_od = np.nan_to_num(educate_od) - educate_od[np.where(sum_educate_od_j == 0)] = 1 - ## other - sum_other_od_j = np.sum(other_od, axis=1) - other_od = other_od / sum_other_od_j[:, np.newaxis] - other_od = np.nan_to_num(other_od) - other_od[np.where(sum_other_od_j == 0)] = 1 - # global variables - ## home_dist, work_od, educate_od, other_od = home_dist, work_od, educate_od, other_od aoi_type2ids = self.aoi_type2ids agent_args = [] a_home_regions = [] @@ -700,6 +716,11 @@ def generate_persons( Returns: - List[Person]: The generated person objects. """ + # init + self.area_shapes = [] + self.aoi2region = defaultdict(list) + self.region2aoi = defaultdict(list) + self.aoi_type2ids = defaultdict(list) self.persons = [] # user input time curve if departure_time_curve is not None: @@ -867,3 +888,140 @@ def _transfer_conn_road_ids( self.persons.append(p) agent_id += 1 return self.persons + + def _generate_schedules(self, input_persons: List[Person], seed: int): + global region2aoi, aoi_map, aoi_type2ids + global home_dist, work_od, other_od, educate_od + aoi_map = {d["id"]: d for d in self.aois} + n_region = len(self.regions) + home_dist, work_od, educate_od, other_od = extract_HWEO_from_od_matrix( + self.aois, + n_region, + self.aoi2region, + self.aoi_type2ids, + self.od_prob, + self.LEN_OD_TIMES, + ) + orig_persons = deepcopy(input_persons) + person_args = [] + bad_person_indexes = set() + rng = np.random.default_rng(seed) + for idx, p in enumerate(orig_persons): + try: + p_home = p.home.aoi_position.aoi_id + p_work = p.work.aoi_position.aoi_id + p_profile = pb2dict(p.profile) + person_args.append( + [ + p.id, + p_home, + self.aoi2region[p_home], # possibly key error + p_work, + self.aoi2region[p_work], # possibly key error + p_profile, + self.modes, + self.p_mode, + rng.integers(0, 2**16 - 1), + ] + ) + except Exception as e: + bad_person_indexes.add(idx) + logging.warning(f"{e} when handling Person {p.id}, Skip!") + to_process_persons = [ + p for idx, p in enumerate(orig_persons) if idx not in bad_person_indexes + ] + # return directly + no_process_persons = [ + p for idx, p in enumerate(orig_persons) if idx in bad_person_indexes + ] + partial_args = ( + len(self.regions), + self.projector, + self.departure_prob, + self.LEN_OD_TIMES, + ) + filled_schedules = [] + _fill_person_schedule_unit_with_arg = partial( + _fill_person_schedule_unit, partial_args + ) + for i in range(0, len(person_args), MAX_BATCH_SIZE): + person_args_batch = person_args[i : i + MAX_BATCH_SIZE] + with Pool(processes=self.workers) as pool: + filled_schedules += pool.map( + _fill_person_schedule_unit_with_arg, + person_args_batch, + chunksize=min(ceil(len(person_args_batch) / self.workers), 500), + ) + for ( + aoi_list, + times, + trip_modes, + trip_models, + activities, + ), orig_p in zip(filled_schedules, to_process_persons): + times = np.array(times) * 3600 # hour->second + p = Person() + p.CopyFrom(orig_p) + p.ClearField("schedules") + for time, aoi_id, trip_mode, activity, trip_model in zip( + times, aoi_list[1:], trip_modes, activities, trip_models + ): + schedule = cast(Schedule, p.schedules.add()) + schedule.departure_time = time + schedule.loop_count = 1 + trip = Trip( + mode=cast( + TripMode, + trip_mode, + ), + end=Position(aoi_position=AoiPosition(aoi_id=aoi_id)), + activity=activity, + model=trip_model, + ) + schedule.trips.append(trip) + self.persons.append(p) + len_filled_persons, len_no_process_persons = len(to_process_persons), len( + no_process_persons + ) + logging.info(f"Filled schedules of {len_filled_persons} persons") + if len_no_process_persons > 0: + logging.warning( + f"Unprocessed persons: {len_no_process_persons}, index range [{len_filled_persons}:] in returned results" + ) + self.persons.extend(no_process_persons) + + def fill_person_schedules( + self, + input_persons: List[Person], + od_matrix: np.ndarray, + areas: GeoDataFrame, + departure_time_curve: Optional[list[float]] = None, + seed: int = 0, + ) -> List[Person]: + """ + Generate person schedules. + """ + # init + self.area_shapes = [] + self.aoi2region = defaultdict(list) + self.region2aoi = defaultdict(list) + self.aoi_type2ids = defaultdict(list) + self.persons = [] + # user input time curve + if departure_time_curve is not None: + assert len(departure_time_curve) >= 24 + sum_times = sum(departure_time_curve) + self.departure_prob = np.array( + [d / sum_times for d in departure_time_curve] + ) + else: + self.departure_prob = None + self.od_matrix = od_matrix + self.areas = areas + self._read_aois() + self._read_regions() + self._read_od_matrix() + self._match_aoi2region() + self._generate_schedules(input_persons, seed) + + return self.persons diff --git a/pyproject.toml b/pyproject.toml index 247d320..4a3e191 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "mosstool" -version = "0.2.21" +version = "0.2.22" description = "MObility Simulation System toolbox " authors = ["Jun Zhang "] license = "MIT"