diff --git a/.gitignore b/.gitignore
index 5b6a065..99cb7da 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,4 @@
-.Rproj.user
-.Rhistory
-.RData
-.Ruserdata
+*.Rproj.user
+*.Rhistory
+*.RData
+*.Ruserdata
diff --git a/R/INAscene.R b/R/INAscene.R
index ddd7221..fd9da6e 100644
--- a/R/INAscene.R
+++ b/R/INAscene.R
@@ -3,7 +3,7 @@
 #'
 #' This function implements and summarizes multiple simulations across a designated range of parameter values
 #'
-#' Updated 2020-10-13
+#' Updated 2021-10-22
 
 #' @param nreals number of realizations to be evaluated
 #' @param ntimesteps number of time steps to be evaluated
@@ -12,6 +12,8 @@
 
 #' @param readgeocoords if T, read in geocoords - otherwise, generate it in each realization
 #' @param geocoords matrix of x,y coordinates of nodes
+#' @param roaddistfilepath filepath to a R Data (.RData extension) file containing an adjacency matrix with the distances via the road network between nodes in geocoords - mutually exclusive with roadtimefilepath
+#' @param roadtimefilepath filepath to a R Data (.RData extension) file containing an adjacency matrix with the travel times via the road network between nodes in geocoords - mutually exclusive with roaddistfilepath
 
 #' @param numnodes (genlocs) the number of nodes, can be a vector of different numbers of nodes for scenario comparisons (numnodes can be entered as a vector to evaluate multiple scenarios)
 #' @param xrange (genlocs) range of x coordinates, e.g., c(0,50)
@@ -107,7 +109,7 @@
 
 
 
