Input:
#Call packages:
pacman::p_load(sf,
rio,
here,
janitor,
tidyverse,
dplyr,
magrittr,
ggplot2,
purrr,
lubridate,
mice,
plotly)
#Import Viet Nam shapefile in level 1 (include 63 provinces):
getwd()
## [1] "C:/Users/locca/Documents/Xuân Lộc/R/Project R"
#Take new workdirectory in clipboard
vn_adm1<-sf::st_read(r"(C:\Users\locca\Documents\Xuân Lộc\R\Data R\Viet Nam shapefile\vnm_admbnda_adm1_gov_20201027.shp)")
## Reading layer `vnm_admbnda_adm1_gov_20201027' from data source
## `C:\Users\locca\Documents\Xuân Lộc\R\Data R\Viet Nam shapefile\vnm_admbnda_adm1_gov_20201027.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 63 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 102.144 ymin: 7.180931 xmax: 117.8355 ymax: 23.39221
## Geodetic CRS: WGS 84
optimize<-read.csv(r"(C:\\Users\\locca\\Documents\\Xuân Lộc\\VILAS\\Final project\\Optimize_df.csv)")
VNnorth<-vn_adm1 %>%
filter(str_detect(str_c(c("Lào Cai,Yên Bái,Điện Biên,Hòa Bình,Lai Châu, Sơn La,Hà Giang,Cao Bằng,Bắc Kạn,Lạng Sơn,Tuyên Quang,Thái Nguyên,Phú Thọ, Bắc Giang,Quảng Ninh,Bắc Ninh,Hà Nam,TP. Hà Nội,Hải Dương,Hưng Yên,Hải Phòng,Nam Định,Ninh Bình,Thái Bình,Vĩnh Phúc,Thanh Hóa,Nghệ An"),"|"),vn_adm1$ADM1_VI)
)
#Add facility positions
manufacter<-data.frame(Lat =c(20.53046974118655),
Long = c(105.90021448808952))
manufacter_sf <- manufacter %>%
st_as_sf(coords = c('Long', 'Lat')) %>%
st_set_crs("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
#Calculate distances:
points<- optimize[,c("Longitude","Latitude")]
library(geosphere)
points<- points %>%
mutate(Dist = distHaversine(manufacter[,c("Long","Lat")],points)/1000)
df<-left_join(points,optimize,by = c("Longitude","Latitude"))
reg1<-lm(Average_Leadtime~Dist,
data =df)
##So we will have y = 0.68404 + 0.005105x refer the correlation between leadtime and distance.
#Create buffer:
manufacter_1day <- st_buffer(manufacter_sf,
units::as_units(61, 'kilometer'))
manufacter_1.5day <- st_buffer(manufacter_sf,
units::as_units(160, 'kilometer'))
manufacter_2day <- st_buffer(manufacter_sf,
units::as_units(257, 'kilometer'))
#New manufacter:
new_manufacter= data.frame(Lat =c(21.12256201),
Long = c(105.9150683))
Output:
#Plot north region map:
library(ggplot2)
library(viridis)
## Loading required package: viridisLite
p<-VNnorth %>%
ggplot()+
##Plot north region map:
geom_sf(data = VNnorth,
aes(geometry = geometry))+
##Highlight facility postion and it's buffer value:
geom_sf(data = manufacter_1day, fill = 'blue', alpha = 0.3) +
geom_sf(data = manufacter_1.5day, fill = 'blue', alpha = 0.2) +
geom_sf(data = manufacter_2day, fill = 'blue', alpha = 0.1) +
geom_sf(data = manufacter_sf, color = 'red')+
##Locate the destination of customers:
geom_point(data = df,
aes(x = Longitude,
y = Latitude,
size = Total.transactions,
color= Average_Leadtime))+
theme_void() +
coord_sf() +
theme(
legend.position = c(0.1, 0.25),
text = element_text(color = '#22211d'),
plot.background = element_rect(fill = '#f5f5f2', color = NA),
panel.background = element_rect(fill = '#f5f5f2', color = NA),
legend.background = element_rect(fill = '#f5f5f2', color = NA),
plot.title = element_text(size= 16, hjust=0.1, color = '#2e2d27',
margin = margin(b = -0.1, t = 0.4, l = 2,
unit = 'cm'),
face = 'bold'),
)+
ggtitle("Mean shipping time")
p <- ggplotly(p)
p
#Correlations:
library(corrplot)
## corrplot 0.92 loaded
M<-cor(df[,c("Average_Leadtime",
"Dist",
"Total.transactions",
"Sum.of.Total.Weight")])
corrplot(M,
method="circle")
#Modeling:
reg<-lm(Average_Leadtime~Dist+Total.transactions+Sum.of.Total.Weight,
data = df)
summary(reg)
##
## Call:
## lm(formula = Average_Leadtime ~ Dist + Total.transactions + Sum.of.Total.Weight,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.67851 -0.21949 -0.07450 0.05458 1.88984
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.424e+00 2.209e-01 10.973 1.96e-14 ***
## Dist -4.532e-04 1.088e-03 -0.417 0.679
## Total.transactions -6.331e-03 6.258e-04 -10.116 2.81e-13 ***
## Sum.of.Total.Weight 1.075e-06 1.627e-07 6.610 3.52e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4903 on 46 degrees of freedom
## Multiple R-squared: 0.7426, Adjusted R-squared: 0.7258
## F-statistic: 44.24 on 3 and 46 DF, p-value: 1.325e-13
library(osrm)
## Data: (c) OpenStreetMap contributors, ODbL 1.0 - http://www.openstreetmap.org/copyright
## Routing: OSRM - http://project-osrm.org/
#Create table:
locs<-data.frame(c("Manufacter"),new_manufacter[,c("Long","Lat")],0,0)
colnames(locs)<-c("Customers",
"Longitude",
"Latitude",
"Sum.of.Total.Weight",
"Average_Leadtime")
locs<-rbind(locs,
optimize %>% select(Customers,
Longitude,
Latitude,
Sum.of.Total.Weight,
Average_Leadtime))
locs<-locs %>%
mutate(Sum.of.Total.Weight = round(Sum.of.Total.Weight,2),
Average_Leadtime = round(Average_Leadtime,2))
#Create labels:
labels<- paste0("<strong> Customers </strong> ",
locs$Customers, "<br/> ",
"<strong> Total weight: </strong> ",
locs$Longitude, "<br/> ",
"<strong> Total weight: </strong> ",
locs$Latitude, "<br/> ",
"<strong> Distance: </strong> ",
locs$Sum.of.Total.Weight, "<br/> ",
"<strong> Ship time </strong> ",
locs$Average_Leadtime, "<br/> ") %>%
lapply(htmltools::HTML)
# Calculate the shortest round trip between the facility locations
trip <- osrmTrip(locs[,c(2:3)],
overview = "full")
# Get the spatial lines dataframe
trip_sp <- trip[[1]]$trip
#Plot map with leaflet:
library(leaflet)
leaflet(data = trip_sp) %>%
addTiles() %>%
addMarkers(lng = locs$Longitude,
lat = locs$Latitude,
popup = locs$Customers) %>%
addAwesomeMarkers(lng = locs$Longitude,
lat = locs$Latitude,
label = ~ labels) %>%
addPolylines()
#Overall:
##Summary output:
library(DT)
datatable(
st_drop_geometry(route <- trip_sp %>%
mutate(duration = round(duration, 1),
distance = round(distance, 1))))
#Total:
trip_summary <- trip[[1]]$summary
trip_summary
## $duration
## [1] 2335.92
##
## $distance
## [1] 2852.003
#Visualize the output
trip2 <- osrmIsodistance(new_manufacter[,c("Long","Lat")],
breaks = c(0,
61*10^3,
160*10^3,
257*10^3,
450*10^3))
trip2$shiptime<-factor(paste(trip2$isomin, "to",trip2$isomax),
levels = c("0 to 61000",
"61000 to 160000",
"160000 to 257000",
"257000 to 450000"),
label = c("0 to 1 day",
"1 to 1.5 day",
"1.5 to 2 day",
"2 to 3 day")
)
pal <- colorFactor(rev(heat.colors(4)), trip2$shiptime)
labels<- paste0("<strong> Customers </strong> ",
locs$Customers, "<br/> ",
"<strong> Total weight: </strong> ",
locs$Longitude, "<br/> ",
"<strong> Total weight: </strong> ",
locs$Latitude, "<br/> ",
"<strong> Distance: </strong> ",
locs$Sum.of.Total.Weight, "<br/> ",
"<strong> Ship time </strong> ",
locs$Average_Leadtime, "<br/> ") %>%
lapply(htmltools::HTML)
#Plot map with leaflet:
library(leaflet)
leaflet(data = trip2) %>%
addTiles() %>%
addMarkers(lng = locs$Longitude,
lat = locs$Latitude,
label = ~labels) %>%
addPolygons(fill = TRUE,
stroke = TRUE,
color = "black",
fillColor = ~pal(trip2$shiptime),
weight = 0.5,
fillOpacity = 0.3,
data = trip2,
popup = trip2$shiptime,
group = "Running Time") %>%
addLegend("bottomright",
pal = pal,
values = trip2$shiptime,
title = "Shipping Time")