The algorithm constructs urban areas using a Travel To Work Area approach. The algorithm starts with a number of city centers and grows commuting zones around these centers. The growth is done by adding the area \(a\) with the strongest commute link \(L(z,a)\) to the zone \(z\) that is being grown. Algorithm breaks when \(Pr(live \in zone|work \in zone)>\alpha\) and \(Pr(work\ \in\ zone|live\ \in\ zone)>\alpha\).

In the first set of commuting zones \(\alpha=.7\), in the second \(\alpha=.8\) and in the last it is \(\alpha=.85\).

rm(list=ls())
path <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//transportzones//"
setwd(path)
library(data.table)
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(ttwa)
## Warning: package 'ttwa' was built under R version 4.0.3
library(tmap)


# First load commute flows
cmu <- read.table("cmu.txt",header=TRUE)
cmu <- setDT(cmu)


# Load a shape file for municipalities
shp <- read_sf("mun_shape.shp")
shp <- shp[order(shp$kode),]

# Find all neighbours
stn <- st_intersects(shp)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
adj <- data.frame(origin=rep(0,98*98),destination=rep(0,98*98),neighbour=rep(0,98*98))
bool <- c()
org <- c()
des <- c()
for (i in 1:98)
    {
        org <- c(org,rep(shp$kode[i],98))
        des <- c(des,shp$kode)
        add <- rep(0,98)
        add[unlist(stn[i])] <- 1
        bool <- c(bool,add)
    }
adj$neighbour <- as.logical(bool)
adj$destination <- des
adj$origin <- org
adj <- setDT(adj)
cmu <- merge(adj,cmu,by.x=c("origin","destination"),by.y=c("residence","workplace"),all.x=TRUE)
cmu[,n:=ifelse(is.na(nn),0,nn)]
cmu[,nn:=NULL]
cmu[,N_o:=sum(n),by=.(origin)]
cmu[,N_d:=sum(n),by=.(destination)]




cities <- c(101,751,461,851,561)
com_zones <- rep(0,98)
for (i in 1:length(cities))
{
    # Growing the zone ... using while loop
    # (1) First choose the center (and set zone to center)
    center <- cities[i]
    zone <- center

    # (2) Define cut_off (share commuting to center)
    cut_off <- 0.70
    test <- 0
    add_to_zone <- c()
    while( test < cut_off)
        {
        
        zone <- c(zone,add_to_zone)

        # Make dummies
        cmu[,I_Lz:=origin%in%zone]
        cmu[,I_Wz:=destination%in%zone]
        

        # Calculate how many live in zone
        cmu[,nL_z:=sum(n*I_Lz)]
        # Calculate how many work in zone
        cmu[,nW_z:=sum(n*I_Wz)]

        # Calculate Pr(work in zone|live in o)
        cmu[,n_wz_lo:=sum(n*I_Wz),by=.(origin)]
        cmu[,p1:=n_wz_lo/N_o]

        # Calculate Pr(live in o|work in z)
        cmu[,p2:=n_wz_lo/nW_z]

        # Calculate Pr(work in o|live in z)
        cmu[,n_wo_lz:=sum(n*I_Lz),by=.(destination)]
        cmu[,p3:=n_wo_lz/nL_z]

        # Calculate Pr(live in z|work in o)
        cmu[,p4:=n_wo_lz/N_d]

        # Calculate weights
        cmu[,w1:=round(p1*p2,8)]
        cmu[,w2:=round(p3*p4,8)]
        cmu[,w:=w1+w2]
        temp <- cmu[origin==destination,.(origin,p1,p2,p3,p4,w)]
        setkey(temp,w)
        temp <- temp[!origin%in%zone,]
        
        cmu[,n_odz:=sum(n*origin%in%zone*destination%in%zone)]
        a1 <- cmu$n_odz[1]/cmu$nL_z[1]
        a2 <- cmu$n_odz[1]/cmu$nW_z[1]

        test <- min(a1,a2)

        add_to_zone <- rev(temp$origin)[1]      
        }

    com_zones[match(zone,sort(unique(cmu$origin)))] <- i

}
shp$zone <- as.factor(com_zones)
temp_shp <- shp #[shp$kode<400,]