-INAscene <- function(nreals, ntimesteps, doplot=F, outputvol='more', readgeocoords, geocoords=NA, numnodes=NA, xrange=NA, yrange=NA, randgeo=NA, readinitinfo, initinfo=NA, initinfo.norp=NA, initinfo.n=NA, initinfo.p=NA, initinfo.dist=NA, readinitbio, initbio=NA, initbio.norp=NA, initbio.n=NA, initbio.p=NA,  initbio.dist=NA, readseam, seam=NA, seamdist=NA, seamrandp=NA, seampla=NA, seamplb=NA, readbpam, bpam=NA, bpamdist=NA, bpamrandp=NA, bpampla=NA, bpamplb=NA, readprobadoptvec, probadoptvec=NA, probadoptmean=NA, probadoptsd=NA, readprobestabvec, probestabvec=NA, probestabmean=NA, probestabsd=NA, maneffdir=NA, maneffmean=NA, maneffsd=NA, usethreshman, maneffthresh=NA, sampeffort=NA) {
+INAscene <- function(nreals, ntimesteps, doplot=F, outputvol='more', readgeocoords, geocoords=NA, roaddistfilepath=NA, roadtimefilepath=NA, numnodes=NA, xrange=NA, yrange=NA, randgeo=NA, readinitinfo, initinfo=NA, initinfo.norp=NA, initinfo.n=NA, initinfo.p=NA, initinfo.dist=NA, readinitbio, initbio=NA, initbio.norp=NA, initbio.n=NA, initbio.p=NA,  initbio.dist=NA, readseam, seam=NA, seamdist=NA, seamrandp=NA, seampla=NA, seamplb=NA, readbpam, bpam=NA, bpamdist=NA, bpamrandp=NA, bpampla=NA, bpamplb=NA, readprobadoptvec, probadoptvec=NA, probadoptmean=NA, probadoptsd=NA, readprobestabvec, probestabvec=NA, probestabmean=NA, probestabsd=NA, maneffdir=NA, maneffmean=NA, maneffsd=NA, usethreshman, maneffthresh=NA, sampeffort=NA) {
 
 ### error message(s)
 
@@ -311,7 +313,7 @@ for (j22 in 1:ncombs) {
   }
 
 
-  temp <- multsame2(nreals2=nreals, ntimesteps2=ntimesteps, doplot2=doplot, readgeocoords2=readgeocoords, geocoords2=geocoords, numnodes2=numnodes, xrange2=xrange, yrange2=yrange, randgeo2=randgeo, readinitinfo2=readinitinfo, initinfo2=initinfo, initinfo.norp2=initinfo.norp, initinfo.n2=initinfo.n, initinfo.p2=initinfo.p, initinfo.dist2=initinfo.dist, readinitbio2=readinitbio, initbio2=initbio, initbio.norp2=initbio.norp, initbio.n2=initbio.n, initbio.p2=initbio.p, initbio.dist2=initbio.dist, readseam2=readseam, seam2=seam,  seamdist2=seamdist, seamrandp2=seamrandp, seampla2=seampla, seamplb2=seamplb, readbpam2=readbpam, bpam2=bpam, bpamdist2=bpamdist, bpamrandp2=bpamrandp, bpampla2=bpampla, bpamplb2=bpamplb, readprobadoptvec2=readprobadoptvec, probadoptvec2=probadoptvec, probadoptmean2=probadoptmean, probadoptsd2=probadoptsd, readprobestabvec2=readprobestabvec, probestabvec2=probestabvec, probestabmean2=probestabmean, probestabsd2=probestabsd, maneffdir2=maneffdir, maneffmean2=maneffmean, maneffsd2=maneffsd, usethreshman2=usethreshman, maneffthresh2=maneffthresh, sampeffort2=sampeffort)
+  temp <- multsame2(nreals2=nreals, ntimesteps2=ntimesteps, doplot2=doplot, readgeocoords2=readgeocoords, geocoords2=geocoords, roaddistfilepath2=roaddistfilepath, roadtimefilepath2=roadtimefilepath, numnodes2=numnodes, xrange2=xrange, yrange2=yrange, randgeo2=randgeo, readinitinfo2=readinitinfo, initinfo2=initinfo, initinfo.norp2=initinfo.norp, initinfo.n2=initinfo.n, initinfo.p2=initinfo.p, initinfo.dist2=initinfo.dist, readinitbio2=readinitbio, initbio2=initbio, initbio.norp2=initbio.norp, initbio.n2=initbio.n, initbio.p2=initbio.p, initbio.dist2=initbio.dist, readseam2=readseam, seam2=seam,  seamdist2=seamdist, seamrandp2=seamrandp, seampla2=seampla, seamplb2=seamplb, readbpam2=readbpam, bpam2=bpam, bpamdist2=bpamdist, bpamrandp2=bpamrandp, bpampla2=bpampla, bpamplb2=bpamplb, readprobadoptvec2=readprobadoptvec, probadoptvec2=probadoptvec, probadoptmean2=probadoptmean, probadoptsd2=probadoptsd, readprobestabvec2=readprobestabvec, probestabvec2=probestabvec, probestabmean2=probestabmean, probestabsd2=probestabsd, maneffdir2=maneffdir, maneffmean2=maneffmean, maneffsd2=maneffsd, usethreshman2=usethreshman, maneffthresh2=maneffthresh, sampeffort2=sampeffort)
 
 
   multdetails[[rowcount]] <- temp
@@ -342,4 +344,4 @@ if (outputvol=='more') {
     list(multout=multout)
 }
 
-}
\ No newline at end of file
+}
diff --git a/R/cities.csv b/R/cities.csv
new file mode 100644
index 0000000..bcdeb6f
--- /dev/null
+++ b/R/cities.csv
@@ -0,0 +1,327 @@
+id,lon,lat
+NewYork[d],-73.93,40.66
+LosAngeles,-118.41,34.01
+Chicago,-87.68,41.83
+Houston,-95.39,29.78
+Phoenix,-112.09,33.57
+Philadelphia[e],-75.13,40
+SanAntonio,-98.52,29.47
+SanDiego,-117.13,32.81
+Dallas,-96.76,32.79
+SanJose,-121.81,37.29
+Austin,-97.75,30.3
+Jacksonville[f],-81.66,30.33
+FortWorth,-97.34,32.78
+Columbus,-82.98,39.98
+Indianapolis[g],-86.14,39.77
+Charlotte,-80.83,35.2
+SanFrancisco[h],-123.03,37.72
+Seattle,-122.35,47.62
+Denver[i],-104.88,39.76
+Washington[j],-77.01,38.9
+Nashville[k],-86.78,36.17
+OklahomaCity,-97.51,35.46
+ElPaso,-106.42,31.84
+Boston,-71.02,42.33
+Portland,-122.65,45.53
+LasVegas,-115.26,36.22
+Detroit,-83.1,42.38
+Memphis,-89.97,35.1
+Louisville[l],-85.64,38.16
+Baltimore[m],-76.61,39.3
+Milwaukee,-87.96,43.06
+Albuquerque,-106.64,35.1
+Tucson,-110.87,32.15
+Fresno,-119.79,36.78
+Sacramento,-121.46,38.56
+KansasCity,-94.55,39.12
+Mesa,-111.71,33.4
+Atlanta,-84.42,33.76
+Omaha,-96.04,41.26
+ColoradoSprings,-104.76,38.86
+Raleigh,-78.64,35.83
+LongBeach,-118.15,33.8
+VirginiaBeach[m],-76.02,36.78
+Miami,-80.2,25.77
+Oakland,-122.22,37.76
+Minneapolis,-93.26,44.96
+Tulsa,-95.9,36.12
+Bakersfield,-119.01,35.32
+Wichita,-97.34,37.69
+Arlington,-97.12,32.7
+Aurora,-104.68,39.68
+Tampa,-82.47,27.97
+NewOrleans[n],-89.93,30.05
+Cleveland,-81.67,41.47
+Honolulu[b],-157.84,21.32
+Anaheim,-117.76,33.85
+Lexington[o],-84.45,38.04
+Stockton,-121.31,37.97
+CorpusChristi,-97.17,27.75
+Henderson,-115.03,36
+Riverside,-117.39,33.93
+Newark,-74.17,40.72
+SaintPaul,-93.1,44.94
+SantaAna,-117.88,33.73
+Cincinnati,-84.5,39.14
+Irvine,-117.77,33.67
+Orlando,-81.27,28.41
+Pittsburgh,-79.97,40.43
+St.Louis[m],-90.24,38.63
+Greensboro,-79.82,36.09
+JerseyCity,-74.06,40.71
+Anchorage[p],-149.28,61.17
+Lincoln,-96.68,40.81
+Plano,-96.74,33.05
+Durham,-78.9,35.98
+Buffalo,-78.85,42.89
+Chandler,-111.85,33.28
+ChulaVista,-117.01,32.62
+Toledo,-83.58,41.66
+Madison,-89.42,43.08
+Gilbert[q],-111.74,33.31
+Reno,-119.84,39.54
+FortWayne,-85.14,41.08
+NorthLasVegas,-115.09,36.28
+St.Petersburg,-82.64,27.76
+Lubbock,-101.88,33.56
+Irving,-96.97,32.85
+Laredo,-99.48,27.56
+WinstonSalem,-80.26,36.1
+Chesapeake[m],-76.3,36.67
+Glendale,-112.18,33.53
+Garland,-96.63,32.9
+Scottsdale,-111.86,33.68
+Norfolk[m],-76.24,36.92
+Boise[r],-116.23,43.6
+Fremont,-121.94,37.49
+Spokane,-117.43,47.66
+SantaClarita,-118.5,34.4
+BatonRouge[s],-91.13,30.44
+Richmond[m],-77.47,37.53
+Hialeah,-80.3,25.86
+SanBernardino,-117.29,34.14
+Tacoma,-122.45,47.25
+Modesto,-121,37.63
+Huntsville,-86.67,34.69
+DesMoines,-93.61,41.57
+Yonkers,-73.86,40.94
+Rochester,-77.61,43.16
+MorenoValley,-117.2,33.92
+Fayetteville,-78.97,35.08
+Fontana,-117.46,34.1
+Columbus[t],-84.87,32.51
+Worcester,-71.8,42.26
+PortSt.Lucie,-80.38,27.28
+LittleRock,-92.35,34.72
+Augusta[u],-82.07,33.36
+Oxnard,-119.2,34.2
+Birmingham,-86.79,33.52
+Montgomery,-86.26,32.34
+Frisco,-96.82,33.15
+Amarillo,-101.83,35.19
+SaltLakeCity,-111.93,40.77
+GrandRapids,-85.65,42.96
+HuntingtonBeach,-118,33.69
+OverlandPark,-94.69,38.88
+Glendale,-118.24,34.18
+Tallahassee,-84.25,30.45
+GrandPrairie,-97.02,32.68
+McKinney,-96.66,33.19
+CapeCoral,-81.99,26.64
+SiouxFalls,-96.73,43.53
+Peoria,-112.3,33.78
+Providence,-71.41,41.82
+Vancouver,-122.59,45.63
+Knoxville,-83.94,35.97
+Akron,-81.52,41.08
+Shreveport,-93.79,32.46
+Mobile,-88.1,30.66
+Brownsville,-97.45,25.99
+NewportNews[m],-76.52,37.07
+FortLauderdale,-80.14,26.14
+Chattanooga,-85.24,35.06
+Tempe,-111.93,33.38
+Aurora,-88.29,41.76
+SantaRosa,-122.7,38.44
+Eugene,-123.11,44.05
+ElkGrove,-121.38,38.41
+Salem,-123.02,44.92
+Ontario,-117.6,34.03
+Cary[v],-78.81,35.78
+RanchoCucamonga,-117.56,34.12
+Oceanside,-117.3,33.22
+Lancaster,-118.17,34.69
+GardenGrove,-117.96,33.77
+PembrokePines,-80.34,26.02
+FortCollins,-105.06,40.54
+Palmdale,-118.1,34.59
+Springfield,-93.29,37.19
+Clarksville,-87.34,36.56
+Salinas,-121.63,36.69
+Hayward,-122.1,37.62
+Paterson,-74.16,40.91
+Alexandria[m],-77.08,38.82
+Macon[w],-83.69,32.8
+Corona,-117.56,33.86
+KansasCity[x],-94.74,39.12
+Lakewood,-105.11,39.69
+Springfield,-72.54,42.11
+Sunnyvale,-122.02,37.38
+Jackson,-90.21,32.31
+Killeen,-97.73,31.07
+Hollywood,-80.16,26.03
+Murfreesboro,-86.41,35.85
+Pasadena,-95.15,29.65
+Bellevue,-122.15,47.59
+Pomona,-117.76,34.05
+Escondido,-117.07,33.13
+Joliet,-88.14,41.51
+Charleston,-79.95,32.81
+Mesquite,-96.58,32.76
+Naperville,-88.16,41.74
+Rockford,-89.06,42.25
+Bridgeport,-73.19,41.18
+Syracuse,-76.14,43.04
+Savannah,-81.15,32
+Roseville,-121.31,38.76
+Torrance,-118.34,33.83
+Fullerton,-117.92,33.88
+Surprise,-112.45,33.67
+McAllen,-98.24,26.23
+Thornton,-104.94,39.91
+Visalia,-119.32,36.32
+Olathe,-94.81,38.88
+Gainesville,-82.34,29.67
+WestValleyCity,-112.01,40.68
+Orange,-117.86,33.78
+Denton,-97.14,33.21
+Warren,-83.02,42.49
+Pasadena,-118.13,34.16
+Waco,-97.18,31.56
+CedarRapids,-91.67,41.96
+Dayton,-84.19,39.77
+Elizabeth,-74.19,40.66
+Hampton[m],-76.29,37.04
+Columbia,-80.89,34.02
+Kent,-122.21,47.38
+Stamford,-73.54,41.07
+Lakewood,-74.21,40.08
+Victorville,-117.35,34.52
+Miramar,-80.33,25.97
+CoralSprings,-80.25,26.27
+SterlingHeights,-83.03,42.58
+NewHaven,-72.92,41.31
+Carrollton,-96.89,32.98
+Midland,-102.11,32.02
+Norman,-97.34,35.24
+SantaClara,-121.96,37.36
+Athens[y],-83.37,33.94
+ThousandOaks,-118.87,34.19
+Topeka,-95.69,39.03
+SimiValley,-118.74,34.26
+Columbia,-92.32,38.95
+Vallejo,-122.26,38.1
+Fargo,-96.82,46.86
+Allentown,-75.47,40.59
+Pearland,-95.32,29.55
+Concord,-122,37.97
+Abilene,-99.73,32.45
+Arvada,-105.15,39.83
+Berkeley,-122.29,37.86
+AnnArbor,-83.73,42.27
+Independence,-94.35,39.08
+Rochester,-92.47,44.01
+Lafayette[z],-92.02,30.2
+Hartford,-72.68,41.76
+CollegeStation,-96.29,30.58
+Clovis,-119.68,36.82
+Fairfield,-122.03,38.25
+PalmBay,-80.66,27.98
+Richardson,-96.7,32.97
+RoundRock,-97.66,30.52
+Cambridge,-71.11,42.37
+Meridian,-116.39,43.61
+WestPalmBeach,-80.12,26.74
+Evansville,-87.53,37.98
+Clearwater,-82.76,27.97
+Billings,-108.54,45.78
+WestJordan,-112,40.6
+Richmond,-122.36,37.95
+Westminster,-105.06,39.88
+Manchester,-71.44,42.98
+Lowell,-71.32,42.63
+Wilmington,-77.88,34.2
+Antioch,-121.79,37.97
+Beaumont,-94.14,30.08
+Provo,-111.64,40.24
+NorthCharleston,-80.06,32.91
+Elgin,-88.32,42.03
+Carlsbad,-117.28,33.12
+Odessa,-102.34,31.88
+Waterbury,-73.03,41.55
+Springfield,-89.64,39.79
+LeagueCity,-95.1,29.49
+Downey,-118.13,33.93
+Gresham,-122.44,45.5
+HighPoint,-79.99,35.99
+BrokenArrow,-95.78,36.03
+Peoria,-89.61,40.75
+Lansing,-84.55,42.71
+Lakeland,-81.95,28.05
+PompanoBeach,-80.13,26.24
+CostaMesa,-117.91,33.66
+Pueblo,-104.61,38.26
+Lewisville,-96.98,33.04
+MiamiGardens,-80.24,25.94
+LasCruces,-106.78,32.32
+SugarLand,-95.61,29.59
+Murrieta,-117.19,33.57
+Ventura[aa],-119.25,34.26
+Everett,-122.19,47.95
+Temecula,-117.13,33.49
+Dearborn,-83.21,42.31
+SantaMaria,-120.44,34.93
+WestCovina,-117.9,34.05
+ElMonte,-118.02,34.07
+Greeley,-104.76,40.41
+Sparks,-119.73,39.55
+Centennial,-104.86,39.59
+Boulder,-105.25,40.02
+SandySprings,-84.36,33.93
+Inglewood,-118.34,33.95
+Edison,-74.41,40.52
+SouthFulton,-84.67,33.59
+GreenBay,-87.98,44.52
+Burbank,-118.32,34.19
+Renton,-122.19,47.47
+Hillsboro,-122.93,45.52
+ElCajon,-116.96,32.8
+Tyler,-95.3,32.31
+Davie[ac],-80.28,26.07
+SanMateo,-122.31,37.56
+Brockton,-71.02,42.08
+Concord,-80.6,35.4
+JurupaValley,-117.46,34
+DalyCity,-122.46,37.7
+Allen,-96.66,33.09
+RioRancho,-106.67,35.29
+Rialto,-117.38,34.11
+Woodbridge,-74.28,40.56
+SouthBend,-86.26,41.67
+SpokaneValley,-117.23,47.67
+Norwalk,-118.08,33.9
+Menifee,-117.19,33.69
+Vacaville,-121.97,38.35
+WichitaFalls,-98.52,33.9
+Davenport,-90.6,41.55
+Quincy,-71,42.25
+Chico,-121.83,39.74
+Lynn,-70.95,42.47
+Lee'sSummit,-94.37,38.92
+NewBedford,-70.93,41.64
+FederalWay,-122.35,47.32
+Clinton,-82.92,42.59
+Edinburg,-98.16,26.3
+Nampa,-116.56,43.57
+Roanoke[m],-79.93,37.27
diff --git a/R/genmovnet.R b/R/genmovnet.R
index 13d80db..0f6454a 100644
--- a/R/genmovnet.R
+++ b/R/genmovnet.R
@@ -3,7 +3,7 @@
 #'
 #' This function generates an adjacency matrix for movement, assumed symmetric (in this version).  It is used by functions including \code{INAscene}. The movement adjacency matrix is composed of 1s and 0s only if lktype="pa" option is used
 #'
