mypath <- "C://Users//np83zg//OneDrive - Aalborg Universitet//Skrivebord//denmark"
setwd(mypath)
library(data.table)
mydata1 <- as.data.table(read.table("dat1.txt"))
mydata2 <- as.data.table(read.table("dat2.txt"))
mydata3 <- as.data.table(read.table("dat3.txt"))
mydata4 <- as.data.table(read.table("dat4.txt"))


setkey(mydata1,workplace,year)
mydata1[,density:=count/size]
mydata_2003 <- mydata1[year==2003,]
mydata_2008 <- mydata1[year==2008,]
mydata_2016 <- mydata1[year==2016,]

mydata_2016[,comsize:=sum(size),by=.(commutingzone)]
mydata_2016[,comcount:=sum(count),by=.(commutingzone)]
mydata_2016[,com_density:=comcount/comsize]
mydata_2016[,com_log_density:=log(comcount/comsize)]
mydata_2016[,com_m_wage:=sum(size*m_wage)/comsize,by=.(commutingzone)]
mydata_2016[,com_s_uni_edu:=sum(size*s_uni_edu)/comsize,by=.(commutingzone)]
mydata_2016[,com_m_jobtenure:=sum(size*m_jobtenure)/comsize,by=.(commutingzone)]
mydata_2016[,com_m_exper:=sum(size*m_exper)/comsize,by=.(commutingzone)]

index <- match(unique(mydata_2016$commutingzone), mydata_2016$commutingzone)
index
##  [1]  1 42 44 46 47 48 55 56 58 59 60 61 70 71 72 74 80 84 86 87 88 90 91 93 98
comtab <- mydata_2016[index,]
comtab <- comtab[order(comtab$com_density),]
comtab <- comtab[,c("commutingzone","workplace","comcount","comsize","com_density","com_m_wage",
                         "com_m_exper","com_s_uni_edu","com_m_jobtenure")]


# Load municipality shapes
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tmap)
mun_shape <- st_read("mun_shape.shp")
## Reading layer `mun_shape' from data source `C:\Users\np83zg\OneDrive - Aalborg Universitet\Skrivebord\denmark\mun_shape.shp' using driver `ESRI Shapefile'
## Simple feature collection with 98 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 8.074458 ymin: 54.55906 xmax: 15.15735 ymax: 57.75238
## geographic CRS: WGS 84
tmap_mode("plot")
## tmap mode set to plotting
myzone_1 <- c(101,147,175,157,185,159,167,173,161,153)
myzone_2 <- c(163,187,151,183,165,155,190,169,253,223)
myzone_3 <- c(201,269,265,219,240,217,259,260,230,250)
mun_shape$zone <- 1*as.numeric(mun_shape$kode%in%myzone_1) + 
    2*as.numeric(mun_shape$kode%in%myzone_2) + 
    3*as.numeric(mun_shape$kode%in%myzone_3)

tm_shape(mun_shape) + tm_fill("seashell3") + tm_borders("white", lwd = 1.2)

tm_shape(mun_shape) + tm_polygons("zone")