Process of Revising SWIM Zone Structure

So You've Decided to Revise the Zones in SWIM

Hopefully you didn't come to this decision lightly. A lot of work and endlessly complicated steps are now in your future. You might be saying, "well, why don't we simplify the process, and make it easier". Good thought my friend, that's what I said too.

Under the Gen2 SWIM version, land use data (zonal data) was spread out over a seemingly endless number of differing input files with differing formats. The very idea of changing the zone structure was laughable and never attempted under Gen2. A lot of thought and energy went into properly capturing all the input files that were related to zone attributes and changing their storage to a single location in SWIM 2.5 (where we stand today). Therefore, a user can currently go to the SWIM version file (typically for the starting year, t0, but any network year will work), and review all ~160-170 zone attributes.

For people (like myself) with travel demand modeling experience, you might think, "Great! I have all the information I need, right in one place. What's the issue?" Unfortunately, as I found out too, SWIM is complicated. Like, really complicated. At the heart of the complication is SWIM's land use model AA. In a travel demand model, when you change the zone structure, you are really changing the land use inputs (givens) to a travel model. In the case of SWIM, you can edit the zones and re-tabulate the ~160 zone attributes, but that's just the beginning. You then need to ensure that all associated files are consistent with your overall changes and that the calibration of AA has not been impacted. That leads us to two levels of complicated that you can head down:

• Only impacting zones at the "alpha" level (transportation model) level
• Making zone changes that impact both the alpha and the beta (land use model) level zones

The steps are still more complicated than you might initially anticipate either way you choose, but there are more details to be aware of if impacting the beta level (land use model) zones in addition to just the alpha. Before you go any further with editing zones, make sure you are aware of some of the SWIM basics. I would start with the Coding Network Revisions Page and then build from there.

Steps to Update Zones When Isolating Changes to Just the Alpha Level

A good way to work towards improving SWIM while trying to avoid opening Pandora's box is to leave beta zones unaltered and just work on refining the zone structure within each Beta zone; editing the alpha zone level. Working at just the alpha zone level, and keeping the beta zone information as-is removes a fair amount of the complication but it doesn't fully simplify the process. Those who attempt will still need to work through:

• Saving existing zone data (for zones that aren't altered)
• Keeping zone totals unchanged at the beta zone level
• Rebuilding the alpha (and possibly beta) zone structure for all years with SWIM network (version) files
• Re-tabulate and load all zone attributes
• Updating the alpha connectors for all years with SWIM networks.
• Updating and iterating seed skims until stable
• Updating other "seed" zonal inputs that deal with information by zone

While this is a lot of work, there are several good reasons to consider these steps:

• Helping urban models nest properly in SWIM zones,
• Finding SWIM zones that are too big or too small to be represented properly,
• Finding a SWIM zone that simply cannot be fixed by connector placement.

The following sections will walk through the steps to modify SWIM zones at the alpha zone level.

Saving Existing Zone Data

Prior to any zone edits the original zone attributes need to be saved, example - “ZoneOrigAttributes.csv”. Additionally, the connector “net” file (example, “Orig_Connector.net” and "Orig_ZoneCoordinates.csv") and the zone centroid coordinates need to be saved out. This is done because a large percentage of the zones will likely not get altered in any zone revision, and saving out the original information will make it easier to add in revisions while retaining much of the network (zones) unchanged.

Here is the top three rows example for "Orig_ZoneCoordinates.csv":

$ZONE:NO	XCOORD	YCOORD
1	754548.8	1381558
2	755783.8	1381248
3	756201.4	1383425

Revising Alphas While Keeping Beta Totals Unchanged

Step one, take a breath in.. and out. Now that you have fresh oxygen, you can start hacking the existing alpha zone structure. For me, I chose to use existing shapes (like urban model zones) to split alpha zones. I used the "Split Polygons" function in ArcMap.

• First you turn the polygons you want to split with into polylines with Data Management Tools -> Features -> Feature to Line. Note that sometimes I had to hand extend the polylines I wanted to split with, to properly extend fully across the alpha zone I intended to split.
• The user then selects the series of polylines that they intend to represent the new boundary for the split alpha zone. and uses the "Split Polygons" button in ArcMap. This is an advanced function and requires a higher end license, for ODOT that means, Programs -> ArcGIS -> ArcGIS administrator -> Desktop -> choose “Advanced (ArcInfo) Concurrent Use. Remember to set the license back to "Basic (ArcView) Concurrent Use" when you are done using the advanced features.

For help in splitting alpha zones, I would typically have census or employment density maps, road networks, or geographic figures (like rivers) in the background to help choose the best lines to split on. If I felt a feature like a river was the best place to split. I would turn that river shape into a polyline (if it wasn't already a line) and select that line to split the alpha zone.

While completing the splits one also needs to renumber the new alpha zones that are created. To do this, I had a full list of all existing alpha zone numbers in use and would choose new alpha zone numbers accordingly (based on numbers that were available). I kept a table on the slide to reference which old alpha zone numbers were split into which new alpha zone numbers. This was very helpful in re-calculating zonal attributes and helped keep a good accounting as the process got more complicated, which it does immediately in the following steps. Additionally, a master list should also be kept, here is an example that was titled "TextRecrodOfZoneChanges.txt" and will appear in script snipets below (this file was used to identify which zones need new information; new zone attributes, new connectors, new zone centroids):

# Alphas to remove
Rm <- c(#4001,4004,4010,4011,4026,4071,4072,4073,4093,4094,4096,4140,4141, # halo cut
        3068,3076,3069,3058,3059,3063,3066,3065,3070,3906,3910,3902,2913,2159, # combined list
        2917,2902,2904,2939,2935,2916,2936,2941,2937,2947,2943,2948,2949,2183,
        3909,2923,2919,2933,2928,2924,2931,2932,2801,2813,2814,2808,2803,2575,
        2809,2081,2079,3276,2083,2651,2649,2629,2623,2626,2606,2635,2636,303,
        2625,2612,2152,2142,2128,2127,2133,2135,2013,1387,2673,2880,2171,
        3848,3819,3816,3898,2844,3847,3241,3193,3291,3273,2864,2645, # Beta merges  
        2818,3055,2338,1470,2670,2861,2846,2779, # Remove the hand moved betas

        # From readme  (shape edited in some way, so data needs to be re-calc-ed)
        3054, #duplicates 2846,1470,2338,2670,2818,2844,3193,3241,3273,3291,3819,3816,3847,3848,3898
        2971,2772,2800,2802,2807,3071,3907,2856,2848,2423,2121, #duplicates 2846,2848,2645,2861,2864
        1319,2190,2297,2674,2273,2972,2867,2598,2759,2756,2781,2765,2758, #duplicates 2772,2781,2780
        2774,2766,1724,1729,2777,2778,2275,2776,2768,2767,2008,1731,2761,1826,1730, #duplicates 2774
        1824,1827,1725,1726,1727,2747,2748,2751,2752,2753,2750,2749,1716,1708, #duplicates 1724,1729
        2717,2718,2713,2757,2006,2007,2721,2727,2725,2711,2729,2730,2715,2716,2760,2676,
        2678,2684,2695,2712,2719,1704,2672,2734,2733,2697,2969,2683,2688, #duplicates 1708,2695,2676
        2700,2701,2702,2699,2706,2708,2732,2735,2709,2710,2707,2704,2703,2690,2691,2692,
        2705,2679,2680,2681,2682,2862,2860,3007,3011,3010,3009,3003,3008,2986,2989,3006,
        3004,3001,2968,3000,2998,2999,2995,2996,2997,2994,2993,2992,2990,2983,2985, #duplicates 2969
        2982,2981,2984,2987,2976,2977,2978,2980,2979,2966,2964,2965,2875,2874, #duplicates 2867,2972
        2866,2870,2314,2312,2313,2311,2080,2075,2613,2622,2621,
        2780, # one more special case (see Rn above), where the zone was re-numbered, but then edited so extremely, that it needed to be re-worked entirely
        3040, # extra zone from the Alpha zone to New Beta and 
        535,521,553,537,538,1290,1291,1292,281,282,283,285,273,2256,2299, #SplitZones
        1559,1588,2896,2838,2830,2877,2881,3860,2878,2696,2893,2973,2974, #SplitZones
        3062,3075,3060,3064,3057,2907,2899,2901,2910,2918,2905,2911,2915,2945,2934, # the following rows identified combined zones
        2946,2942,2950,2953,2921,2920,2926,2925,2929,2930,2820,2806,2804,2082,2077,
        2084,2630,2627,2632,2624,2605,2607,2610,2150,2134,2126,2132,2136,1392,2879,
        2172,2158,2182,2576,302,
        2285:2290,2292:2296, # woodburn zones 2297 already covered in beta zone edits above
        3840,3056,3067,3277) # additional zones impacted in Pendelton edits

Spiting Alpha zones is the easiest edit that can be made. The total information is known for all zonal attributes within that alpha, the only thing left to do is determine the percentage of each attribute that goes to each new zone. The problem starts to become increasingly complicated if split zones need to be merged into other splits of other alpha zones, which is a pretty likely case given the size and existing shaping of the SWIM zones.

The process is similar as is laid out above, however now two alphas (or more) are getting split and then pieces of each are getting merged together. A separate table was kept to keep a record of how each alpha zone sections were being combined versus split. As opposed to trying to keep track of what percentage needed to come from which zone to another zone, these zones were marked as needing a total re-calculation of attributes, versus a simple split. The full totals for the total model area for each attribute were still held true, but all manipulated zones (beyond just getting split) were tagged with a full re-allocation, to be distributed within the model area (untouched zones retained all their attributes).

Besides tracking merged zones in a different way; additional zone attributes also need to be considered when merging sections of alpha zones:

• Beta
• PUMA
• County
• ALD regions

Basically, one needs to verify that the process of merging alpha zones isn't going to accidentally create new regions, districts, or aggregations working at higher levels of the model... unless you really mean to.

Rebuilding the Alpha (and possibly beta) Zone Structure for All Years with SWIM Network (version) Files

After all the zone edits have been made and finalized in ArcGIS (or whatever software is being used to manage the zone boundary revisions), the changes will need to be brought back into Visum (the network and input management software for SWIM). It's typically not just as simple as deleting the previous zones (alpha) and main zones (beta) layers in Visum and importing the revised shapefiles; although those are the basic steps. Before that delete and import can occur (for both the alpha and beta zone layers), the alpha layer needs to be checked, the beta layer needs to be re-built from the revised alpha zone layer, and the resulting beta zone layer needs to be checked. The example script below does this layer creation and basic checking (cleaning and checking is also completed in ArcGIS manually). Here are some of the steps that the provided script below completes:

• Ensuring that all zones have numbers
• Ensuring that no zone numbers are repeated
• Verifying that the total number of zones is the anticipated amount
• Ensuring that none of the alpha (and then the creation of the beta layer) have mutli-part elements in the polygons - no overlaps or holes.
• Pulls over a handful of geographic fields like County, PUMA, ALDRegion...
• The script also provides examples of edits that can be documented and made in R as opposed to by hand in ArcGIS.

After these steps are completed either manually or with something like the script provided below, then the old zones and main zones can be deleted in Visum and the new alpha and beta shapefiles can be imported as zones and main zones into the Visum version files for all the reference networks (however many different networks exist for the different time periods; currently just 2010 and then a 2016+, but soon there will at least be a third future network).

Important reminder - the external station series is the 5000 series in the alpha zone layer. Unless the intention is to alter the external stations, when zones are being deleted from the given Version file, filter the zones so that the 5000 series (zones greater than or equal to 5000) are kept (not deleted). The development and import process laid out in these steps is only focused on the internal zones (4999 and less).

Also note that the script pulls in a SWIM shapefile and an attribute csv that are not included here - however the path to where they sit on ODOT's "6420Only" drive is provided for those within ODOT. An example of the "SplitZonesRecord.csv" that is called into the script is shown down in: Updating and iterating seed skims until stable. However, note that the example provided just shows the format. The actual file used in the example script has many more zones identified.

# 3-29-17
# This script is a record for how the zone revisions (in the alpha zone layer) were used to create a new Beta zone layer.

library(rgeos)
library(rgdal)

a <- readOGR(dsn="//s6000e/6420only/Statewide/Tools/SWIM/SWIM_2_5/ModelBuild/CalibrationWork/2017_reducedHalo_NetworkAndSkims/Revised_Alpha", layer="NewAlpha")

# Update fields taken from GIS
a@data <- a@data[,c("OBJECTID","NO","Shape_Area")]
names(a@data)[3] <- "AREASQFT"

# verify that no zone numbers are repeated, max should be = 1
max(table(a@data$NO))
# verify no NA's, sum should be zero
sum(is.na(as.numeric(a@data$NO)))

# verify no multi point polgons exist
table(unlist(lapply(a@polygons,function(x) length(x@Polygons)))) # all zones should have one part

# read in pre-existing data and linkages to new zones
atts <- read.csv("ZoneOrigAttributes.csv", as.is=T)
# get rid of combined zones:
atts <- atts[!atts$NO %in% c(3068,3076,3069,3058,3059,3063,3066,3065,3070,3906,3910,3902,2913,2917,2902,2904,2939,2935,2916,2936,2941,2937,2947,2943,2948,2949,3909,2923,2919,2933,2928,2924,2931,2932,2801, 
                             2813,2814,2808,2803,2809,2081,2079,3276,2083,2651,2649,2629,2623,2626,2606,2635,2636,2625,2612,2152,2142,2128,2127,2133,2135,2013,1387,2673,2880,2171,2159,2183,2575,303),]  
# beta zone merge list
atts <- atts[!atts$NO %in% c(3848,3819,3816,3898,2844,3847,3241,3193,3291,3273,2864,2645),]  

# remove numbers for zones hand moved
atts <- atts[!atts$NO %in% c(2818,3055,2338,1470,2670,2861,2846,2779),]
for(i in 1:8) atts[atts$NO==c(3818,3855,3838,3870,3873,3861,3846,2780)[i], "NO"] = c(2818,3055,2338,1470,2670,2861,2846,2779)[i]

# shift "NO" to rownames
rownames(atts) <- atts$NO

# need to hand switch a handful of Beta zone tags so that we don't create multi-part beta zone shapes
atts["2138","BZONE"] <- 2139 # was 2134
atts["2122","BZONE"] <- 2123 # was 2122
atts[atts$BZONE==2122,"BZONE"] <- 2145 # shifting to the centered zone on the correct side of the bay
atts["1043","BZONE"] <- 1000 # was 1025
# I don't like this solution for Beta 2009 (Astoria-Warrenton), but it's the best I have right now
atts[c("2009","2010"),"BZONE"] <- 2011 # was 2009
atts[atts$BZONE==2009,"BZONE"] <- 2019 # shifting to main Warrenton Zone

# Woodburn hand shift
atts["2285","BZONE"] <- 2286 # was 2283
atts["2289","BZONE"] <- 2287 # was 2288

# check Beta count, should be 517 (518 - 13 = 505, but plus 12 externals = 517
length(unique(atts$BZONE))  #517

# linkage from old zones to the new splits
Splits <- read.csv("SplitZonesRecord.csv", as.is=T)
# Only capture new zone numbers
Splits <- Splits[Splits[,1] != Splits[,2],1:2]
# verify that none of the new zone numbers exist in the previous zone system
Splits[as.character(Splits[,2]) %in% rownames(atts),]
# duplicate all split rows
sum(!as.character(Splits[,1]) %in% rownames(atts)) # needs to be zero, meaning that all zones are avaliable in the listing
temp <- atts[as.character(Splits[,1]),]
rownames(temp) <- temp$NO <- Splits[,2]
length(unique(rownames(temp)))==nrow(Splits) # needs to be true 
# combine the split list with the atts list
atts <- rbind(atts, temp)
rm(temp,i, Splits)

# find any zones that don't belong and hunt them down with murderess vengence.
atts$NO[!atts$NO %in% a@data$NO]   # 5000 series (external zones) can stay)
a@data$NO[!a@data$NO %in% atts$NO ]

# merge additional zone attributes to the 
NewAtts <- cbind(a@data[,2:3],atts[as.character(a@data$NO),names(atts)[names(atts)!= "NO"]]) 
write.csv(NewAtts[order(NewAtts[,1]),], "NewAlphaAttributes.csv", row.names=F)

# Add Beta zone data to shape, to aggregate up
a@data$BZONE <- atts[as.character(a@data$NO),"BZONE"]

# write out updated Alpha
writeOGR(a, "//s6000e/6420only/Statewide/Tools/SWIM/SWIM_2_5/ModelBuild/CalibrationWork/2017_reducedHalo_NetworkAndSkims/Revised_Alpha", "Alpha", driver="ESRI Shapefile")

# create a new Beta zone shape
b <- gUnaryUnion(a, a@data$BZONE)
DATA <- as.data.frame(unlist(lapply(b@polygons,function(x)x@ID)))
rownames(DATA) <- unlist(lapply(b@polygons,function(x)x@ID))
names(DATA) <- "NO"
b <- SpatialPolygonsDataFrame(b, DATA)
 plot(b)
 length(sort(table(b@data$NO)) )
table(unlist(lapply(b@polygons,function(x) length(x@Polygons))))  # all zones should have one part
b@data[which(unlist(lapply(b@polygons,function(x) length(x@Polygons)))>1),]
writeOGR(b, "//s6000e/6420only/Statewide/Tools/SWIM/SWIM_2_5/ModelBuild/CalibrationWork/2017_reducedHalo_NetworkAndSkims/Revised_Alpha", "NewBeta", driver="ESRI Shapefile")

# did manual cleaning, now checking the beta layer

d <- readOGR(dsn="//s6000e/6420only/Statewide/Tools/SWIM/SWIM_2_5/ModelBuild/CalibrationWork/2017_reducedHalo_NetworkAndSkims/Revised_Alpha", layer="NewBeta")
# verify that no zone numbers are repeated, max should be = 1
max(table(d@data$NO))
# verify no NA's, sum should be zero
sum(is.na(as.numeric(d@data$NO)))

# verify no multi point polgons exist
table(unlist(lapply(d@polygons,function(x) length(x@Polygons))))  # all zones should have one part

Re-Tabulate and Load All Zone Attributes

At this point, a number of zone boundaries have been re-drawn, the shapes and zone numbers have been checked to verify that they are correct. Old zones structures have been removed from the Visum version files and the revised zone boundaries have been loaded back in. Now, the version files should have clean zone boundaries, but largely empty zone attribution data.

When dealing with any number of zone edits beyond a small handful (like less than 5 zones altered), it is highly recommended (basically required) that scripting be brought in to properly keep track of all the 160+ zone attributes. An example script is provided below. For those within ODOT, the script ("AttributeUpdate.R") and related files can be found on the 6420Only at: \s6000e\6420only\Statewide\Tools\SWIM\SWIM_2_5\ModelBuild\CalibrationWork\2017_reducedHalo_NetworkAndSkims

Overall the script deals with several major components:

• Ensures that the number of base (2010) households remains the same and works to reallocate households to revised zones using "PopSyn_2_SWIM.R" (18 fields)
• Ensures that base dollars by employment categories remains the same and works to reallocate employment dollars in revised zones using a statewide 2010 employment layer (52 Fields).
• For altered zones, available land by 55 zoning categories is re tabulated from "ZoningByAlphaZone.txt" file (55 fields)
• Ensures that the base level of sqft by floorspace type remains the same and reallocates floorspace by household and employment allocation (18 fields)

Before getting into the script it's important to note that in the bullets above a series of significant prep steps are completed or assumed complete.

1. The process of reallocating households relies on a census block layer existing for the state, and that a synthetic population has been generated for the State of Oregon (assuming no zone changes outside of Oregon). The synthetic population result is relied on to tabulate the joint distribution of base year households by revised alpha zone for the 18 income and household size categories.
2. The process of re-tabulating build able land (zoned land) requires that the user has split the latest statewide zoning layer with the revised zone structure. The process to do this is laid out in "Steps_2build_SWIM_ZoningLayer.docx" but can be summarized as:
   • Union the Statewide Zoning Layer with the Statewide Water Layer
   • Union the Resulting Database with the Revised SWIM alpha zones
   • Output the resulting union for processing in R (script below)
3. The process for re-tabulating employment dollars assumes that you have an employment point file for the state (or at least the region where the zones are being revised in SWIM).

# Alex Bettinardi
# 5-22-17

# A script to create the new Alpha Zone level attributes for the Visum file (for internal zones)
 library(maptools)
 library(mgcv)

# Read in original attribution information:
atts <- read.csv("ZoneOrigAttributes.csv", as.is=T) # original 2998 zones minus 13 removed Halo zones (total = 2985)
rownames(atts) <- atts$NO

# read in completed zone attributes under "Revised_Alpha/alpha_to_beta.R"
NewAtts <- read.csv("Revised_Alpha/NewAlphaAttributes.csv", as.is=T)
# quickly add GRIDACRES to this table off of the SQFT calc'd
NewAtts$GRIDACRES <- NewAtts$AREASQFT / 43560
rownames(NewAtts) <- NewAtts$NO
#sum(NewAtts$AREASQFT)/sum(atts$AREASQFT )

# Identify the Fields to update (to append to the fields already calculated
FullUpdateList <- names(atts)[!names(atts) %in% names(NewAtts)] # 143
# Break the update list into four categories, HH, FLR, Zoning, and Emp
HHnames <- FullUpdateList[grep("^HH", FullUpdateList)] #18 
FLRnames <- FullUpdateList[grep("^FLR", FullUpdateList)]  #18
ZoningNames <- sort(c(FullUpdateList[nchar(FullUpdateList)<6],"EFU160"))  #55
EmpNames <- FullUpdateList[!FullUpdateList %in% c(HHnames,FLRnames,ZoningNames)]  #52
#sum(atts[,ZoningNames])/sum(atts$AREASQFT)
#sum(atts[atts$STATE=="OR",HHnames])/1518938

# zoning Layer incomplete
zoningIncomplete <- T

  if(zoningIncomplete){  # only do this step if there are adjusted zones that are missing zoning data.
     # new zones                                                     # these were all perfect, just missing roads... gives confidence, that we reached a good threashold
                                                                     #,1827,2690,2679,2681,1824,2688)                                                                2945,2946,2950,2942
     noZoningZones        <- c(2079,2820,2256,2649,2075,2629,2084,2182,2134,3064,2607,2623,2136,2606,3059,2159,2622,2133,2256,2081,2807,2077,3055,2801,3065,3068,3071,2152,3062,3267,3277,3277,3075,3067)
     #switching order - now the bottom (names) are the new zones, the upper (vector) are the old zones - this allows the old zone to get more than one zone
     names(noZoningZones) <- c(2080,2820,2256,2627,2075,2632,2084,2126,2134,3064,2607,2624,2136,2605,3060,2158,2622,2132,2246,2082,2807,2077,3054,2802,3057,3062,3071,2150,3056,3276,3277,3076,3075,3067)
    
     # update remove zones with plenty of coverage in the latest file    
     noZoningZones <- noZoningZones[!names(noZoningZones) %in% as.character(c(2136,2134,2075,2080,2150,2084,2082,2256,2126,2077,2246,2132))]   
  
     # pull out the zoning detail for the empty-ish zones
     zoningGaps.. <- atts[as.character(noZoningZones),ZoningNames]
     rownames(zoningGaps..) <- names(noZoningZones)
     rm(noZoningZones)
  }

# Calc and save total col sums for all fields, so we can get roughtly the same total
TotalsSave <- colSums(atts[,FullUpdateList])

# Assemble list of all zones to remove 
# Alphas re-numbered
Rn <-        c(2818,3055,2338,1470,2670,2861,2846)#,2779) Special case, zone 2780 handled below
names(Rn) <- c(3818,3855,3838,3870,3873,3861,3846)#,2780)
# switch the Alpha numbering
atts[names(Rn),"NO"] <- Rn
atts <- atts[!rownames(atts) %in% as.character(Rn),]
rownames(atts) <- atts$NO

# Alphas to remove
source("TextRecordOfZoneChanges.txt")  # creates the object "Rm" (for zones to remove)      
        
# remove all zones that have been shifted in some way
as.character(Rm)[!as.character(Rm) %in% rownames(atts)]
atts <- atts[!rownames(atts) %in% as.character(Rm),]
rm(Rm, Rn)

# separate externals
ext <- atts[as.character(5001:5012),]
atts <- atts[!rownames(atts) %in% rownames(ext),]

# verify all "keep" zones exist in the new structure
rownames(atts)[!rownames(atts) %in% rownames(NewAtts)]
# Develop the list of zones to update
update.Zn <- rownames(NewAtts)[!rownames(NewAtts) %in% rownames(atts)] # length(update.Zn) is 339 which is good - correct

# Complete New zone data "container"
for(Name in FullUpdateList){
   NewAtts$temp <- 0
   names(NewAtts)[length(names(NewAtts))] <- Name
}
rm(Name)
NewAtts <- NewAtts[,names(atts)]
# populate the "keep" zones and append the external zones
NewAtts[rownames(atts),FullUpdateList] <- atts[,FullUpdateList]
NewAtts <- rbind(NewAtts,ext[,names(NewAtts)])

# a couple of checks
sum(NewAtts[rownames(atts),"AREASQFT"])/sum(atts$AREASQFT )
sum(NewAtts[,HHnames])/sum(atts[,HHnames])
sum(NewAtts[,EmpNames])/sum(atts[,EmpNames])
sum(NewAtts[,ZoningNames])/sum(atts[,ZoningNames])
sum(NewAtts[,FLRnames])/sum(atts[,FLRnames])

rm(atts, ext)

# for each "update" zone:
# Read shapes to be used in the for loop
# 1. SWIM zone boundaries
  SWIM_ <- readShapePoly("//s6000e/6420only/Statewide/Tools/SWIM/SWIM_2_5/ModelBuild/CalibrationWork/2017_reducedHalo_NetworkAndSkims/Revised_Alpha/NewAlpha.shp")
# 2. Statewide population 
  pop_ <- readShapePoints("//s6000e/6420only/Data/Census/Census2010/Censusshapes/OregonBlocks/stateblock_points.shp")
  # read statewide pop syn..
  source("PopSyn_2_SWIM.R")  
  # create an empty vector of all hh categories to fill by zone
  HHVec <- rep(0, length(HHnames))
  names(HHVec) <- HHnames 

# 3. Something with zoning...
# read my instructions -  "Steps_2build_SWIM_ZoningLayer.docx"  
  #zoning <- read.csv("MidtermShapes/ZoningByAlphaZone.txt", as.is=T)
  # updated this section 5-4-17 to read the latest zoning data
  # zoning <- read.csv("MidtermShapes/ZoningByAlphaZone_2ndUpdate.txt", as.is=T)
  # updated this section 9-22-17 to align with Pendleton edits
  zoning <- read.csv("MidtermShapes/ZoningByAlphaZone_3rdUpdate.txt", as.is=T)
  sqft.AzZc <- tapply(as.numeric(gsub(",","",zoning$Shape_Area)), list(as.numeric(gsub(",","",zoning$NO)),zoning$orZCode),sum)
  sqft.AzZc[is.na(sqft.AzZc)] <- 0
  sqft.AzZc <- sqft.AzZc[rownames(sqft.AzZc) != "0",colnames(sqft.AzZc) != ""]
  sqft.AzZc <- sqft.AzZc[,colnames(sqft.AzZc) %in% ZoningNames]
  #colnames(sqft.AzZc)[!colnames(sqft.AzZc) %in% ZoningNames]
  #ZoningNames[!ZoningNames %in% colnames(sqft.AzZc)]
  rm(zoning)
  
  # check to see which zones that need zoning information don't have any, or don't have enough:
  #temp<-sort(round(100*rowSums(sqft.AzZc[update.Zn,])/NewAtts[update.Zn,"AREASQFT"]))
  #names(temp)[temp<80][!names(temp)[temp<80] %in% rownames(zoningGaps..)]
  #sort(round(rowSums(sqft.AzZc)[rownames(zoningGaps..)]/NewAtts[rownames(zoningGaps..),"AREASQFT"],2))
  #rm(temp)
  # zones with less than 80% coverage:
  # insert the percenatages from the surrogate zones, here...
  
  if(zoningIncomplete){  # only do this step if there are adjusted zones that are missing zoning data.
     
     # now add the old zone data to the zones that have zoning data gaps
     sqft.AzZc[rownames(zoningGaps..),] <- sqft.AzZc[rownames(zoningGaps..),] + as.matrix(zoningGaps..[,colnames(sqft.AzZc)])
     rm(zoningGaps..)
  }
  # after investigation (zone 3877) it was noted that the zoning layer has some duplicate polygons
  # so in the small handful of cases where the sqft totals to greater than 100% 
  # factor all sq-footage back to total and trusted AREASQFT, from the original zones
  sqft.AzZc[update.Zn,] <- sqft.AzZc[update.Zn,]*(NewAtts[update.Zn,"AREASQFT"]/rowSums(sqft.AzZc[update.Zn,]))

  # update zoning fields
  NewAtts[update.Zn,ZoningNames] <- sqft.AzZc[update.Zn,ZoningNames]
  #sum(NewAtts[,ZoningNames])/sum(TotalsSave[ZoningNames])
  #sum(NewAtts[,ZoningNames])/sum(NewAtts$AREASQFT)
  
  #rm(sqft.AzZc)

# 4. Statewide employment
  emp_ <- readShapePoints("//s6000e/6420only/Statewide/Tools/SWIM/SWIM_2_5/ModelBuild/CalibrationWork/2017_reducedHalo_NetworkAndSkims/MidtermShapes/EmpReference.shp")

# Emp crosswalk processing
# read in crosswalk, expand and look for any numbers that aren't 6 digits, overlap and missing
  empCW <- read.csv("cw_NAICS.csv", as.is=T)
  table(c(nchar(empCW$MinNAICS),nchar(empCW$MaxNAICS)) )  # all equal 6
  naicsCW <- sapply(1:nrow(empCW), function(i) seq(empCW[i,1],empCW[i,2],1))
  length(unlist(naicsCW)) ==  length(unique(unlist(naicsCW)))
  naicsCW. <- unlist(naicsCW)
  naicsCW <- unlist(sapply(1:nrow(empCW), function(i) rep(empCW[i,3],length(naicsCW[[i]]))))
  names(naicsCW) <- naicsCW.
  rm(empCW, naicsCW.)
  sort(unique(emp_$NAICS)[!unique(emp_$NAICS) %in% names(naicsCW)] )  # all NAICS accounted for 
  
  # create an empty vector of all employment categories to fill by zone
  EmpVec <- rep(0, length(EmpNames))
  names(EmpVec) <- EmpNames 
  
# non-res floorspace cw
  nresFLRcw <- read.csv("cw_EMP2FLR.csv", as.is=T)
  nresFLRcw[,3] <- toupper(gsub(" ", "_",nresFLRcw[,3]))
  nresFLRcw[,3][!nresFLRcw[,3] %in% FLRnames] 
  nresFLRcw[,1][!nresFLRcw[,1] %in% EmpNames]
  #FLRnames[!FLRnames %in% nresFLRcw[,3]]   

# res floorspace cw
  far <- as.matrix(read.csv("far.csv",as.is=T, row.names=1)[1:6,] )
  rownames(far) <- gsub(" ", "_", rownames(far))

for(Zn in as.numeric(update.Zn)){   #rownames(NewAtts)[NewAtts$STATE=="OR" & NewAtts$NO<4000]){#
# Recalc HH based on block points and syn pop
     blocks <- as.character(pop_$GEOID10[in.out(SWIM_@polygons[[match(Zn,SWIM_$NO)]]@Polygons[[1]]@coords, pop_@coords)])
     hh.Zn <- table(toupper(hh[hh$GEOID10 %in% blocks,"ind"]))
     # zero out HH Vector
     HHVec[] = 0
     HHVec[names(hh.Zn)] <- hh.Zn
     NewAtts[as.character(Zn),HHnames] <- HHVec[HHnames]
    
# Recalc Emp based on emp points and NACIS crosswalk  (see "MidtermShapes/EmpReference")
     Eind <- in.out(SWIM_@polygons[[match(Zn,SWIM_$NO)]]@Polygons[[1]]@coords, emp_@coords)
     emp.Zn <- tapply(emp_$AVGEMP[Eind],naicsCW[as.character(emp_$NAICS[Eind])],sum)
     #names(emp.Zn)[!names(emp.Zn) %in% EmpNames]
     if(is.na(emp.Zn["CNST"])){ 
        emp.Zn <- c(0,emp.Zn)
        names(emp.Zn)[1] <- "CNST"
     }
     # zero out Employment Vector
     EmpVec[] = 0
     # fill Employment Vector
     EmpVec[names(emp.Zn)[names(emp.Zn) %in% EmpNames]] <- emp.Zn[names(emp.Zn)[names(emp.Zn) %in% EmpNames]]
     
     #special cases fills
     # note that employment is duplicated, but we are only using employment to allocate emp dollars by category, 
     # so these are just relatative numbers by zone within a given employment category, so there actual value doesn't matter
     # just their share in relation to the other zones being editied.
     # CNST
     EmpVec[c("CNST_MAIN_XXX","CNST_NRES_XXX","CNST_OFFC_OFF","CNST_RES_XXX")] <- emp.Zn["CNST"]
     # ENGY_OFFC_OFF just gets  ENG emp in the area
     EmpVec[c("ENGY_OFFC_OFF")] <- sum(EmpVec[grep("ENGY",names(EmpVec))])
     # GOV_OFFC_OFF just gets GOV_ADMN_GOV
     EmpVec[c("GOV_OFFC_OFF")] <- EmpVec["GOV_ADMN_GOV"]
     # INFO_INFO_OFF_LI gets INFO_INFO_OFF
     EmpVec[c("INFO_INFO_OFF_LI")] <- EmpVec["INFO_INFO_OFF"]
     # K12_K12_OFF gets K12_K12_K12
     EmpVec[c("K12_K12_OFF")] <- EmpVec["K12_K12_K12"]
     # MFG_HTEC_LI gets MFG_HTEC_HI
     EmpVec[c("MFG_HTEC_LI","MFG_HTEC_HI")] <- EmpVec["MFG_HTEC_HI"]/2
     # MFG_FOOD_LI gets MFG_FOOD_HI
     EmpVec[c("MFG_FOOD_LI","MFG_FOOD_HI")] <- EmpVec["MFG_FOOD_HI"]/2
     # MFG_HVTW_LI gets MFG_HVTW_HI
     EmpVec[c("MFG_HVTW_LI","MFG_HVTW_HI")] <- EmpVec["MFG_HVTW_HI"]/2
     # MFG_OFFC_OFF just the total MFG emp in the area
     EmpVec[c("MFG_OFFC_OFF")] <- sum(EmpVec[grep("MFG",names(EmpVec))])
     # RES_OFFC_OFF gets  RES_FORST_LOG
     EmpVec[c("RES_OFFC_OFF")] <- EmpVec["RES_FORST_LOG"]
     # RET_STOR_OFF and RET_NSTOR_OFF get RET_STOR_RET
     EmpVec[c("RET_STOR_OFF","RET_NSTOR_OFF")] <- EmpVec["RET_STOR_RET"]
     # TRNS_TRNS_OFF gets TRNS_TRNS_WARE
     EmpVec[c("TRNS_TRNS_OFF")] <- EmpVec["TRNS_TRNS_WARE"]
     # UTL_OTHR_OFF gets UTL_OTHR_OFF_LI
     EmpVec[c("UTL_OTHR_OFF")] <- EmpVec["UTL_OTHR_OFF_LI"]
     # WHSL_OFFC_OFF gets WHSL_WHSL_WARE
     EmpVec[c("WHSL_OFFC_OFF")] <- EmpVec["WHSL_WHSL_WARE"]
     
     # fill employment values into the full table - after the for loop these will be factored up to dollars as opposed to employees 
     NewAtts[as.character(Zn),EmpNames] <- EmpVec[EmpNames]

# think about away to fumble together FLR starting data. (in sqft)  
# Will likely need to be from previous zone total and then aligning with coverage data 
# in some intellgent way - good luck

# answer see employment-> no-res and zoning-> res crosswalks above this for loop
# use those to assign FLR like-ish activity, then after the for loop use "activity" to propotion displaced FLR sqft
     for(flr in unique(nresFLRcw$Floorspace.options)) NewAtts[as.character(Zn),flr] <- sum(NewAtts[as.character(Zn),nresFLRcw[nresFLRcw[,3] == flr,1]])
     NewAtts[as.character(Zn),rownames(far)] <- rowSums(sweep(far,2,as.numeric(NewAtts[as.character(Zn),colnames(far)]),"*"))

} # end for loop

# clean up
rm(hh.Zn, hh, blocks, naicsCW, emp.Zn, Eind, HHVec, EmpVec, Zn, flr, zoningIncomplete)

# factor employment to dollars
emp.Dol <- TotalsSave[EmpNames] - colSums(NewAtts[!rownames(NewAtts) %in% update.Zn,EmpNames])
NewAtts[update.Zn,EmpNames] <- sweep(NewAtts[update.Zn,EmpNames],2, emp.Dol/colSums(NewAtts[update.Zn,EmpNames]),"*")
summary(colSums(NewAtts[,EmpNames])- TotalsSave[EmpNames])

rm(emp.Dol, EmpNames)

# factor update built FLR for udpated zones
flr.Sq <- TotalsSave[FLRnames] - colSums(NewAtts[!rownames(NewAtts) %in% update.Zn,FLRnames])
NewAtts[update.Zn,FLRnames] <- sweep(NewAtts[update.Zn,FLRnames],2, flr.Sq/colSums(NewAtts[update.Zn,FLRnames]),"*")
summary(colSums(NewAtts[,FLRnames])- TotalsSave[FLRnames])

rm(flr.Sq, FLRnames,nresFLRcw, far)

summary(colSums(NewAtts[,HHnames])-TotalsSave[HHnames] )
colSums(NewAtts[,HHnames])-TotalsSave[HHnames]
sum(NewAtts[update.Zn,HHnames])   
sum(colSums(NewAtts[,HHnames])-TotalsSave[HHnames])/ sum(NewAtts[update.Zn,HHnames])  

# final check on any zones with wackey space
temp <- cbind(as.numeric(rownames(NewAtts)[abs(rowSums(NewAtts[,ZoningNames])-NewAtts$AREASQFT)>1000000]),abs(rowSums(NewAtts[,ZoningNames])-NewAtts$AREASQFT)[abs(rowSums(NewAtts[,ZoningNames])-NewAtts$AREASQFT)>1000000],(rowSums(NewAtts[,ZoningNames])/NewAtts$AREASQFT)[abs(rowSums(NewAtts[,ZoningNames])-NewAtts$AREASQFT)>1000000])
#temp[order(temp[,3]),]
#temp[,1][as.character(temp[,1]) %in% update.Zn]
temp <- cbind(NewAtts[(100*rowSums(NewAtts[,ZoningNames])/NewAtts$AREASQFT) < 90, c("NO","MPO","AREASQFT")] , ZoningSQFT=rowSums(NewAtts[,ZoningNames])[(100*rowSums(NewAtts[,ZoningNames])/NewAtts$AREASQFT) < 90 ])
temp <- temp[!is.na(temp$NO),]
temp$Per <- round(100*temp[,4]/temp[,3]) 
temp <- temp[order(temp$Per),]
temp$newTot <- rowSums(sqft.AzZc[rownames(temp),])
temp$NewPer <-  round(100*temp$newTot/temp[,3])
#temp[order(abs(temp$NewPer-temp$Per)),]
fill.Zn <- as.character(sort(temp[abs(temp$NewPer-temp$Per)>8,"NO"]))

sqft.AzZc[fill.Zn,] <- sqft.AzZc[fill.Zn,]*(NewAtts[fill.Zn,"AREASQFT"]/rowSums(sqft.AzZc[fill.Zn,]))
NewAtts[fill.Zn,ZoningNames] <- sqft.AzZc[fill.Zn,ZoningNames]

#when the updated zoning layer comes from DLCD - we will need to udpate all Oregon zones that have 
# ~60-70% complete coverage under the new layer (found some pretty bad wholes (see jeffereson and brookings) in the current zoning data.

# smashed all back together, export and the attribution should be done, 
write.csv(NewAtts,"RevisedAlphaAttributes.csv",row.names=F)

# then move on to connector / network update for the new zones 

Once all this is complete, the process creates "RevisedAlphaAttriutes.csv" The user can alter the header to add an initial row to the file and alter the first two cells to:

Cell A1 = $VISION
Cell A2 = $ZONE:NO

These changes will allow the table to be easily copied and pasted into the revised SWIM Version files (networks). In Visum this is completed by copying the altered attribute table to the clipboard, opening the Visum version file of interest, listing the zones table in Visum and then pressing the paste button and following the dialogue.

After the copy and paste is complete, the version file can be saved and updating the zone attribution is complete.

Updating the Alpha Connectors for All Years with SWIM Networks.

So far the edited zone boundaries and and revised (re-tabulated) zone attributes have been added back into the Visum version files. Now the new zones have to be re-connected to the highway network. At this point all of the zones that were originally deleted from Visum (if following these instructions, that's all zones except for the 5000 series) are no longer connected to the network. As part of the first step, the instructions had the user save out the original connectors and zone centriods. The following script quickly process those files, removing the information for edited zones, so that all un-edited zones can be quickly reconnected:

#Alex Bettinardi
# 4-27-17

# A script to help process the connectors to be used... (kept)
# new connectors will be built by hand.

# Alphas to remove (to create new connectors for:
source("TextRecordOfZoneChanges.txt")  # creates the object "Rm" (for zones to remove)

# look for any additional zone names that I might be missing in my list
New.Zn <- as.numeric(read.csv("Revised_Alpha/NewAlphaAttributes.csv", as.is=T)[,1]) 
Orig.Zn <- as.numeric(read.csv("ZoneOrigAttributes.csv", as.is=T)[,1]) 
NewNumbs <- New.Zn[!New.Zn %in% Orig.Zn]

# Add beta zone number shifts to the Rm list
Rm <-   c(Rm,2818,3055,2338,1470,2670,2861,2846,3818,3855,3838,3870,3873,3861,3846)

# there area new zone numbers not in the Rm record that need to have new connectors
Rm <- unique(c(Rm,NewNumbs[!NewNumbs %in% Rm]))
rm(NewNumbs, New.Zn, Orig.Zn)

#The external stations will not be removed, so those connectors need to be removed from the import list - since they will still be in Visum

# read in original connector net file to be edited.  
conn <- readLines("Orig_Connector.net")
conn.Zn <- as.numeric(unlist(lapply(strsplit(substring(conn,1,5),";"), function(x) x[1]))[14:length(conn)])
writeLines(conn[c(1:13,13+(1:length(conn.Zn))[!conn.Zn %in% Rm])],"Orig_Connectors_WithNoEdits.net")
rm(conn.Zn)

# read in zone centriods and remove zones that need to be re-addressed (any zone that has been edited, the RM list)
centers <- read.csv("Orig_ZoneCoordinates.csv",as.is=T)
names(centers)[1] <- "$ZONE:NO"

write.csv(centers[!centers[,1] %in% Rm,],"Orig_ZoneCoordinates_WithNoEdits.csv",row.names=F)

With the two files that are output from these lines "WithNoEdits", the user can drag the "Orig_Connectors_WithNoEdits.net" file into an open Visum SWIM network and have the old connectors re-established. And similarly, can open "Orig_ZoneCoordinates_WithNoEdits.csv", add a header row, alter cells A1 and A2:

Cell A1 = $VISION
Cell A2 = $ZONE:NO

And then paste the connector centriod update into the Visum zone listing table to bring those original zone centriods back to their original positions. The remaining work of connecting revised zones requires manual review. In theory the connectors could be automatically generated or scripted, but best practice is to review each zone and hand place the connectors as a lot of judgement goes into properly placing connectors. Judgment on elements like:

• How is the activity in the zone centered?
• What are the main access points to the zone?
• Are there any physical boundaries (rivers, rail, mts.) blocking certain ways into the zone?

Additionally for SWIM, since there is a truck and a personal vehicle network, special care and consideration needs to be taken to ensure that all zones are connected so that both cars (vehicle class "a") and trucks (vehicle class "d") have access to all other zones (the full network). In some rural and forest zones, the analyst may need to display the network that is available to trucks ("d") to ensure that they are building connections to for all zones that have access to both cars and trucks.

Another tip, applying a filter to the zones is a good way to see how many zones have not yet been connected to the highway network. In this case the zones can be tagged (by applying a filter of zones with zero connectors and giving all those zones an “AddVal1” value of 1). Then each zone identified is reviewed by hand and is reconnected and the center shifted as needed. A record of the changes can then be exported using the AddVal1 flag and ported over to the future year (like 2016). This is done using a similar methodology as the first step. In this case the connectors are filtered to only flagged zones, then those connectors can be exported to a net file and the edited centriods can be saved into a csv. Then the steps above can be followed to copy those edits into other future versions of the SWIM network, so the hand edits / review only has to occur once.

However, before moving the connectors over, an important step is to recalculate the lengths and travel time on all connectors. The way SWIM has been established, any zone with only one connector will uses the length of the connector to calculate the travel time for all the modes. Zones with more than one connector have the connectors travel time zeroed out. This is because zones with multiple connectors are typically representing a zone that has two significantly different population centers or access points and having a travel time to the center of that zone does not correctly represent the real travel time to that population (activity center). Noting, for the most part zones with only one connector have had their zone centroid moved to the activity center of that zone, so having a travel time on the connector is correct in those cases. Before the connectors are finalized (and potentially moved to a future version file) ensure that zones that are supposed to have travel time and those that aren't are properly coded.

Updating and Iterating Seed Skims Until Stable

One of the major "boot-strapped" inputs to SWIM is the large set of skims required by the land use model AA. In the outputs/t19 folder there are 84 matrices that are needed to get the model running. Of these 84 matrices, 58 are at the alpha zone level and the remaining 26 are the beta zone level. So if only alpha zones have been changed, only 58 matrices need to be updated (said with a smile). One could approach this by using Visum (the traffic assignment software in SWIM) to create new free flow travel matrices for the area, and using those to seed. The issue with this method is that the land use model sees no congestion either for the whole system or just those altered zones. No congestion for the entire model area would create a very poor starting point for the model, so the analyst will need to keep track of which zones were altered so that only those cells in the old matrices get added or updated, leaving all unchanged zones with their original congested travel times. The other aspect that complicates this, is that there are many matrices for many different modes, and it would be very difficult to try and develop free flow versions for all 58 different matrices, and then switch out the information for only the new zones from new matrices to the old original matrices.

Because of these complications, the method I chose (and am recommending) is to associate each new or altered zone with a near by zone. The matrices are then updated to duplicate rows and cells for new or altered zones with information from near by or similar zones. This information will be rough, but it makes for an easier process and (I would argue) a better starting point for the land use model.

Either way you choose, one still needs to iterate the first two years of the model (currently 2009 and 2010). The way it currently works is that:

• the updated matrices are placed into t19 (2009),
• the model is run until the end of t20 (2010) which generates a complete set of new matrices with correct travel times and distances based on the land use provided,
• that land use wasn't perfect, so the 2010 skims (which are close enough to congested travel in 2009) are input into t19 and the model is run again through t20.

This process is repeated several times, usually 3-5 depending on how much time is available. No specific convergence criteria has been developed at this point given 84 total matrices to be reviewed, so the practice is just to run this loop for as many times as can be done in a week or so (or a reasonable number given available time constraints).

Here is the R coded provided to update the 58 matrices:

source("readWriteZMXFile.R")

##################
# Alpha zone loop
##################

# Matrix files to work on 
Files <- c(paste(apply(expand.grid(c("op","pk"),c("auto","trk1"),c("fftime","dist","time","toll")),1,paste,collapse=""),".zmx", sep=""),
           paste(apply(expand.grid(c("op","pk"),c("dair"),c("drv","far","fwt","ivt")),1,paste,collapse=""),".zmx", sep=""),
           paste(apply(expand.grid(c("op","pk"),c("wicr"),c("awk","ewk","far","fwt","ivt","twt","xwk")),1,paste,collapse=""),".zmx", sep=""),
           paste(apply(expand.grid(c("op","pk"),c("wt"),c("awk","brd","ewk","far","fwt","ivt","twt","xwk")),1,paste,collapse=""),".zmx", sep=""),
           paste(apply(expand.grid(c("op","pk"),c("wltf"),c("ivt","ovt")),1,paste,collapse=""),".zmx", sep=""))

# Alphas to remove
Rm <- c(4001,4004,4010,4011,4026,4071,4072,4073,4093,4094,4096,4140,4141, # halo cut
        3068,3076,3069,3058,3059,3063,3066,3065,3070,3906,3910,3902,2913,2159, # combined list
        2917,2902,2904,2939,2935,2916,2936,2941,2937,2947,2943,2948,2949,2183,
        3909,2923,2919,2933,2928,2924,2931,2932,2801,2813,2814,2808,2803,2575,
        2809,2081,2079,3276,2083,2651,2649,2629,2623,2626,2606,2635,2636,303,
        2625,2612,2152,2142,2128,2127,2133,2135,2013,1387,2673,2880,2171,
        3848,3819,3816,3898,2844,3847,3241,3193,3291,3273,2864,2645, # Beta merges  
        2818,3055,2338,1470,2670,2861,2846,2779) # Remove the hand moved betas

# Alphas re-numbered
Rn <- as.character(c(2818,3055,2338,1470,2670,2861,2846,2779))
names(Rn) <-       c(3818,3855,3838,3870,3873,3861,3846,2780)

# linkage from old zones to the new splits
Splits <- read.csv("/SplitZonesRecord.csv", as.is=T)
# Only capture new zone numbers
Splits <- Splits[Splits[,1] != Splits[,2],1:2]


for(f in Files){
   x <- readZipMat(paste("t19_Update/",f,sep=""))

   # clear deleted Betas:
   x <- x[!rownames(x) %in% as.character(Rm),!colnames(x) %in% as.character(Rm)]

   # create new zone name vector
   nameUpdate <- rownames(x)
   names(nameUpdate) <- rownames(x)
   nameUpdate[names(Rn)] <- Rn
   
   # switch beta zone numbering to updated vector
   rownames(x) <- colnames(x) <- nameUpdate
   
   # verify that none of the new zone numbers exist in the previous zone system
   Splits[as.character(Splits[,2]) %in% rownames(x),]
   # duplicate all split rows
   sum(!as.character(Splits[,1]) %in% rownames(x))
   
   # Create a full set of new names
   FullNames <- c(nameUpdate, as.character(Splits[,1]))
   names(FullNames) <- c(nameUpdate, as.character(Splits[,2]))   
   FullNames <- FullNames[order(as.numeric(names(FullNames)))]
   
   # build new Full Matrix and rename
   x <- x[FullNames,FullNames]
   rownames(x) <- colnames(x) <- names(FullNames)
   
   # Write out final full matrix
   writeZipMat(x,paste("t19_Update2/",f,sep=""))
}

Here's the "SplitZonesRecord.csv" table referenced in the code above:

SplitZone	NewNumber	Area
3267	3267	Pendleton
3267	3065	Pendleton
3277	3277	Pendleton
3277	3276	Pendleton
3277	3076	Pendleton

Read / write functions for ZMX format can be found here: https://github.com/tlumip/tlumip/tree/master/root/scenario/model/helper_scripts/zmx

Updating Other "Seed" Zonal Inputs that Deal with Information by Zone

In addition to the skims in the t19 (2009) input folder, there are also three other "boot strap" inputs that are related to the specific alpha zones. These files help kick-off the land use model; they're supposed to represent previous (base) year information. So even though they are related to specific alpha zones, it wouldn't make sense to try and put this information into the Visum version file with the other zone fields. It would make the SWIM version file unnecessarily complicated. These version files just need to be edited enough to get the land use model to run initially, and then the resulting t20 files can be used as the new t19 files. Given that the beta zones have not been altered, these files will simply be summed up to the beta zone level. So even if the alpha zone information is off, as long as the beta zone information still sums correctly, the land use model should correctly run for the kick-off (t19) year, and then roll along in future years with fully updated information. Here are example steps to update the three files (two in t19 and one in parameters):

# Alphas to remove
Rm <- c(4001,4004,4010,4011,4026,4071,4072,4073,4093,4094,4096,4140,4141, # halo cut
        3068,3076,3069,3058,3059,3063,3066,3065,3070,3906,3910,3902,2913,2159, # combined list
        2917,2902,2904,2939,2935,2916,2936,2941,2937,2947,2943,2948,2949,2183,
        3909,2923,2919,2933,2928,2924,2931,2932,2801,2813,2814,2808,2803,2575,
        2809,2081,2079,3276,2083,2651,2649,2629,2623,2626,2606,2635,2636,303,
        2625,2612,2152,2142,2128,2127,2133,2135,2013,1387,2673,2880,2171,
        3848,3819,3816,3898,2844,3847,3241,3193,3291,3273,2864,2645, # Beta merges  
        2818,3055,2338,1470,2670,2861,2846,2779) # Remove the hand moved betas

# Three "boot strap" alpha zone files that need to be updated in addition to the alpha skims

x <- read.csv("t19_Update/Employment.csv")  
x <- x[!x$Azone %in% Rm,]
names(nameUpdate)[!names(nameUpdate) %in% as.character(unique(x$Azone))]
x$Azone <- as.numeric(nameUpdate[as.character(x$Azone)])
rownames(x) <- x$Azone
FullNames <- FullNames[FullNames %in% rownames(x)]
x <- x[FullNames,]
x$Azone <- as.numeric(names(FullNames))
write.csv(x,"t19_Update2/Employment.csv", row.names=F, quote=F)


x <- read.csv("t19_Update/Increments.csv")
x <- x[!x$AZone %in% Rm,]
names(nameUpdate)[!names(nameUpdate) %in% as.character(unique(x$AZone))]
x$AZone <- as.numeric(nameUpdate[as.character(x$AZone)])
rownames(x) <- paste(x$AZone,x$FLRName)
FullNames.Inc <- FullNames[FullNames %in% unique(x$AZone)]
FullNames.Inc <- apply(expand.grid(FullNames.Inc, unique(x$FLRName)),1, paste, collapse=" ")
x <- x[FullNames.Inc,]
x$AZone <- as.numeric(rep(names(FullNames),length(unique(x$FLRName))))
write.csv(x,"t19_Update2/Increments.csv", row.names=F, quote=F)
rm(FullNames.Inc)

x <- read.csv("parameters_Update/ActivityConstraintsIWorldMarkets.csv", as.is=T)
x <- x[x$azone>5999,]
write.csv(tapply(x$quantity, list(x$azone,x$activity),sum), "ActivityConstraintsIWorldMarkets_2UpdateVisumPOIs.csv")

After the model has been run once, "Employment.csv" and "Increments.csv" should be copied from the successful t20 outputs folder back as the starting point for t19. "ActivityConstraintsIWorldMarkets.csv" is zero for all internal zones and is actually an old legacy file that will at some point be deleted. The World Market Activity Constraint data is now stored in Visum under the World Market POI table (#7). The code above help produce the file - "ActivityConstraintsIWorldMarkets_2UpdateVisumPOIs.csv" which can then be altered and copied into the Visum POI table. To paste into Visum first insert a blank header row and a blank column and update the cells A1 and A2 accordingly:

Cell A1 = $VISION
Cell A2 = $POIOFCAT_7:CATNO

Change the remaining first column (A) values to all 7 (should be 6 7s, one for each of the 6001-6006 world markets), and then change Cell B2 to "NO" and change value 6001 to 1, 6002, to 2, and so on (basically subtract 6000 from the world market value).

Then this altered World Market data gets pasted into the POI World Market (7) listing.

And you've done it, "Now that wasn't such a chore now was it". Hopefully you have checked and triple checked all your steps. If so, you just copy your updated version files, t19 files, and parameters into the appropriate folders in a new SWIM application and press run. Here are the locations just for clarity:

t19 files -> outputs\t19
parameters file -> inputs\parameters
2010 version file -> inputs\t0
20XX version files -> inputs\tXX+10 (t26 is 2016, 16+10)

It Wasn't Enough to Just Edit Alpha, You Need to Alter Beta as Well.

Don't say I didn't warn you. On the face of it, altering Beta's probably won't seem that much worse when reading these instructions, and it doesn't have to be. Minor edits and alterations will likely have little to no impact. However, any sizable removals or re-establishing of beta zone boundaries will likely kick your land use model out of calibration, which is a whole process in itself, not covered on this page.

So your going to cross your fingers and hope that your calibration holds after altering the beta zones. In addition to the intensive process that you have to go through to edit the alpha zone layer (which needs to be completed as well for beta zone edits), here are the minimal steps that need to be completed if the beta zone layer is altered. There actually aren't that many of them:

• Update the 26 beta zones matrices in t19 (in addition to the 58 alpha zone matrices)
• Update two "boot strap" files in t18 (2008), the only two files in t18
• Update one additional file in "parameters"

The example script below shows the steps needed to process each of these files. As was described when discussing alpha zone matrices, read / write functions for ZMX format can be found here: https://github.com/tlumip/tlumip/tree/master/root/scenario/model/helper_scripts/zmx

# Alex Bettinardi
# 3-31-17

# This script is used to re-establish all the skims for the new zone structure
# first load the zmx reader functions
source("readWriteZMXFile.R")

# optional - read new zone structure alpha-beta crosswalk for checking and reference
a2b <- read.csv("Revised_Alpha/NewAlphaAttributes.csv", as.is=T)
newBs <- sort(unique(a2b$BZONE))


##########################
# Beta zone loop
#####################

# Matrix files to work on 
Files <- c(paste("beta", apply(expand.grid(c("op","pk"),c("auto","trk1"),c("fftime","dist","time","toll")),1,paste,collapse=""),".zmx", sep=""),
           paste(c("b4","b5","b8","c4","o4","s4","w1","w4","w7","w8"),"mcls_beta.zmx",sep=""))

# Betas to remove
Rm <- c(4001,4004,4010,4011,4026,4071,4072,4073,4093,4094,4096,4140,4141)

# Betas re-numbered
Rn <- as.character(c(1330,1523,2062,2101,2104,2480,2757,2145,2019))
names(Rn) <-       c(1318,1505,1736,1774,1775,2475,2727,2122,2009)

for(f in Files){
   x <- readZipMat(paste("t19/",f,sep=""))

   # clear deleted Betas:
   x <- x[!rownames(x) %in% as.character(Rm),!colnames(x) %in% as.character(Rm)]

   # create new zone name vector
   nameUpdate <- rownames(x)
   names(nameUpdate) <- rownames(x)
   nameUpdate[names(Rn)] <- Rn
   
   # switch beta zone numbering to updated vector
   rownames(x) <- colnames(x) <- nameUpdate
   
   # quick check to make sure no new Betas are missing
   if(length(newBs[!newBs %in% nameUpdate])>0) stop("You have new betas missing from your beta matrix")
   #nameUpdate[!nameUpdate %in% newBs]
   
   writeZipMat(x,paste("t19_Update/",f,sep=""))
}

# Three "boot strap" beta zone files that need to be updated inaddition to the beta skims

x <- read.csv("t18/ActivityLocations.csv")
# save total activity to factor NED
Tots <- tapply(x$Quantity, x$Activity,sum)

# keep processing Activity Locations
x <- x[!x$ZoneNumber %in% Rm,]
names(nameUpdate)[!names(nameUpdate) %in% as.character(unique(x$ZoneNumber))]
x$ZoneNumber <- as.numeric(nameUpdate[as.character(x$ZoneNumber)])
# speical format for -Infinity
x$TechnologyLogsum <- gsub("-Inf","-Infinity",x$TechnologyLogsum)
x$LocationUtility <- gsub("-Inf","-Infinity",x$LocationUtility)
# special processing for fields that require a double
TotsRm <- tapply(x$Quantity, x$Activity,sum)
x$Quantity[!(1:nrow(x) %in% grep("\\.",as.character(x$Quantity)))] <- paste(x$Quantity[!(1:nrow(x) %in% grep("\\.",as.character(x$Quantity)))],".0",sep="")
x$ConstraintValue[!(1:nrow(x) %in% grep("\\.",as.character(x$ConstraintValue)))] <- paste(x$ConstraintValue[!(1:nrow(x) %in% grep("\\.",as.character(x$ConstraintValue)))],".0",sep="")
x$SizeUtility <- paste(x$SizeUtility,".0",sep="")
x$ZoneConstant <- paste(x$ZoneConstant,".0",sep="")
x$Size <- paste(x$Size,".0",sep="")
write.csv(x,"t18_Update/ActivityLocations.csv", row.names=F, quote=F)

x <- read.csv("t18/ExchangeResults.csv")
x <- x[!x$ZoneNumber %in% Rm,]
names(nameUpdate)[!names(nameUpdate) %in% as.character(unique(x$ZoneNumber))]
x$ZoneNumber <- as.numeric(nameUpdate[as.character(x$ZoneNumber)])
write.csv(x,"t18_Update/ExchangeResults.csv", row.names=F, quote=F)

x <- read.csv("parameters/ExchangeImportExportI.csv")
x <- x[!x$ZoneNumber %in% Rm,]
names(nameUpdate)[!names(nameUpdate) %in% as.character(unique(x$ZoneNumber))]
x$ZoneNumber[x$ZoneNumber>0] <- as.numeric(nameUpdate[as.character(x$ZoneNumber[x$ZoneNumber>0])])
write.csv(x,"parameters_Update/ExchangeImportExportI.csv", row.names=F, quote=F)

Other Elements (files) that May Need to Be Revised

Again, for minor edits Beta zone changes might not be that much larger of a lift than alpha zones edits, just a small handful of additional files that need editing. However, if the changes are large enough, other changes may be needed. And the bigger the change, the more likely the land use model needs to be re-calibrated, which is the extra work that one would typically want to avoid.

If the edits to the beta zone level include removal of sections of the model, which the examples throughout include, then the overall economic inputs to the model (economic inputs at the model region level) need to be changed (in this case reduced). There are 5 NED files in inputs\t0 that will need to be edited. Below is a section of script that re tabulates the five NED files. Note that this script builds off of the beta zone section of script above, so one would add this to that script if NED files need to be edited. The "Tots" and "TotsRm" activity regional total R objects are created in those steps above. Those R objects are used to scale the overall economic inputs down or up depending on how the Beta's have been revised.

####
# update NED files:
x <- read.csv("t0_NED/activity_forecast.csv",as.is=T)
x$employment <- round(x$employment*TotsRm[x$activity]/Tots[x$activity])
x$output <- round(x$output*TotsRm[x$activity]/Tots[x$activity])
write.csv(x,"t0_Update/activity_forecast.csv", row.names=F, quote=F)

x <- read.csv("t0_NED/construction_forecast.csv",as.is=T)
x$dollars[x$res_or_non_res=="residential"] <- round(x$dollars[x$res_or_non_res=="residential"]*TotsRm["CNST_res_xxx"]/Tots["CNST_res_xxx"])
x$dollars[x$res_or_non_res=="non_residential"] <- round(x$dollars[x$res_or_non_res=="non_residential"]*TotsRm["CNST_nres_xxx"]/Tots["CNST_nres_xxx"])
write.csv(x,"t0_Update/construction_forecast.csv", row.names=F, quote=F)

x <- read.csv("t0_NED/government_forecast.csv",as.is=T)
x$dollars <- round(x$dollars*TotsRm[x$activity]/Tots[x$activity])
write.csv(x,"t0_Update/government_forecast.csv", row.names=F, quote=F)

x <- read.csv("t0_NED/population_forecast.csv",as.is=T)
x$population <- round(x$population*sum(TotsRm)/sum(Tots))
write.csv(x,"t0_Update/population_forecast.csv", row.names=F, quote=F)

x <- read.csv("t0_NED/trade_forecast.csv",as.is=T)
TraAct <- x$trade_activity
TraAct <- gsub(" ", "_", TraAct)
unique(TraAct)[!unique(TraAct) %in% names(Tots)]
x$dollars <- round(x$dollars*TotsRm[TraAct]/Tots[TraAct])
write.csv(x,"t0_Update/trade_forecast.csv", row.names=F, quote=F)

Similar to the process for just alpha revision, in the case of impacting both beta zones and alpha zones the analyst needs to check and triple check all their steps. And then just copy the updated version files, t19 files, parameter files, t18, and t0 files into the appropriate folders in a new SWIM application and press run. Here are the locations again for clarity:

t18 files -> outputs\t18
t19 files -> outputs\t19
parameters file -> inputs\parameters
NED files (t0) -> inputs\t0
2010 version file -> inputs\t0
20XX version files -> inputs\tXX+10 (t26 is 2016, 16+10)

But Wait, Something is Still Not Working

I Warned You. I tried to tell you that editing Beta zones is a thread you really don't want to pull. Enough "I told you so"; on to solutions.

Unanticipated Issues with ALD

One thing that may have been impacted is changing or removing an ALD Regions. If you reduced the SWIM boundary enough you may have removed an ALD region. If that's the case you will need to identify which region was removed and then remove it from these 6 ALD inputs (all of these files are found in inputs\parameters):

• ALDAsc1.RgFn.csv
• ALDAsc1.RgFr.csv
• ALDAsc2.RgFn.csv
• ALDAsc2.RgFr.csv
• ALDAsc3.RgFc.csv
• ALDAsc4.RgFc.csv

However, worse than a cold hard model crash because a region no longer exists, is not having the model crash. The coefficients and calibration factors in these ALD parameters were carefully adjusted based on the original information from the land use model and land inventory from the unaltered zone structure. If the edits of beta zones impacted the ALD regions as well, then there will likely need to be a step to revisit the parameters in these files and make sure that the model is still producing sensible results. The simplest way to do this is to generate an html reference summary using swimr, and try and assess if rents or land development seems to be behaving as anticipated. Any further instruction beyond that will require another wiki page, probably just as long as this one, as well as better documentation on ALD in general.

Unanticipated Issues with AA

Unanticipated probably isn't the right word. If you are impacting the beta zone boundaries, you should anticipate that some issues in the AA model will occur. The team has not gone through this process enough to have a full accounting of all the issues and check lists to consider, but here is a least one that we are aware of. The following steps relate to an issue where the assumed internal office workforce for employment types like construction is no longer in balance. These office employment totals are calculated inside AA (they are not direct inputs), so it can be tricky to figure out what to adjust.

Here is an error similar to the one you might find in "model_run_output.txt":

27-Sep-2018 13:24, INFO, 27-Sep-2018 13:24, ERROR, All zones constrained for RES_offc_off but constraints do not match total within 0.1%, this should be fixed so that ActivityConstraintsI total up to match ActivityTotalsI
27-Sep-2018 13:24, INFO, 27-Sep-2018 13:24, ERROR, All zones constrained for ENGY_offc_off but constraints do not match total within 0.1%, this should be fixed so that ActivityConstraintsI total up to match ActivityTotalsI
27-Sep-2018 13:24, INFO, 27-Sep-2018 13:24, ERROR, All zones constrained for CNST_offc_off but constraints do not match total within 0.1%, this should be fixed so that ActivityConstraintsI total up to match ActivityTotalsI
27-Sep-2018 13:24, INFO, 27-Sep-2018 13:24, ERROR, Allocating 1.5429881270963593E9 of MFG_offc_off amongst unconstrained zones, but utility of every unconstrained zone is negative infinity
27-Sep-2018 13:24, INFO, 27-Sep-2018 13:24, ERROR, This can happen if the unconstrained zones have no suitable space, in which case the solution is to enter a constraint of zero in any zone with no suitable space
27-Sep-2018 13:24, INFO, 27-Sep-2018 13:24, ERROR, The first unconstrained zone is MFG_offc_off in PECASZone:

Here are the steps if this is occurring:

• The files ActivityConstrainsI.csv and ActivityTotalsW.csv are needed to find the imbalance, however, you can't get these files without running AA. Therefore, first go through all the steps above (which you likely already have) and generate all the input files needed and run them through the model and allow it to crash.
• The model will crash in AA t19 and generate the two files above. Using ActivityConstrainsI.csv the t0 NED files need to be recreated using the controls out of ActivityConstraints (script history missing, but from the beta changes script above, ActivityConstrainsI.csv needs to be used to calculate TotsRm, at least that's what my notes say).
• The model is allowed to run and crash again, which again creates the two files above. Both of those new files are used to build the file - t0_VisumFieldUpdate.csv - and the data in that files is pasted into the version files (using the copy and paste process described above).

Here is the script that helps execute that last step:

# AA failed in t19 -  "ActivityConstraintsI.csv" and "ActivityTotalsW.csv" are copied over 
# Or this can be run in the t19 folder of a failed run.

Fields2Balance <- c("RES_offc_off","ENGY_offc_off","CNST_offc_off","MFG_offc_off","WHSL_offc_off",
                    "RET_stor_off","TRNS_trns_off","INFO_info_off","GOV_offc_off","K12_k12_off", "UTL_othr_off")

y <- read.csv("ActivityConstraintsI.csv",as.is=T)
TotsUpdate <- tapply(y$quantity, y$activity,sum)
TotsUpdate <- TotsUpdate[Fields2Balance]

x <- read.csv("ActivityTotalsW.csv",as.is=T)
x <- x[x[,1] %in% Fields2Balance,]
rownames(x) <- x[,1]
TotsUpdate <- x[names(TotsUpdate),2]/TotsUpdate

Zn <- sort(unique(y$Azone))
Zn <- as.character(Zn[Zn<6000])

Out <- sapply(Fields2Balance, function(x) { z=y[y$activity == x,3]; names(z) = y[y$activity == x,1]; z = z[Zn]; z*TotsUpdate[x]   })
write.csv(Out,"t0_VisumFieldUpdate.csv")

The updated version files are pasted into the SWIM directory and the model is run again. This time it should run successfully - although AA may not converge - still working on how to correct that.

One additional step to try is helping AA get to converged prices. If you are able to get AA to run through but not converge, you can look through the model run sequence and find the closest year where AA did converge. As an example, at the time of writing this t19 (the first year) was not converging, but t20 (the following year) did successfully converge to a solution. In this case the t20 ExchangeResults.csv output file from AA will contain the prices that AA t20 ended on. These prices should be close to the successful prices for the nearby year that didn't converge, in this case t19. The t20 ExchangeResults.csv file was copied to t19 and renamed ExchangeResultsI.csv so that AA knows to use that file as the starting point for price adjustment. With that starting point, ideally AA will converge. In this case t19 did converge after a handful of iterations. The user then has the option of taking the converged t19 ExchangeResults.csv file and coping it and renaming it in the t19 folder as ExchangeResultsI.csv. This will give t19 the end price results for that specific year, so AA should reach convergence almost instantly, saving time for future model runs.

New NED Files Added into the Mix

When a zone structure change is made that impacts the overall size of the model (shrinks or grows the model boundary) than, in addition to all the steps above, the five NED files should be regenerated. The specifics on how to regenerate the five NED "forecast" files and specifics on NED in general can be found here.

So the five NED files have been regenerated. There's a good chance that activity forecasts totals for the model area (from NED) will not match up with the dollars (activity) by zone stored in Visum. The zonal totals for the base year (currently 2009) need to be summed up compared to the NED activity totals, and then factored through all the zones, so that the zone totals match the model totals from NED in activity_forecast.csv. Below is some example code to do this comparison. It assumes that the user has exported the 52 activity fields out of Visum into the file referenced as "PreviousAlphaDollars.csv":

# Alex Bettinardi
# 9-27-18

# Script to re-factor the 2009 zonal totals in SWIM to align with the NED inputs

# read in existing dollars by zone
Dol.ZnCo <- read.csv("PreviousAlphaDollars.csv", as.is=T,row.names=1)
# total up dollars by commodity for the model region
Dol.Co <- colSums(Dol.ZnCo)

# read in the NED dollar forcast
Act.. <- read.csv("activity_forecast.csv", as.is=T)
Act.. <- Act..[Act..$year==2009,]
Act. <- Act..$output
names(Act.) <- toupper(Act..$activity)
rm(Act..)

# Create set to adjust
Act. <- Act.[Act.>0]

# Create the adjusted file back to Visum and write out answer
write.csv(round(sweep(Dol.ZnCo[,names(Act.)], 2, Act./Dol.Co[names(Act.)], "*")), "Adj_AlphaDollars_back2Visum.csv")

The result from that script can be pasted back into the Visum version files and all of the activity totals will be imbalance. However, just as is discussed in Unanticipated Issues with AA above, the office categories will likely be imbalance and will need to be adjusted as well. You will know this is the case if you get a warning similar to the one captured in Unanticipated Issues with AA when trying to run SWIM with the new NED files and the adjusted zonal totals in Visum.

The solution is very similar to the one discussed above. However, above we assumed that NED files needed to be adjusted, in this case we just recalculated the NED files, so we want to retain those and instead adjust the zonal totals. So the only step is the final one, the coded required is provided above.

That code produces "t0_VisumFieldUpdate.csv". The result in that file can be pasted into Visum after add a header row is added and the cells A1 and A2 are updated accordingly:

Cell A1 = $VISION
Cell A2 = $ZONE:NO

With any luck your model is now running after those significant (or minor) zone improvements, and hopefully AA is converging. If not, maybe another wiki page has been added that can help you, otherwise... contact HBA.