-#' Updated 2020-09-05
+#' Updated 2021-10-22
 
 #' @param amdist4 the function of distance used to estimate movement probability - 'random' (not related to distance) or 'powerlaw' (inverse power law) or 'exp' (negative exponential, to be added)
 #' @param iplot if T, generates igraph plot of adjacency matrix
@@ -11,6 +11,8 @@
 #' @param amplb4 inverse power law parameter b in ad^(-b)
 #' @param amrandp4 random matrix with entries binomial with probability p
 #' @param geocoords4n the matrix of xy coordinates for node locations, used when the probability of a link is a function of distance (note that the distance between each pair of locations is assumed to be greater than 1)
+#' @param roaddistfilepath4n filepath to a R Data (.RData extension) file containing an adjacency matrix with the distances via the road network between nodes in geocoords - mutually exclusive with roadtimefilepath4n
+#' @param roadtimefilepath4n filepath to a R Data (.RData extension) file containing an adjacency matrix with the travel times via the road network between nodes in geocoords - mutually exclusive with roaddistfilepath4n
 #' @keywords dispersal
 #' @export 
 #' @import igraph
@@ -28,13 +30,24 @@
 #' x9 <- genmovnet(j <- genlocs(numnodes4=300, xrange4 = c(0, 10), yrange4 = c(0, 100)), amdist4='powerlaw', ampla4=2, amplb4=1, iplot=T)
 
 
