Constructing Urban Areas (part 2)

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)

# Select cities 
cities <- c(101,751,461,851,561)

# 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_oc:=sum(n),by=.(origin)]
cmu[,N_dc:=sum(n),by=.(destination)]


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)
        # (2) For the zone compute importance of other areas
        # Calculate how many work in zone
        cmu[,N_d_zone:=sum(n*destination%in%zone)]
        # Calculate how many live in zone
        cmu[,N_o_zone:=sum(n*origin%in%zone)]
        # Calculate the importance of the area - zone link
        # For all areas calculate the Pr(work in zone|live in j)
        cmu[,N_dZCo:=sum(n*destination%in%zone)/N_oc,by=.(origin)]
        # For all areas calculate the Pr(live in area|work in zone)
        cmu[,N_oZCd:=sum(n*destination%in%zone)/N_d_zone,by=.(origin)]

        temp <- cmu[,.(d=N_dZCo[1],s=N_oZCd[1],m=round(N_dZCo[1]*N_oZCd[1],3)),.(origin)]
        temp <- temp[order(m),]
        


        # Test if zone is self-containing/importance of zone - zone link
        # Probability of working in zone given that you live in zone
        cmu[,N_odz:=sum(n*destination%in%zone*origin%in%zone)]
        a1 <- cmu$N_odz[1]/cmu$N_d_zone[1]
        # Probability of living in zone given that you work in zone
        a2 <- cmu$N_odz[1]/cmu$N_o_zone[1]
        
        test <- min(a1,a2)

        add_to_zone <- rev(temp[!origin%in%zone,]$origin)[1]        
        }

    com_zones[match(zone,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)