survey modeshare code

survey/champ/HwySkimUtil.py

@@ -0,0 +1,54 @@
+'''
+Created on Jan 25, 2010
+
+@author: Lisa Zorn
+
+Generic trip record class, plus some extra functions that will likely come up.
+
+Modified in 2014 by Bhargava Sana to just isolate a skim query functionality.
+'''
+
+from tables import open_file
+import os
+
+# Functionally constants
+TIMEPERIODS     = { 1:"EA", 2:"AM", 3:"MD", 4:"PM", 5:"EV" }
+
+class SkimUtil2:
+    """
+    Helper class to read Skim files and lookup time/cost/distance for given O-D pairs.
+    This class is written for low-memory, not for speed.  So it'll take forever to go
+    through a trip file and do the skim lookups but you won't be hitting memory limits.
+    """
+
+    def __init__(self, skimdir, timeperiods=[2], skimprefix=""):
+
+        self.skimdir         = skimdir
+        self.hwyskims        = { 1:{}, 2:{}, 3:{}, 4:{}, 5:{} }
+        for tkey in timeperiods:
+            self.hwyskims[tkey] = open_file(os.path.join(skimdir,skimprefix+"HWYALL" + TIMEPERIODS[tkey] + ".h5"), mode="r")
+        self.termtime        = open_file(os.path.join(skimdir,"OPTERM.h5"), mode="r")
+
+        print("SkimUtil2 initialized for " + skimdir)
+
+    def getDASkims(self, otaz, dtaz, timeperiod=2):
+        """ Returns distance, time, out-of-pocket cost (fares, bridge & value tolls)
+            Units: miles, minutes, 1989 dollars.
+        """
+
+        (t,d,f) = (0,0,0)
+        termtime                = 0
+
+        # this is ok because of the PNR zones
+        if (otaz >= self.termtime.get_node('/', '1').shape[0] or
+            dtaz >= self.termtime.get_node('/', '1').shape[0]):
+            termtime = 0
+        else:
+            termtime = self.termtime.get_node('/', '1')[otaz-1][dtaz-1]
+
+        t = self.hwyskims[timeperiod].get_node('/', '1')[otaz-1,dtaz-1] + termtime
+        d = self.hwyskims[timeperiod].get_node('/', '2')[otaz-1,dtaz-1]
+        f = self.hwyskims[timeperiod].get_node('/', '3')[otaz-1,dtaz-1]/100.0
+        return (t,d,f)
+
+

survey/champ/Lookups.py