-genmovnet <- function(geocoords4n, amdist4, iplot=F, amrandp4, ampla4, amplb4){
+genmovnet <- function(geocoords4n, amdist4, roaddistfilepath4n, roadtimefilepath4n, iplot=F, amrandp4, ampla4, amplb4){
 
   dimam <- dim(geocoords4n)[1]
 
   if (amdist4 == 'powerlaw') { # ad^(-b)
 
-    tdist <- as.matrix(dist(geocoords4n, method = "euclidean", diag=T, upper=T))
+    tdist = NA
+    if(!is.na(roaddistfilepath4n)) {
+      #Load Road Distance Data
+      tdist <- roaddata(geocoords5=geocoords4n, roaddistfilepath5=roaddistfilepath4n)
+    } else if(!is.na(roadtimefilepath4n)) {
+      #Load Road Time Data
+      tdist <- roaddata(geocoords5=geocoords4n, roadtimefilepath5=roadtimefilepath4n)
+    } else {
+      tdist <- as.matrix(dist(geocoords4n, method = "euclidean", diag=T, upper=T))
+    }
+    #tdist <- as.matrix(dist(geocoords4n, method = "euclidean", diag=T, upper=T))
+
     linkmat <- ampla4*tdist^(-amplb4)  
   }
 
@@ -77,4 +90,4 @@ genmovnet <- function(geocoords4n, amdist4, iplot=F, amrandp4, ampla4, amplb4){
   }
 
   linkmat
-}
\ No newline at end of file
+}
diff --git a/R/load.R b/R/load.R
new file mode 100644
index 0000000..c36bef0
--- /dev/null
+++ b/R/load.R
@@ -0,0 +1,161 @@
+src = function() {
+  source("INAscene.R");
+  source("estab.R");
+  source("estinfo.R");
+  source("genlocs.R");
+  source("genmovnet.R");
+  source("initvals.R");
+  source("layerplot.R");
+  source("load.R");
+  source("makedec.R");
+  source("multistart.R");
+  source("multsame2.R");
+  source("ntsteps2.R");
+  source("onestart.R");
+  source("reduce_nodes_to_focus.R");
+  source("roaddata.R");
+  source("setup2.R");
+  source("smartsurv.R");
+  source("smartsurv.weight.R");
+  source("spreadstep.R");
+  source("squish.R");
+
+  test();
+}
+
+test = function() {
+  #cities = read.csv("cities.csv");
+  #count = 326;
+  #coords = t(matrix(nrow=2,c(cities$lon[1:count], cities$lat[1:count]), byrow=T));
+
+  #Use Ethiopia data
+  cities = read.csv("farms.csv");
+  count = 1259;
+  coords = t(matrix(nrow=2,c(cities$lon[1:count], cities$lat[1:count]), byrow=T));
+
+  initinfo = round(runif(count, 0, 1));
+  initbio = round(runif(count, 0, 1));
+
+#j25.readgeocoords <- INAscene(nreals=3,
+#  ntimesteps=3,
+#  doplot=F,
+#  readgeocoords=T,
+#  #geocoords=matrix(c(1, 1, 1, 2, 1, 3, 2, 1, 2, 2, 2, 3), byrow=T, ncol=2),
+#  geocoords=coords,
+#  roaddatafilepath="roaddata.RData",
+#  numnodes=NA,
+#  xrange=NA,
+#  yrange=NA,
+#  randgeo=F,
+#  readinitinfo=T,
+#  initinfo=initinfo,
+#  initinfo.norp=NA,
+#  initinfo.n=NA,
+#  initinfo.p=NA,
+#  initinfo.dist=NA,
+#  readinitbio=T,
+#  initbio=initbio,
+#  initbio.norp=NA,
+#  initbio.n=NA,
+#  initbio.p=NA,
+#  initbio.dist=NA,
+#  readseam=F,
+#  seam=NA,
+#  seamdist='powerlaw',
+#  seamrandp=c(0.01, 0.05, 0.1, 0.5),
+#  seampla=1,
+#  seamplb=1,
+#  readbpam=F,
+#  bpam=NA,
+#  bpamdist='random',
+#  bpamrandp=0.1,
+#  bpampla=NA,
+#  bpamplb=NA,
+#  readprobadoptvec=F,
+#  probadoptvec=NA,
+#  probadoptmean=0.1,
+#  probadoptsd=0.1,
+#  readprobestabvec=F,
+#  probestabvec=NA,
+#  probestabmean=0.1,
+#  probestabsd=0.1,
+#  maneffdir='decrease_estab',
+#  maneffmean=0.5,
+#  maneffsd=0.1,
+#  usethreshman=F,
+#  maneffthresh=NA,
+#  sampeffort=NA
+#)
+
+output <- INAscene(
+	nreals=1,
+	ntimesteps=100,
+	doplot=F,
+	outputvol='more',
+	readgeocoords=T,
+	geocoords=coords,
+	roaddistfilepath="roaddata.RData",
+	#roadtimefilepath=NA,
+	numnodes=count,
+	xrange=NA,
+	yrange=NA,
+	randgeo=NA,
+	readinitinfo=F,
+	initinfo=NA,
+	initinfo.norp='num',
+	initinfo.n=0, #May really want to change...
+	initinfo.p=NA,
+	initinfo.dist='random',
+	readinitbio=F,
+	initbio=NA,
+	initbio.norp='num',
+	initbio.n=1,
+	initbio.p=NA,
+	initbio.dist='random',
+	readseam=F,
+	seam=NA,
+	seamdist='powerlaw',
+	seamrandp=NA,
+	seampla=1, #May really want to change later
+	seamplb=1, #May really want to change later
+	readbpam=F,
+	bpam=NA,
+	bpamdist='powerlaw',
+	bpamrandp=NA,
+	bpampla=100, #May really want to change later
+	bpamplb=1, #May really want to change later
+#	readprobadoptvec=F,
+#	probadoptvec=NA,
+#	probadoptmean=1, #May really want to change later
+#	probadoptsd=0, #May really want to change later
+#	readprobestabvec=F,
+#	probestabvec=NA,
+#	probestabmean=1, #May really want to change later
+#	probestabsd=0, #May really want to change later
+#	maneffdir='decrease_estab',
+#	maneffmean=.5, #May really want to change later
+#	maneffsd=0, #May really want to change later
+#	usethreshman=F,
+#	maneffthresh=NA,
+#	sampeffort=NA)
+  readprobadoptvec=F,
+  probadoptvec=NA,
+  #probadoptmean=0.1,
+  #probadoptsd=0.1,
+  probadoptmean=0,
+  probadoptsd=0,
+  readprobestabvec=F,
+  probestabvec=NA,
+  #probestabmean=0.1,
+  probestabmean=.99,
+  #probestabsd=0.1,
+  probestabsd=0,
+  maneffdir='decrease_estab',
+  maneffmean=0.5,
+  maneffsd=0.1,
+  usethreshman=F,
+  maneffthresh=NA,
+  sampeffort=NA
+  )
+return(output);
+}
diff --git a/R/map.R b/R/map.R
new file mode 100644
index 0000000..454adf0
--- /dev/null
+++ b/R/map.R
@@ -0,0 +1,121 @@
+library(rgl)
+library(raster)
+library(viridis)
+library(tiff)
+library(rworldmap)
+
+rotateccw = function(x) t(apply(t(x),2,rev));
+rotatecw = function(x) t(apply(x,2,rev));
+
+rendermap = function(datatype="", infected) {
+	filename = "img/ethiopia";
+	if(datatype != "") {
+		filename = paste(filename, "_", datatype, sep="")
+
+	}
+
+	filename = paste(filename, ".tif", sep="")
+
+	readimg = function(raster, filename) {
+		col_rgb = readTIFF(filename);
+
+		col = matrix(nrow=raster@nrows,ncol=raster@ncols);
+
+		for (x in 1:raster@nrows) {
+			for (y in 1:raster@ncols) {
+				col[x,y] = rgb(col_rgb[x,y,1], col_rgb[x,y,2], col_rgb[x,y,3]);
+			}
+		}
+
+		return(rotatecw(col));
+	}
+
+	africa_map_altitude <- raster("img/ethiopia.tif")
+	ethiopia<-africa_map_altitude/max(getValues(africa_map_altitude),na.rm=TRUE) #rescale between 0 and 1 (otherwise scale is in minutes)
+
+	scale = 5;
+
+	heights = rotateccw(scale * matrix(ethiopia@data@values,nrow=ethiopia@nrows));
+
+	col = readimg(ethiopia, filename);
+
+
+	rgl.clear();
+
+	str(heights);
+
+	surface3d(1:nrow(heights),1:ncol(heights),heights,color=col)
+
+	data("countriesLow") #load map border data
+	countriesLow = subset(countriesLow,LAT>0 & LAT<20 & LON>30 & LON<50) #Filter to only local countries
+
+	for(j in 1:length(countriesLow)) {
+		points = countriesLow@polygons[[j]]@Polygons[[1]]@coords;
+
+		x = ((points[,1] - 30) / 20 * 240);
+		y = ((1 - (points[,2] - 0) / 20) * 240);
+
+		filter = (x > 0 & x < 240 & y > 0 & y < 240);
+		x = x[filter];
+		y = y[filter];
+
+		z = c();
+		for(i in 1:length(x)) {
+			z = c(z, 1.1 * heights[round(x[i]), 240 - round(y[i])]);
+		}
+
+		lines3d(x,240 - y,z,col="black",lwd=4);
+
+	}
+
+	rgl.viewpoint(phi=5, zoom=.75);
+	return(heights);
+}
+
+points = function(infected, heights) {
+	visited = rep(FALSE, length(infected[[1]]));
+
+	for(i in 1:length(infected)) {
+		for(j in 1:length(infected[[i]])){
+			count = 0;
+			if(infected[[i]][[j]] & !visited[[j]]) {
+				x = ((cities$lon[[j]] - 30) / 20 * 240);
+				y = ((1 - (cities$lat[[j]]- 0) / 20) * 240);
+
+				z = 1.1 * heights[round(x), 240 - round(y)];
+
+				spheres3d(x,240 - y,z,col="red",radius=1);
+				count = count + 1;
+
+				visited[[j]] = TRUE;
+			}
+		}
+
+
+
+		filename=paste("/mnt/ram/frame",i,".png",sep="");
+		rgl.snapshot(filename=filename,fmt="png");
+	}
+}
+
+INAPlot = function() {
+	cities = read.csv("farms.csv");
+	count = 1259;
+	coords = t(matrix(nrow=2,c(cities$lon[1:count], cities$lat[1:count]), byrow=T));
+
+	initinfo = round(runif(count, 0, 1));
+	initbio = round(runif(count, 0, 1));
+
+	#res = src();
+	#rgl.open();
+	rgl.clear();
+
+	heights=rendermap();
+	points(infected=res$multdetails[[1]]$multout[[1]]$estabvecL, heights=heights);
+
+}
+
+ll = function() {
+	source("map.R");
+	INAPlot();
+}
diff --git a/R/multsame2.R b/R/multsame2.R
index 099603a..25cd796 100644
--- a/R/multsame2.R
+++ b/R/multsame2.R
@@ -3,7 +3,7 @@
 #'
 #' This function runs and summarizes multiple INA simulations with the same parameter values.  (It uses INA functions estinfo, genlocs, initvals, setup2, genmovnet, spreadstep, makedec, estab, and ntsteps2.)
 #'
