diff --git a/activitysim/abm/models/initialize_from_usim.py b/activitysim/abm/models/initialize_from_usim.py index cef9dbf77a..f6532852f2 100644 --- a/activitysim/abm/models/initialize_from_usim.py +++ b/activitysim/abm/models/initialize_from_usim.py @@ -1,7 +1,6 @@ import os import numpy as np import pandas as pd -import os import matplotlib.pyplot as plt import geopandas as gpd import orca @@ -10,6 +9,7 @@ from urbansim.utils import misc import requests import openmatrix as omx +from shapely import wkt import logging from activitysim.core import config @@ -29,6 +29,62 @@ def get_zone_geoms_from_h3(h3_ids): return polygon_shapes +def assign_taz(df, gdf): + ''' + Assigns the gdf index (TAZ ID) for each index in df + Input: + - df columns names x, and y. The index is the ID of the object(blocks, school, college) + - gdf: Geopandas DataFrame with TAZ as index, geometry and area value. + Output: + A series with df index and corresponding gdf id + ''' + + df = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.x, df.y), crs = "EPSG:4326") + gdf.geometry.crs = "EPSG:4326" + + assert df.geometry.crs == gdf.geometry.crs + + # Spatial join + df = gpd.sjoin(df, gdf, how = 'left', op = 'intersects') + + #Drop duplicates and keep the one with the smallest H3 area + df = df.sort_values('area') + index_name = df.index.name + df.reset_index(inplace = True) + df.drop_duplicates(subset = [index_name], keep = 'first', inplace = True) + df.set_index(index_name, inplace = True) + + #Check if there is any assigined object + if df.index_right.isnull().sum()>0: + + #Buffer unassigned ids until they reach a hexbin. + null_values = df[df.index_right.isnull()].drop(columns = ['index_right','area']) + + result_list = [] + for index, value in null_values.iterrows(): + buff_size = 0.0001 + matched = False + geo_value = gpd.GeoDataFrame(value).T + geo_value.crs = "EPSG:4326" + while matched == False: + geo_value.geometry = geo_value.geometry.buffer(buff_size) + result = gpd.sjoin(geo_value, gdf, how = 'left', op = 'intersects') + matched = ~result.index_right.isnull()[0] + buff_size = buff_size + 0.0001 + result_list.append(result.iloc[0:1]) + + null_values = pd.concat(result_list) + + # Concatenate newly assigned values to the main values table + df = df.dropna() + df = pd.concat([df, null_values], axis = 0) + + return df.index_right + + else: + return df.index_right + + # ** 1. CREATE NEW TABLES ** # Zones @@ -92,6 +148,7 @@ def schools(blocks): enrollment = enrollment[[ 'ncessch', 'county_code', 'latitude', 'longitude', 'enrollment']].set_index('ncessch') + enrollment.rename(columns = {'longitude':'x', 'latitude':'y'}, inplace = True) return enrollment.dropna() @@ -130,46 +187,9 @@ def colleges(blocks): @orca.column('blocks', cache = True) def TAZ(blocks, zones): - - # Tranform blocks to a Geopandas dataframe - blocks_df = blocks.to_frame(columns=['x', 'y']) - zones_df = zones.to_frame(columns=['geometry', 'area']) - h3_gpd = gpd.GeoDataFrame(zones_df, crs='EPSG:4326') - - blocks_df = gpd.GeoDataFrame( - blocks_df, geometry=gpd.points_from_xy(blocks_df.x, blocks_df.y), - crs="EPSG:4326") - - # Spatial join - blocks_df = gpd.sjoin(blocks_df, h3_gpd, how='left', op = 'intersects') - - # Drop duplicates and keep the one with the smallest H3 area - blocks_df = blocks_df.sort_values('area') - blocks_df.drop_duplicates(subset = ['x', 'y'], keep = 'first', inplace = True) - - # Buffer unassigned blocks until they reach a hexbin. - null_blocks = blocks_df[blocks_df.index_right.isnull()].drop(columns = ['index_right','area']) - - result_list = [] - for index, block in null_blocks.iterrows(): - buff_size = 0.0001 - matched = False - geo_block = gpd.GeoDataFrame(block, crs='EPSG:4326').T - while matched == False: - geo_block.geometry = geo_block.geometry.buffer(buff_size) - result = gpd.sjoin(geo_block, h3_gpd, how = 'left', op = 'intersects') - matched = ~result.index_right.isnull()[0] - buff_size = buff_size + 0.0001 - result_list.append(result.iloc[0:1]) - - null_blocks = pd.concat(result_list) - - # Concatenate newly assigned blocks to the main blocks table - blocks_df = blocks_df.dropna() - blocks_df = pd.concat([blocks_df, null_blocks], axis = 0) - - return blocks_df.index_right - + blocks_df = blocks.to_frame(columns = ['x', 'y']) + h3_gpd = zones.to_frame(columns = ['geometry', 'area']) + return assign_taz(blocks_df, h3_gpd) @orca.column('blocks') def CI_employment(jobs, blocks): @@ -202,47 +222,9 @@ def RESACRE(blocks): @orca.column('schools', cache = True) def TAZ(schools, zones): - - #Tranform blocks to a Geopandas dataframe - zones_df = zones.to_frame(columns=['geometry', 'area']) - h3_gpd = gpd.GeoDataFrame(zones_df, crs='EPSG:4326') - - school_gpd = schools.to_frame(columns = ['ncessch','longitude', 'latitude']) - school_gpd = gpd.GeoDataFrame( - school_gpd, - geometry=gpd.points_from_xy(school_gpd.longitude, school_gpd.latitude), - crs="EPSG:4326") - # Spatial join - school_gdf = gpd.sjoin(school_gpd, h3_gpd, how = 'left', op = 'intersects') - - #Drop duplicates and keep the one with the smallest H3 area - school_gdf = school_gdf.sort_values('area') - school_gdf.reset_index(inplace = True) - school_gdf.drop_duplicates(subset = ['ncessch'], keep = 'first', inplace = True) - - #Buffer unassigned blocks until they reach a hexbin. - null_schools = school_gdf[school_gdf.index_right.isnull()].drop(columns = ['index_right','area']) - - result_list = [] - for index, school in null_schools.iterrows(): - buff_size = 0.0001 - matched = False - geo_school = gpd.GeoDataFrame(school, crs='EPSG:4326').T - while matched == False: - geo_school.geometry = geo_school.geometry.buffer(buff_size) - result = gpd.sjoin(geo_school, h3_gpd, how = 'left', op = 'intersects') - matched = ~result.index_right.isnull().iloc[0] - buff_size = buff_size + 0.0001 - result_list.append(result.iloc[0:1]) - - null_school = pd.concat(result_list) - - # Concatenate newly assigned blocks to the main blocks table - school_gdf = school_gdf.dropna() - school_all = pd.concat([school_gdf, null_school], axis = 0) - school_all.set_index('ncessch', inplace = True) - return school_all.index_right - + h3_gpd = zones.to_frame(columns = ['geometry', 'area']) + school_gpd = orca.get_table('schools').to_frame(columns = ['x', 'y']) + return assign_taz(school_gpd, h3_gpd) # Colleges Variables @@ -294,26 +276,11 @@ def part_time_enrollment(): return s -@orca.column('colleges', cache=True) +@orca.column('colleges', cache = True) def TAZ(colleges, zones): - #Tranform blocks to a Geopandas dataframe colleges_df = colleges.to_frame(columns = ['x', 'y']) - zones_df = zones.to_frame(columns = ['geometry', 'area']) - h3_gpd = gpd.GeoDataFrame(zones_df, crs="EPSG:4326") - - colleges_df = gpd.GeoDataFrame( - colleges_df, geometry=gpd.points_from_xy(colleges_df.x, colleges_df.y), - crs="EPSG:4326") - - # Spatial join - colleges_df = gpd.sjoin(colleges_df, h3_gpd, how = 'left', op = 'intersects') - - #Drop duplicates and keep the one with the smallest H3 area - colleges_df = colleges_df.sort_values('area') - colleges_df.drop_duplicates(subset = ['x', 'y'], keep = 'first', inplace = True) - - return colleges_df.index_right - + h3_gpd = zones.to_frame(columns = ['geometry', 'area']) + return assign_taz(colleges_df, h3_gpd) # Households Variables @@ -660,10 +627,17 @@ def COLLPTE(colleges, zones): @orca.column('zones') -def area_type(): - # Integer, 0=regional core, 1=central business district, 2=urban business, - # 3=urban, 4=suburban, 5=rural - return 0 # Assuming all regional core +def area_type(mpo_taz, zones): + + mpo = mpo_taz.to_frame(columns = ['geometry','area_type','ACRES']) + h3_gpd = zones.to_frame(columns = ['geometry', 'area']) + + join = gpd.sjoin(h3_gpd, mpo, how = 'left',op='intersects') + join.area_type.fillna(5, inplace = True) #Fill non-matched areas with 5 (Rural areas) + join = join.groupby(['TAZ', 'area_type'])['ACRES'].sum().reset_index() + join = join.sort_values(['TAZ', 'ACRES'], ascending = True) + s = join.groupby('TAZ')['area_type'].last() + return s @orca.column('zones') @@ -694,6 +668,8 @@ def load_usim_data(data_dir, settings): persons = hdf['/persons'] blocks = hdf['/blocks'] jobs = hdf['/jobs'] + mpo_taz = hdf['/mpo_taz'] + hdf.close() # add home x,y coords to persons table @@ -705,6 +681,10 @@ def load_usim_data(data_dir, settings): left_on='block_id', right_index=True) persons['home_x'] = persons_w_xy['x'] persons['home_y'] = persons_w_xy['y'] + + #Tranform mpo_taz to a geoDataFrame + mpo_taz['geometry'] = mpo_taz['geometry'].apply(wkt.loads) + mpo_taz = gpd.GeoDataFrame(mpo_taz, geometry='geometry', crs ='EPSG:4326') del persons_w_res_blk del persons_w_xy @@ -713,8 +693,8 @@ def load_usim_data(data_dir, settings): orca.add_table('usim_persons', persons) orca.add_table('blocks', blocks) orca.add_table('jobs', jobs) - - + orca.add_table('mpo_taz', mpo_taz) + # Export households tables @inject.step() def create_inputs_from_usim_data(data_dir): diff --git a/activitysim/abm/models/initialize_skims_from_beam.py b/activitysim/abm/models/initialize_skims_from_beam.py index 52083ff362..190a91d807 100644 --- a/activitysim/abm/models/initialize_skims_from_beam.py +++ b/activitysim/abm/models/initialize_skims_from_beam.py @@ -102,6 +102,7 @@ def create_skims_from_beam(raw_beam_skims, data_dir): # Adding car distance skims vals = auto_df[beam_asim_hwy_measure_map['DIST']].values mx = vals.reshape((num_taz, num_taz)) + skims['DIST'] = mx # active skims diff --git a/austin_mp/area_type_impute.py b/austin_mp/area_type_impute.py new file mode 100644 index 0000000000..902c54f9fe --- /dev/null +++ b/austin_mp/area_type_impute.py @@ -0,0 +1,42 @@ +import pandas as pd +import numpy as np +import geopandas as gpd +import h3 +import matplotlib.pyplot as plt +from shapely import wkt +import orca + +# ## Preprocessing the original MPO .shp files with TAZ +# Objective: Get area type as: +# - 0: Regional core +# - 1: CBD +# - 2: Urban Business +# - 3: Urban +# - 4: Suburban +# - 5: Rural + +#Load MPO TAZs shapefiles +mpo_taz = gpd.read_file('data/tazs_austin/2015_2045 CAMPO TAZ SHAPE.shp') +mpo_taz = mpo_taz.to_crs('EPSG:4326') + +#Transformation values: +mpo_taz = mpo_taz[mpo_taz.SMTDNAME != 'OutofArea'] +area_type_dict = {'CBD': 1, 'UrbIntTravis': 2, 'UrbTravis': 3, + 'SubTravis': 4, 'RurTravis':5,'UrbIntWilliamson': 2, + 'UrbWilliamson': 3, 'SubWilliamson': 4,'RurWilliamson': 5, + 'UrbIntHays': 2, 'UrbHays': 3, 'SubHays': 4, 'RurHays': 5, + 'UrbIntBastrop': 2, 'UrbBastrop': 3, 'SubBastrop': 4, + 'RurBastrop': 5, 'UrbCaldwell': 3, 'SubCaldwell': 4, + 'RurCaldwell': 5, 'UrbBurnet': 3,'SubBurnet':4, 'RurBurnet':5} + +mpo_taz['area_type'] = mpo_taz.SMTDNAME.replace(area_type_dict) + + +#Transform geopandas to dataframe +mpo_taz = pd.DataFrame(mpo_taz) +mpo_taz['geometry'] = mpo_taz.geometry.astype('str') + +#Save it to .H5 file +hdf = pd.HDFStore('model_data.h5') +hdf.append(key = 'mpo_taz', value = mpo_taz) +hdf.close() diff --git a/austin_mp/configs/settings.yaml b/austin_mp/configs/settings.yaml index 3df0dfd664..c10823ffcb 100644 --- a/austin_mp/configs/settings.yaml +++ b/austin_mp/configs/settings.yaml @@ -48,20 +48,20 @@ use_shadow_pricing: False ## - example sample -households_sample_size: 0 -chunk_size: 400000000 -num_processes: 24 -stagger: 2 +households_sample_size: 10000 +chunk_size: 5000000000 +num_processes: 20 +stagger: 0 # - tracing trace_hh_id: trace_od: -#trace_hh_id: 1482966 +trace_hh_id: 195809 #trace_od: [5, 11] # to resume after last successful checkpoint, specify resume_after: _ -# resume_after: trip_mode_choice +# resume_after: _ models: ### mp_initialize step @@ -71,7 +71,7 @@ models: ### mp_households step - school_location - workplace_location - - auto_ownership_simulate +# - auto_ownership_simulate - free_parking - cdap_simulate - mandatory_tour_frequency diff --git a/austin_mp/simulation_mp.py b/austin_mp/simulation_mp.py index 32977f1303..2400a2893b 100644 --- a/austin_mp/simulation_mp.py +++ b/austin_mp/simulation_mp.py @@ -15,7 +15,6 @@ from activitysim.core import mp_tasks from activitysim.core import chunk - logger = logging.getLogger('activitysim') @@ -34,12 +33,12 @@ def cleanup_output_files(): def run(run_list, injectables=None): - + # Create a new skims.omx file from BEAM (http://beam.lbl.gov/) skims # if skims do not already exist in the input data directory if config.setting('create_skims_from_beam'): pipeline.run(models=['create_skims_from_beam']) - pipeline.close_pipeline() + pipeline.close_pipeline() # Create persons, households, and land use .csv files from UrbanSim # data if these files do not already exist in the input data directory