@@ -0,0 +1,127 @@
+# Initial revision 2011 Sept 14 by lmz
+# From Y:\champ\util\pythonlib\champUtil
+
+
+
+MAX_MTC_ZONE = 1475
+MAX_SF_ZONE = 2475
+MAX_SF_COUNTY_ZONE = 981
+
+
+class Lookups:
+    """
+    This class is just for those lookups that don't really belong anywhere else.
+    """
+
+    TIMEPERIODS_NUM_TO_STR  = {1:"EA", 2:"AM", 3:"MD", 4:"PM", 5:"EV" }
+    TIMEPERIODS_STR_TO_NUM  = dict((v,k) for k,v in TIMEPERIODS_NUM_TO_STR.items())
+    TIMEPERIOD_BINS = [0,180,360,750,930,1440]
+    TIMEPERIOD_LABELS = ['EA','AM','MD','PM','EV']
+    OPCOST                  = 0.12 # dollars/mile
+    TNC_FIXED               = 3.22 # dollars per trip
+    TNC_PER_MILE            = 1.04 # dollars/mile
+    WALKSPEED               = 3.0  # mph
+    BIKESPEED               = 10.0 # mph
+    PSEG                    = {1:"Worker", 2:"AdultStudent", 3:"Other", 4:"ChildStudent"}
+    PURPOSE_NUM_TO_STR      = {1:"Work",    2:"GradeSchool", 3:"HighSchool",
+                               4:"College", 5:"Other",       6:"WorkBased" }
+    PURPOSE_STR_TO_NUM      = dict((v,k) for k,v in PURPOSE_NUM_TO_STR.items())
+
+    #IMPORTANT - THIS ORDER IS "set" and shouldn't be changed unless changes are made to src/sftripmc/define.h
+    CHAMP_TRIP_MODES        =["DA",     "SR2",     "SR3",                         
+                              "DA_TOLL","SR2_TOLL","SR3_TOLL",
+                              "DA_PAID","SR2_PAID","SR3_PAID",
+                              "WALK","BIKE",
+                              "WLOC","WLRT","WPRE","WFER","WBAR",
+                              "DLOCW","DLRTW","DPREW","DFERW","DBARW",
+                              "TAXI",
+                              "WLOCD","WLRTD","WPRED","WFERD","WBARD",
+                              "TNC1","TNC2","TNC3"] # Add TNC by occupancy  
+    CHAMP_TRIP_MODES_NUM_TO_STR = dict(zip(list(range(1,len(CHAMP_TRIP_MODES)+1)), CHAMP_TRIP_MODES))
+    CHAMP_TRIP_MODES_STR_TO_NUM = dict(zip(CHAMP_TRIP_MODES, list(range(1,len(CHAMP_TRIP_MODES)+1))))
+
+
+    TRANSITMODES    = ["WLW", "ALW", "WLA",
+                       "WMW", "AMW", "WMA",
+                       "WPW", "APW", "WPA", 
+                       "WFW", "AFW", "WFA",
+                       "WBW", "ABW", "WBA"]
+
+    #IMPORTANT - THIS ORDER IS "set" and shouldn't be changed unless changes are made to src/sfchamp/define.h
+    CHAMP_TOUR_MODES        =["DA",     "SR2",      "SR3",
+                              "DA_TOLL","SR2_TOLL", "SR3_TOLL",
+                              "WALK",   "BIKE",
+                              "WTRN",   "DTRN",
+                              "TAXI"]
+    CHAMP_TOUR_MODES_NUM_TO_STR = dict(zip(list(range(1,len(CHAMP_TOUR_MODES)+1)), CHAMP_TOUR_MODES))
+    CHAMP_TOUR_MODES_STR_TO_NUM = dict(zip(CHAMP_TOUR_MODES, list(range(1,len(CHAMP_TOUR_MODES)+1))))
+
+    # from sftripmc/persdata.cpp
+    #                   //----------+-------------------------------------------------
+    #                   //          |                     TDDEPART
+    #                   // TODEPART |        1         2         3         4         5
+    #                   //----------+-------------------------------------------------
+    DURATION_TRIP = [
+                                    [      0.3,       1.2,     8.4,      10.5,    14.1],
+                                    [      1.2,       0.3,     4.8,       8.9,    11.5],
+                                    [      8.4,       4.8,     0.8,       2.7,     7.7],
+                                    [     10.5,       8.9,     2.7,       0.4,     2.0],
+                                    [     14.1,      11.5,     7.7,       2.0,     1.1] ]
+
+    # from sfourmc/persdata.cpp
+    DURATION_TOUR = [
+                                    [      0.3,       1.5,     8.2,      10.2,    13.1],
+                                    [      1.5,       0.4,     5.1,       8.7,    10.9],
+                                    [      8.2,       5.1,     1.0,       3.1,     7.6],
+                                    [     10.2,       8.7,     3.1,       0.5,     2.1],
+                                    [      2.4,       6.8,     9.4,      13.8,     1.4] ]
+
+    TIMEPERIODS_TO_SUBTIMES  = {"EA":[300, 500],
+                                "AM":[600, 630, 700, 730, 800, 830],
+                                "MD":[900, 1000, 1100, 130, 230,],
+                                "PM":[330, 400, 430, 500, 530, 600],
+                                "EV":[630, 730] }
+
+    TIMEPERIODS_TO_SUBTIME_DURATIONS = {"EA":[2.0, 1.0],
+                                        "AM":[0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
+                                        "MD":[1.0, 1.0, 2.5, 1.0, 1.0],
+                                        "PM":[0.5, 0.5, 0.5, 0.5, 0.5, 0.5],
+                                        "EV":[1.0, 7.5] }
+
+    @classmethod
+    def readSubTimeVolumeFactors(self):
+        """
+        Returns dict of dicts.  E.g.  { "EA":{300:0.405, 500:0.595}, ... }
+        """
+        import re
+        WHITESPACE_RE   = re.compile(r"^\s*$")
+        VOLFAC_RE       = re.compile(r"^\s*VOLFAC_(EA|AM|MD|PM|EV)([0-9]+)\s*=\s*([0-9\.]+)\s*$")
+        ret_dict        = {}
+
+        volfacfile = open('VolumeFactors.ctl', 'r')
+        for line in volfacfile:
+            # skip comments
+            if line[0] == ";": continue
+            # skip whitespace
+            if WHITESPACE_RE.match(line) != None: continue
+            match = VOLFAC_RE.match(line)
+            if match == None:
+                print("Do not understand line: [%s]" % line)
+                continue
+            timeperiod  = match.group(1)
+            subtime     = int(match.group(2))
+            volfac      = float(match.group(3))
+            if timeperiod not in ret_dict:
+                ret_dict[timeperiod] = {}
+            ret_dict[timeperiod][subtime] = volfac
+
+        volfacfile.close()
+        # verify they sum to 1 per main time period
+        for timeperiod in list(ret_dict.keys()):
+            total = 0.0
+            for subtime in list(ret_dict[timeperiod].keys()):
+                total += ret_dict[timeperiod][subtime]
+            if abs(total - 1.0) > 0.01:
+                print("Total for timeperiod %s is %f not 1.0: %s" % (timeperiod, total, str(ret_dict)))
+                exit(2)
+        return ret_dict

survey/champ/Skim.py

@@ -0,0 +1,50 @@
+import os
+from tables import open_file
+
+class SkimException(Exception):
+    """
+    This class is used to communicate Skim errors.
+    """
+    pass
+
+
+class Skim:
+    """
+    Base skim class.  Not sure what code will go in here or if it's just an API.
+    """
+
+    #: the time period codes
+    TIMEPERIOD_NUM_TO_STR = { 1:"EA", 2:"AM", 3:"MD", 4:"PM", 5:"EV" }
+
+    #: the purpose codes
+    PURPOSE_NUM_TO_STR = { 1:"Work",    2:"GradeSchool", 3:"HighSchool",
+                           4:"College", 5:"Other",       6:"WorkBased" }
+
+    def __init__(self, file_dir, file_names):
+        """
+        Opens the skim table file[s] in *file_dir*.  
+        *file_names* should be a list.
+        """
+
+        # mapping of filename -> skim file
+        self.skim_table_files = {}
+
+        for file_name in file_names:
+            full_file = os.path.join(file_dir, file_name)
+            if not os.path.exists(full_file):
+                raise SkimException("Skim: %s file doesn't exist" % full_file)
+
+            self.skim_table_files[file_name] = open_file(full_file, mode="r")
+
+    def __del__(self):
+        """
+        Closes the skim files
+        """
+        # print "Destructing Skim"
+        filenames = list(self.skim_table_files.keys())
+
+        for filename in filenames:
+            self.skim_table_files[filename].close()
+            del self.skim_table_files[filename]
+
+

survey/champ/TransitTourSkim.py

@@ -0,0 +1,78 @@
+from champ.Skim import Skim
+from tables import open_file
+
+class TransitTourSkim(Skim):
+    """
+    Transit Tour Skim class.
+    """
+
+    #: Matching of matrix in h5 file and attribute name
+    TABLE_NUMBER_TO_NAME = {1:"LIVT", 2:"RIVT", 3:"MIVT", 4:"PIVT", 5:"FIVT", 6:"BIVT", 
+                            7:"ACC_TIME", 8:"EGR_TIME", 9:"ACC_DIST", 10:"EGR_DIST",
+                            11:"IWAIT", 12:"XWAIT", 
+                            13:"TrDIST", 14:"DrDIST", 
+                            15:"FUNTIME", 16:"XWKTIME", 17:"NUM_LINKS", 18:"TOT_FARE",
+                            19:"ACCNODE", 20:"EGRNODE"}
+
+    #: All variables returned
+    ALL_VARS = list(TABLE_NUMBER_TO_NAME.values())
+    ALL_VARS.append("TOT_TIME")
+
+    #: Skims related to time (for converting hundredths of mins to mins)
+    TIME_SKIMS = ["LIVT", "RIVT", "MIVT", "PIVT", "FIVT", "BIVT",
+                  "ACC_TIME", "EGR_TIME", "IWAIT", "XWAIT", "FUNTIME", "XWKTIME"]
+    #: Skims related to distances (for converting hundredths of miles to miles)
+    DIST_SKIMS = ["TrDIST", "DrDIST", "ACC_DIST", "EGR_DIST"]
+    #: Skims related to cost
+    FARE_SKIMS = ["TOT_FARE"]
+
+    #: Tour Skim types, e.g. WTW, etc
+    TOUR_SKIM_TYPES = ["WTW"] # ["WTW", "ATW", "WTA"]
+
+    def __init__(self, file_dir, timeperiod="AM"):
+        """
+        Opens the given skim
+        """
+        self.timeperiod = timeperiod
+        self.trn_skim_files = list("TRN%s%s.h5" % (tourtype, timeperiod)
+                                   for tourtype in TransitTourSkim.TOUR_SKIM_TYPES)
+
+        Skim.__init__(self, file_dir, self.trn_skim_files)
+
+    def getTourAttributes(self, otaz, dtaz, tour_type="WTW"):
+        """
+        Returns a tuple of (time, distance, fare)
+
+        `tour_type` is one of :py:attr:`TransitTourSkim.TOUR_SKIM_TYPES`
+
+        Units are minutes, miles and 1989 cents.
+
+        Currently this only returns outbound OR return (not the sum) depending on how called.
+
+        A value for `TOT_TIME` is also included for convenience.
+
+        """
+        skim_file = "TRN%s%s.h5" % (tour_type, self.timeperiod)
+
+        transitAttributes = {}
+        tot_time = 0
+        tot_dist = 0
+
+        for tablenum,tablename in TransitTourSkim.TABLE_NUMBER_TO_NAME.items():
+            # convert hundredths of minutes to minutes
+            if tablename in TransitTourSkim.TIME_SKIMS:
+                transitAttributes[tablename] = 0.01 * self.skim_table_files[skim_file].root._f_get_child("%d" % tablenum)[otaz-1][dtaz-1]
+                tot_time += transitAttributes[tablename]
+
+            # convert hundredths of miles to miles
+            elif tablename in TransitTourSkim.DIST_SKIMS:
+                transitAttributes[tablename] = 0.01 * self.skim_table_files[skim_file].root._f_get_child("%d" % tablenum)[otaz-1][dtaz-1]
+                tot_dist += transitAttributes[tablename]
+            # FAREs are in the correct units already
+            else:
+                transitAttributes[tablename] = self.skim_table_files[skim_file].root._f_get_child("%d" % tablenum)[otaz-1][dtaz-1]
+
+        transitAttributes["TOT_TIME"] = tot_time
+
+        return (tot_time, tot_dist, transitAttributes['TOT_FARE'])
+

survey/champ/WalkSkim.py

@@ -0,0 +1,75 @@
+from champ.Skim import Skim
+from tables import open_file
+import numpy as np
+from champ.Lookups import MAX_SF_ZONE
+
+class WalkSkim(Skim):
+    r"""
+    Walk Skim class.
+
+    For now, test with Y:\champ\networks\RTP2009_CHAMP4.3plus\2000\hwy\addPedAttributes\walkSkim.h5
+    """
+
+    #: Maps the table number to the name of the skim table
+    TABLE_NUMBER_TO_NAME = { 
+        1 :"DISTANCE",     # link sum.  (miles)
+        2 :"INDIRECTNESS",   # distance divided by rock dove distance
+        3 :"RISE",         # link sum   (feet)
+        4 :"PER_RISE",     # Percent rise, or [rise / distance]
+        5 :"ABS_RISE",     # link sum when rise>0 (feet)
+        6 :"ABS_PER_RISE", # Percent rise [abs_rise / distance]
+        # the following are weighted by link distance
+        7:"AVGCAP",       # average road capacity (vph)
+        8:"AVGLANEAM",    # average road lanes
+        9:"AVGFFSPEED",   # average freeflow roadway speed
+        # the following are TAZ-based.  Also weighted by link distance
+        10:"AVGPOPDEN",    # average pop/acre
+        11:"AVGEMPDEN",    # average employment/acre
+        12:"AVGENTROPY",   # average entropy
+        13:"AVGENTROPYNW", # average non-work entropy
+        14:"AVGAREATYPE",  # average AREATYPE
+    }
+    TABLE_NAME_TO_NUMBER =  dict((v,k) for k,v in TABLE_NUMBER_TO_NAME.items())
+
+    # TABLE_NAMES          =  list(TABLE_NUMBER_TO_NAME[i] for i in range(1,len(TABLE_NUMBER_TO_NAME)+1))
+
+    ALL_VARS = list(TABLE_NUMBER_TO_NAME.values())
+
+    def __init__(self, file_dir, file_name = "walkSkim.h5"):
+        self.walk_skim_file = file_name
+        Skim.__init__(self, file_dir, [self.walk_skim_file])        
+
+
+    def getWalkSkimAttribute(self, orig_taz, dest_taz, attribute_name):
+        """
+        Returns the given walk skim attribute
+        """    
+        attribute_num = "%d" % WalkSkim.TABLE_NAME_TO_NUMBER[attribute_name]
+        return self.skim_table_files[self.walk_skim_file].root._f_get_child(attribute_num)[orig_taz-1][dest_taz-1] 
+
+
+    def getWalkSkimAttributes(self, orig_taz, dest_taz):
+        """
+        Returns all of the walk skim attributes in a dictionary (attribute name -> value)
+
+        If you want to access ``DISTANCE``::
+
+          walkSkimAttrs = walkSkim.getWalkSkimAttributes(otaz, dtaz)
+          dist = walkSkimAttrs["DISTANCE"]
+
+        """
+        walkAttributes = {}
+        for tableNum in list(WalkSkim.TABLE_NUMBER_TO_NAME.keys()):
+            walkAttributes[WalkSkim.TABLE_NUMBER_TO_NAME[tableNum]] = \
+                self.skim_table_files[self.walk_skim_file].root._f_get_child("%d" % tableNum)[orig_taz-1][dest_taz-1]
+
+        return walkAttributes
+
+    def getSkimTable(self, variable):
+        if variable.upper() not in WalkSkim.ALL_VARS: 
+            print("Requested Variable %s not available" % (variable))
+            exit(2)
+        table    = np.zeros((MAX_SF_ZONE, MAX_SF_ZONE))
+        tablenum = WalkSkim.TABLE_NAME_TO_NUMBER[variable.upper()]
+        table[:,:] = self.skim_table_files[self.walk_skim_file].root._f_get_child("%d" % tablenum).read()[:MAX_SF_ZONE,:MAX_SF_ZONE]
+        return table