-#' Updated 2020-09-05
+#' Updated 2021-10-22
 
 #' @param nreals2 number of realizations to be evaluated
 #' @param usethreshman2 if T, use the threshold for the management effect size estimate; thus, information is never present anywhere unless the management effect estimate exceeds the threshold
@@ -28,6 +28,8 @@
 
 #' @param readgeocoords2 if T, read in geocoords2 - otherwise, generate it in each realization 
 #' @param geocoords2 matrix of x,y coordinates of nodes, read in if readgeocoords2=T
+#' @param roaddistfilepath2 filepath to a R Data (.RData extension) file containing an adjacency matrix with the distances via the road network between nodes in geocoords - mutually exclusive with roadtimefilepath2
+#' @param roadtimefilepath2 filepath to a R Data (.RData extension) file containing an adjacency matrix with the travel times via the road network between nodes in geocoords - mutually exclusive with roaddistfilepath2
 
 #' @param xrange2 (genlocs) range of x coordinates
 #' @param yrange2 (genlocs) range of y coordinates
@@ -85,7 +87,7 @@
 
 
 
-multsame2 <- function(nreals2, usethreshman2, ntimesteps2, geocoords2, readgeocoords2, maneffdir2, maneffmean2, maneffsd2, maneffthresh2, sampeffort2, xrange2, yrange2, numnodes2, randgeo2, readinitinfo2, initinfo2, initinfo.dist2, initinfo.n2, initinfo.norp2, initinfo.p2, readinitbio2, initbio2, initbio.dist2, initbio.n2, initbio.norp2, initbio.p2, readseam2, seam2, seamdist2, seampla2, seamplb2, seamrandp2, readbpam2, bpam2, bpamdist2, bpampla2, bpamplb2, bpamrandp2, probadoptmean2, probadoptsd2, probadoptvec2, readprobadoptvec2, probestabmean2, probestabsd2, readprobestabvec2, probestabvec2, doplot2=F){
+multsame2 <- function(nreals2, usethreshman2, ntimesteps2, geocoords2, roaddistfilepath2, roadtimefilepath2, readgeocoords2, maneffdir2, maneffmean2, maneffsd2, maneffthresh2, sampeffort2, xrange2, yrange2, numnodes2, randgeo2, readinitinfo2, initinfo2, initinfo.dist2, initinfo.n2, initinfo.norp2, initinfo.p2, readinitbio2, initbio2, initbio.dist2, initbio.n2, initbio.norp2, initbio.p2, readseam2, seam2, seamdist2, seampla2, seamplb2, seamrandp2, readbpam2, bpam2, bpamdist2, bpampla2, bpamplb2, bpamrandp2, probadoptmean2, probadoptsd2, probadoptvec2, readprobadoptvec2, probestabmean2, probestabsd2, readprobestabvec2, probestabvec2, doplot2=F){
 
 
   # determine the number of nodes, 
@@ -114,12 +116,11 @@ multsame2 <- function(nreals2, usethreshman2, ntimesteps2, geocoords2, readgeoco
 
   for(j in 1:nreals2) {
 
-
     tempo <- setup2(maneffmean3s=maneffmean2, maneffsd3s=maneffsd2, maneffthresh3=maneffthresh2, sampeffort3=sampeffort2, usethreshman3=usethreshman2, readgeocoords3s=readgeocoords2, geocoords3s=geocoords2, xrange3=xrange2, yrange3=yrange2, numnodes3=numnodes2, randgeo3=randgeo2, readinitinfo3=readinitinfo2, initinfo3=initinfo2, initinfo.dist3=initinfo.dist2, initinfo.n3=initinfo.n2, initinfo.norp3=initinfo.norp2, initinfo.p3=initinfo.p2, readinitbio3=readinitbio2, initbio3=initbio2, initbio.dist3=initbio.dist2, initbio.n3=initbio.n2, initbio.norp3=initbio.norp2, initbio.p3=initbio.p2, plotmp=doplot2)
 
 # note that the following uses the output object tempo
 
-    temp2 <- ntsteps2(nsteps=ntimesteps2, infon=tempo$com.yes, geocoords3n=tempo$geocoords3s, vect1cn=tempo$infovec, vect1dn=tempo$estabvec, readseam3=readseam2, seam3=seam2, seamdist3=seamdist2, seampla3=seampla2, seamplb3=seamplb2, seamrandp3=seamrandp2, readbpam3=readbpam2, bpam3=bpam2, bpamdist3=bpamdist2, bpampla3=bpampla2, bpamplb3=bpamplb2, bpamrandp3=bpamrandp2, readprobadoptvec3=readprobadoptvec2, probadoptvec3=probadoptvec2, probadoptmean3=probadoptmean2, probadoptsd3=probadoptsd2, probestabmean3=probestabmean2, probestabsd3=probestabsd2, maneffdir3=maneffdir2, maneffmean3n=maneffmean2, maneffsd3n=maneffsd2, readprobestabvec3=readprobestabvec2, probestabvec3=probestabvec2, plotmpn=doplot2)
+    temp2 <- ntsteps2(nsteps=ntimesteps2, infon=tempo$com.yes, geocoords3n=tempo$geocoords3s, roaddistfilepath3n=roaddistfilepath2, roadtimefilepath3n=roadtimefilepath2, vect1cn=tempo$infovec, vect1dn=tempo$estabvec, readseam3=readseam2, seam3=seam2, seamdist3=seamdist2, seampla3=seampla2, seamplb3=seamplb2, seamrandp3=seamrandp2, readbpam3=readbpam2, bpam3=bpam2, bpamdist3=bpamdist2, bpampla3=bpampla2, bpamplb3=bpamplb2, bpamrandp3=bpamrandp2, readprobadoptvec3=readprobadoptvec2, probadoptvec3=probadoptvec2, probadoptmean3=probadoptmean2, probadoptsd3=probadoptsd2, probestabmean3=probestabmean2, probestabsd3=probestabsd2, maneffdir3=maneffdir2, maneffmean3n=maneffmean2, maneffsd3n=maneffsd2, readprobestabvec3=readprobestabvec2, probestabvec3=probestabvec2, plotmpn=doplot2)
 
     # save the output from setup2 and ntsteps2
     multout[[j]] <- temp2
@@ -152,4 +153,4 @@ multsame2 <- function(nreals2, usethreshman2, ntimesteps2, geocoords2, readgeoco
 
 
   list(multout=multout, setocom=setocom, setodec=setodec, setodisp=setodisp, setoestab=setoestab, meancom=meancom, meandec=meandec, meandisp=meandisp, meanestab=meanestab, com5=com5, dec5=dec5, disp5=disp5, estab5=estab5, com95=com95, dec95=dec95, disp95=disp95, estab95=estab95)
-}
\ No newline at end of file
+}
diff --git a/R/ntsteps2.R b/R/ntsteps2.R
index aa2f4aa..91ec250 100644
--- a/R/ntsteps2.R
+++ b/R/ntsteps2.R
@@ -4,7 +4,7 @@
 #' This function gives the info status and establishment status at each node after multiple time steps. Its use follows use of the function \code{setup2}, or another source of node locations and initial conditions.  (Note that communication and dispersal status are given for the beginning of the time step (rather than for the end), so the corresponding output for these has length equal to the number of time steps plus 1.  In contrast, adoption and establishment output is of length equal to the number of time steps.) It uses the following functions, as needed: genmovnet, spreadstep, makedec, and estab. It draws on output from function setup2 in terms of whether communication occurs (estinfo), generated geographic locations of nodes as needed (genloc), and the initial locations for management information and the bioentity.
 
 #'
-#' Updated 2020-09-05
+#' Updated 2021-10-22
 
 #' @param nsteps number of time steps to be evaluated
 
@@ -15,6 +15,8 @@
 #' @param seam3 communication adjacency matrix, read in if readseam3=T (rows as sources and columns as sinks)
 #' @param bpam3 dispersal adjacency matrix, read in if readbpam3=T (rows as sources and columns as sinks)
 #' @param geocoords3n matrix of x,y coordinates of nodes
+#' @param roaddistfilepath3n filepath to a R Data (.RData extension) file containing an adjacency matrix with the distances via the road network between nodes in geocoords - mutually exclusive with roadtimefilepath3n
+#' @param roadtimefilepath3n filepath to a R Data (.RData extension) file containing an adjacency matrix with the travel times via the road network between nodes in geocoords - mutually exclusive with roaddistfilepath3n
 
 #' @param seamdist3 com (genmovnet) the function of distance used to estimate movement probability - 'random' (not related to distance) or 'powerlaw' 
 #' @param seampla3 com (genmovnet) inverse power law parameter a in ad^(-b)
@@ -69,7 +71,7 @@
 
 
 
-ntsteps2 <- function(nsteps, geocoords3n, infon, vect1cn, vect1dn, readseam3, seam3, seamdist3, seampla3, seamplb3, seamrandp3, readbpam3, bpam3, bpamdist3, bpampla3, bpamplb3, bpamrandp3, readprobadoptvec3, probadoptvec3, probadoptmean3, probadoptsd3, maneffdir3, maneffmean3n, maneffsd3n, readprobestabvec3, probestabvec3, probestabmean3, probestabsd3, plotmpn=F){
+ntsteps2 <- function(nsteps, geocoords3n, roaddistfilepath3n, roadtimefilepath3n, infon, vect1cn, vect1dn, readseam3, seam3, seamdist3, seampla3, seamplb3, seamrandp3, readbpam3, bpam3, bpamdist3, bpampla3, bpamplb3, bpamrandp3, readprobadoptvec3, probadoptvec3, probadoptmean3, probadoptsd3, maneffdir3, maneffmean3n, maneffsd3n, readprobestabvec3, probestabvec3, probestabmean3, probestabsd3, plotmpn=F){
 
 # Numbering of steps refers to the numbering used in the INA User's Manual (1-5)
 
@@ -85,7 +87,7 @@ if(infon){
   # if generating the communication adj mat (not reading it in)
 
   if (!readseam3){ 
-    seam3 <- genmovnet(geocoords4n=geocoords3n, amdist4=seamdist3, amrandp4=seamrandp3, ampla4=seampla3, amplb4=seamplb3)
+    seam3 <- genmovnet(geocoords4n=geocoords3n, roaddistfilepath4n=roaddistfilepath3n, roadtimefilepath4n=roadtimefilepath3n, amdist4=seamdist3, amrandp4=seamrandp3, ampla4=seampla3, amplb4=seamplb3)
   }
 
   # 4.1. decision making - underlying probabilities of adoption
@@ -111,9 +113,10 @@ if(infon){
 # 3.2. bioentity dispersal adjacency matrix - underlying probabilities
 
 # if generating the bioentity dispersal network (rather than reading it in)
+# pathogens probably do not travel by road, so do not use road data
 
 if (!readbpam3){ 
-  bpam3 <- genmovnet(geocoords4n=geocoords3n, amdist4=bpamdist3, amrandp4=bpamrandp3, ampla4=bpampla3, amplb4=bpamplb3)  
+  bpam3 <- genmovnet(geocoords4n=geocoords3n, roaddistfilepath4n=roaddistfilepath3n, roadtimefilepath4n=roadtimefilepath3n, amdist4=bpamdist3, amrandp4=bpamrandp3, ampla4=bpampla3, amplb4=bpamplb3)
 }
 
 # 4.2. establishment - underlying probabilities of establishment
diff --git a/R/roaddata.R b/R/roaddata.R
new file mode 100644
index 0000000..1fccb4d
--- /dev/null
+++ b/R/roaddata.R
@@ -0,0 +1,219 @@
+#' Evaluates scenarios in an impact network analysis (INA)
+#'
+#' This function implements and summarizes multiple simulations across a designated range of parameter values
+#'
+#' Updated 2021-10-22
+
+#' @param geocoords5 matrix of x,y coordinates of nodes
+#' @param roaddistfilepath5n filepath to a R Data (.RData extension) file containing an adjacency matrix with the distances via the road network between nodes in geocoords - mutually exclusive with roadtimefilepath5n
+#' @param roadtimefilepath5n filepath to a R Data (.RData extension) file containing an adjacency matrix with the travel times via the road network between nodes in geocoords - mutually exclusive with roaddistfilepath5n
+#' @keywords road
+#' @export 
+#' @import osrm
+
+#' @examples
+
+#NEED TO WRITE EXAMPLES
+
+
+
+roaddata = function(geocoords5, roaddistfilepath5=NA, roadtimefilepath5=NA) {
+	filepath = "";
+	isDist = F;
+    if(!is.na(roaddistfilepath5)) {
+	filepath = roaddistfilepath5;
+	isDist = T;
+    } else if(!is.na(roadtimefilepath5))  {
+	filepath = roadtimefilepath5;
+	isDist = F;
+    } else {
+	stop("Please specify either roaddistfilepath or roadtimefilepath");
+    }
+
+	if(file.exists(filepath)) {
+		load(filepath);
+
+		if(isDist) {
+			return(distMat); # Loaded from file
+		} else {
+			return(timeMat); # Loaded from file
+		}
+
+	} else {
+
+		message("The data file containing the adjacency matrix with distances between locations has not been found.");
+		message("It will be downloaded and generated using the OSRM API.");
+		message("If you have this file, ensure that the correct filepath has been entered.");
+		message("This download may take a while.");
+
+    continue = readline(prompt="Are you sure you would like to continue (Y/n): ");
+
+    if(!(continue == "Y" || continue == "y")) {
+      stop("Either do not specify roaddatafilepath, specify a valid roaddatafilepath, or allow download to occur for simulation to proceed");
+    }
+
+		# Load osrm library if not already loaded
+		# Fail if osrm is not installed
+		allLibs = library()$results[,1];
+
+		if(!is.element("osrm", allLibs)) {
+			message("osrm library not installed, exiting");
+      message("Please install osrm package to use road distance functionality");
+			return;
+		}
+
+		loadedLibs = (.packages());
+
+		if(!is.element("osrm", loadedLibs)) {
+			library(osrm);
+		}
+
+		# With osrm loaded, proceed to get data
+
+
+    #BIG ASSUMPTION: geocoords (x,y) data coorresponds to coordinate (lon, lat) values
+    #Convert coordinates into list, needed for function
+    #Also need to add id for OSRM package
+    lon = geocoords5[,1]
+    lat = geocoords5[,2]
+    id = 1:length(lon)
+    #List varible order is important, OSRM breaks if out of order
+    locations = data.frame(id=id, lon=lon, lat=lat)
+
+	count = length(locations$id);
+	distMat = matrix(nrow=count, ncol=count, dimnames=list(locations$id,locations$id));
+	timeMat = matrix(nrow=count, ncol=count, dimnames=list(locations$id,locations$id));
+
+	#The OSRM API's default server does not allow the generation of a matrix with greater than or equal to
+	#10000 entries (a 100x100 square matrix, 50x200 rectangular matrix, etc.)
+	#Thus, for lists with 100 points or more, the matrix needs to be loaded in 99x99 submatrices, and stitched together
+	#into the final matrix
+
+	repCount = ceiling(count / 99);
+
+	loadedCount = 0;
+	for(y in 1:repCount) {
+		for(x in 1:repCount) {
+			#Declare lower and upper bounds for locations to be downloaded
+			xl = 1 + 99 * (x - 1);
+			yl = 1 + 99 * (y - 1);
+			xu = min(99 * x, count); #Don't exceed maximum matrix dimension
+			yu = min(99 * y, count); #Don't exceed maximum matrix dimension
+
+			#Get result from OSRM
+			result = osrmTable(
+				src=locations[yl:yu, c("id", "lon", "lat")],
+				dst=locations[xl:xu, c("id", "lon", "lat")],
+				measure=c("distance", "duration"), osrm.profile="car"
+			);
+
+			distMat[yl:yu,xl:xu] = result$distances #Copy from result into large matrix
+			timeMat[yl:yu,xl:xu] = result$durations #Copy from result into large matrix
+			loadedCount = loadedCount + (xu-xl+1) * (yu-yl+1); #Update progress counter
+
+			#Output progress
+			message(paste("Downloaded ", loadedCount, " of ", count^2, " distances, ", trunc(100 * loadedCount / count^2, digits=5), "% done", sep=""));
+		}
+	}
+
+    #Check for NA values, should only need to do for either distMat or timeMat, but not both (would be redundant)
+    temp = c(distMat);
+    #idx = 1:length(temp)
+    if(length(temp[is.na(temp)]) > 0) {
+      message("There are some NA entries in the adjacency matrix.");
+      message("This is because a road network route could not be found between at least one pair of locations.");
+      message("This may occur if, for instance, one location is on an island.")
+      message("These can either be replaced with Inf (infinity), these locations can be removed from the dataset.");
+      continue = readline(prompt="Substitute for infinity, or terminate and manually filter dataset (I/r): ");
+
+      if(!(continue == "I" || continue == "i")) {
+        stop("Please rerun with an approppriate dataset");
+      }
+
+      #Replace all NA with Inf
+      for(i in 1:count) {
+        for(j in 1:count) {
+          if(is.na(distMat[i,j])) {
+            distMat[i,j] = Inf;
+          }
+        }
+      }
+
+
+      for(i in 1:count) {
+        for(j in 1:count) {
+          if(is.na(distMat[i,j])) {
+            timeMat[i,j] = Inf;
+          }
+        }
+      }
+
+      #Probably more efficient way
+      #temp[is.na(temp)] = Inf;
+      #asdf=matrix(nrow=nrow(distMat), temp);
+
+    }
+
+		save( distMat, timeMat, file=filepath);
+
+		if(isDist) {
+			return(distMat); # Loaded from file
+		} else {
+			return(timeMat); # Loaded from file
+		}
+	}
+}
+
+validate = function(locations, distMat) {
+	x = round(runif(1, 1, 228));
+	y = round(runif(1, 1, 228));
+	#x=228;
+	#y=228;
+	width = round(runif(1, 1, 98));
+	height = round(runif(1, 1, 98));
+
+	print(x)
+	print(y)
+	print(width)
+	print(height)
+
+
+	#waste = file.remove("roaddist.RData");
+
+	#distMat = genDist(locations, "roaddist.RData");
+
+	download = osrmTable(
+			src=locations[y:(y+98), c("id", "lon", "lat")],
+			dst=locations[x:(x+98), c("id", "lon", "lat")],
+			measure="distance", osrm.profile="car"
+		);
+
+	download2 = osrmTable(
+			src=locations[y:(y+height - 1), c("id", "lon", "lat")],
+			dst=locations[x:(x+width - 1), c("id", "lon", "lat")],
+			measure="distance", osrm.profile="car"
+		);
+
+	#str(download$distances);
+	#str(distMat[y:(y+98),x:(x+98)]);
+
+	#print(download$distances - distMat[y:(y+98),x:(x+98)]);
+
+	#print(sum(abs(distMat[y:(y+98),x:(x+98)] - download$distances)));
+	#print(sum((distMat[y:(y+98),x:(x+98)] - download$distances)));
+
+	print(sum(abs(download$distances[1:(height),1:(width)] - download2$distances)));
+	print(sum((download$distances[1:(height),1:(width)] - download2$distances)));
+
+	#print(download$distances);
+	#print(download2$distances);
+
+	write.table(download$distances, file="download.txt");
+	write.table(download2$distances, file="download2.txt");
+}
+
+#test = function() {
+	#cities = read.csv("cities.csv");
+	#filepath = "roaddist.RData";
+	#waste = genDist(cities,filepath);
+#}