## tmap mode set to plotting
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(temp_shp) + tm_polygons("zone",palette="Accent") +
tm_shape(temp_shp) + tm_borders("white", lwd = 1.2)

# First load commute flows
cmu <- read.table("cmu.txt",header=TRUE)
cmu <- setDT(cmu)


# Load a shape file for municipalities
shp <- read_sf("mun_shape.shp")
shp <- shp[order(shp$kode),]

# Find all neighbours
stn <- st_intersects(shp)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
adj <- data.frame(origin=rep(0,98*98),destination=rep(0,98*98),neighbour=rep(0,98*98))
bool <- c()
org <- c()
des <- c()
for (i in 1:98)
    {
        org <- c(org,rep(shp$kode[i],98))
        des <- c(des,shp$kode)
        add <- rep(0,98)
        add[unlist(stn[i])] <- 1
        bool <- c(bool,add)
    }
adj$neighbour <- as.logical(bool)
adj$destination <- des
adj$origin <- org
adj <- setDT(adj)
cmu <- merge(adj,cmu,by.x=c("origin","destination"),by.y=c("residence","workplace"),all.x=TRUE)
cmu[,n:=ifelse(is.na(nn),0,nn)]
cmu[,nn:=NULL]
cmu[,N_o:=sum(n),by=.(origin)]
cmu[,N_d:=sum(n),by=.(destination)]




cities <- c(101,751,461,851,561)
com_zones <- rep(0,98)
for (i in 1:length(cities))
{
    # Growing the zone ... using while loop
    # (1) First choose the center (and set zone to center)
    center <- cities[i]
    zone <- center

    # (2) Define cut_off (share commuting to center)
    cut_off <- 0.80
    test <- 0
    add_to_zone <- c()
    while( test < cut_off)
        {
        
        zone <- c(zone,add_to_zone)

        # Make dummies
        cmu[,I_Lz:=origin%in%zone]
        cmu[,I_Wz:=destination%in%zone]
        

        # Calculate how many live in zone
        cmu[,nL_z:=sum(n*I_Lz)]
        # Calculate how many work in zone
        cmu[,nW_z:=sum(n*I_Wz)]

        # Calculate Pr(work in zone|live in o)
        cmu[,n_wz_lo:=sum(n*I_Wz),by=.(origin)]
        cmu[,p1:=n_wz_lo/N_o]

        # Calculate Pr(live in o|work in z)
        cmu[,p2:=n_wz_lo/nW_z]

        # Calculate Pr(work in o|live in z)
        cmu[,n_wo_lz:=sum(n*I_Lz),by=.(destination)]
        cmu[,p3:=n_wo_lz/nL_z]

        # Calculate Pr(live in z|work in o)
        cmu[,p4:=n_wo_lz/N_d]

        # Calculate weights
        cmu[,w1:=round(p1*p2,8)]
        cmu[,w2:=round(p3*p4,8)]
        cmu[,w:=w1+w2]
        temp <- cmu[origin==destination,.(origin,p1,p2,p3,p4,w)]
        setkey(temp,w)
        temp <- temp[!origin%in%zone,]
        
        cmu[,n_odz:=sum(n*origin%in%zone*destination%in%zone)]
        a1 <- cmu$n_odz[1]/cmu$nL_z[1]
        a2 <- cmu$n_odz[1]/cmu$nW_z[1]

        test <- min(a1,a2)

        add_to_zone <- rev(temp$origin)[1]      
        }

    com_zones[match(zone,sort(unique(cmu$origin)))] <- i

}
shp$zone <- as.factor(com_zones)
temp_shp <- shp #[shp$kode<400,]

## tmap mode set to plotting
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(temp_shp) + tm_polygons("zone",palette="Accent") +
tm_shape(temp_shp) + tm_borders("white", lwd = 1.2)

# First load commute flows
cmu <- read.table("cmu.txt",header=TRUE)
cmu <- setDT(cmu)


# Load a shape file for municipalities
shp <- read_sf("mun_shape.shp")
shp <- shp[order(shp$kode),]

# Find all neighbours
stn <- st_intersects(shp)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
adj <- data.frame(origin=rep(0,98*98),destination=rep(0,98*98),neighbour=rep(0,98*98))
bool <- c()
org <- c()
des <- c()
for (i in 1:98)
    {
        org <- c(org,rep(shp$kode[i],98))
        des <- c(des,shp$kode)
        add <- rep(0,98)
        add[unlist(stn[i])] <- 1
        bool <- c(bool,add)
    }
adj$neighbour <- as.logical(bool)
adj$destination <- des
adj$origin <- org
adj <- setDT(adj)
cmu <- merge(adj,cmu,by.x=c("origin","destination"),by.y=c("residence","workplace"),all.x=TRUE)
cmu[,n:=ifelse(is.na(nn),0,nn)]
cmu[,nn:=NULL]
cmu[,N_o:=sum(n),by=.(origin)]
cmu[,N_d:=sum(n),by=.(destination)]




cities <- c(101,751,461,851,561)
com_zones <- rep(0,98)
for (i in 1:length(cities))
{
    # Growing the zone ... using while loop
    # (1) First choose the center (and set zone to center)
    center <- cities[i]
    zone <- center

    # (2) Define cut_off (share commuting to center)
    cut_off <- 0.85
    test <- 0
    add_to_zone <- c()
    while( test < cut_off)
        {
        
        zone <- c(zone,add_to_zone)

        # Make dummies
        cmu[,I_Lz:=origin%in%zone]
        cmu[,I_Wz:=destination%in%zone]
        

        # Calculate how many live in zone
        cmu[,nL_z:=sum(n*I_Lz)]
        # Calculate how many work in zone
        cmu[,nW_z:=sum(n*I_Wz)]

        # Calculate Pr(work in zone|live in o)
        cmu[,n_wz_lo:=sum(n*I_Wz),by=.(origin)]
        cmu[,p1:=n_wz_lo/N_o]

        # Calculate Pr(live in o|work in z)
        cmu[,p2:=n_wz_lo/nW_z]

        # Calculate Pr(work in o|live in z)
        cmu[,n_wo_lz:=sum(n*I_Lz),by=.(destination)]
        cmu[,p3:=n_wo_lz/nL_z]

        # Calculate Pr(live in z|work in o)
        cmu[,p4:=n_wo_lz/N_d]

        # Calculate weights
        cmu[,w1:=round(p1*p2,8)]
        cmu[,w2:=round(p3*p4,8)]
        cmu[,w:=w1+w2]
        temp <- cmu[origin==destination,.(origin,p1,p2,p3,p4,w)]
        setkey(temp,w)
        temp <- temp[!origin%in%zone,]
        
        cmu[,n_odz:=sum(n*origin%in%zone*destination%in%zone)]
        a1 <- cmu$n_odz[1]/cmu$nL_z[1]
        a2 <- cmu$n_odz[1]/cmu$nW_z[1]

        test <- min(a1,a2)

        add_to_zone <- rev(temp$origin)[1]      
        }

    com_zones[match(zone,sort(unique(cmu$origin)))] <- i

}
shp$zone <- as.factor(com_zones)
temp_shp <- shp #[shp$kode<400,]

## tmap mode set to plotting
tmap_mode("plot")
## tmap mode set to plotting
tm_shape(temp_shp) + tm_polygons("zone",palette="Accent") +
tm_shape(temp_shp) + tm_borders("white", lwd = 1.2)

```