Compositional spline of tmax- all province in 2015

In this report, I have used all OUR MODIFIED functions which are based on robCompositions package

1. Import climate data

## ===== Packages =================================

library(readstata13)        ## Importing Stata data
library(ggplot2)            ## Figures

## ==== Other packages ============================
## Source = https://swilke-geoscience.net/post/spatial_interpolation/


# Some packages for (spatial) data processing
library(tidyverse) # wrangling tabular data and plotting
library(sf) # processing spatial vector data
library(sp) # another vector data package necessary for continuity
library(raster) # processing spatial raster data. !!!overwrites dplyr::select!!!

# Different packages to test their interpolation functions
library(gstat)  # inverse distance weighted, Kriging
library(fields) # Thin Plate Spline
library(interp) # Triangulation
library(mgcv)   # Spatial GAM
library(automap)# Automatic approach to Kriging

# Finally, some packages to make pretty plots
library(patchwork)
library(viridis)
require(readxl)
require(stringr)
##--------------------Transform hist using clr
require(robCompositions)
require(splines)
## ===== Load data ================================

setwd("D:/Compositional regression/Climate")
tmax_2015 <-read.dta13("tmax_2015.dta")

tentinh <- read_excel("Tentinh.xlsx")

source("FunctionCRF.R") # function fcenLRF
## ===== Data manipulations =======================
# Change name to capita


tmax_2015 <- tmax_2015 %>% mutate(province = str_to_upper(province))
setdiff(tmax_2015$province, tentinh$PROVINCE_NAME2)

##  [1] "YUNNAN"             "GUANGXI ZHUANG"     "GUANGXI"           
##  [4] "PHONGSALY"          "HANOI"              "STATE NOT FOUND"   
##  [7] "OUDOMXAY"           "LUANG PRABANG"      "HUAPHANH"          
## [10] "HAIPHONG"           "XIENGKHUANG"        "HAINAN"            
## [13] "VIENTIANE"          "XAYSOMBOON"         "BORIKHAMXAY"       
## [16] "VIENTIANE CAPITAL"  "NONG KHAI"          "BUENG KAN PROVINCE"
## [19] "UDON THANI"         "SAKON NAKHON"       "NAKHON PHANOM"     
## [22] "KHAMMUANE"          "NONG BUA LAMPHU"    "KHON KAEN"         
## [25] "KALASIN"            "MUKDAHAN"           "SAVANNAKHET"       
## [28] "CHAIYAPHUM"         "MAHA SARAKHAM"      "ROI ET"            
## [31] "SARAVANE"           "NAKHON RATCHASIMA"  "YASOTHON"          
## [34] "AMNAT CHAROEN"      "UBON RATCHATHANI"   "SEKONG"            
## [37] "BURIRAM"            "SURIN"              "SISAKET"           
## [40] "CHAMPASACK"         "ATTAPEU"            "ODDAR MEANCHEY"    
## [43] "PREAH VIHEAR"       "STUNG TRENG"        "RATANAKIRI"        
## [46] "SA KAEO"            "BANTEAY MEANCHEY"   "SIEM REAP"         
## [49] "CHANTHABURI"        "BATTAMBANG"         "KAMPONG THOM"      
## [52] "KRATIE"             "MONDULKIRI"         "PAILIN"            
## [55] "POUTHISAT"          "TRAT PROVINCE"      "KAMPONG CHHNANG"   
## [58] "KAMPONG CHAM"       "TBOUNG KHMUM"       "KOH KONG"          
## [61] "KAMPONG SPEU"       "KANDAL"             "PREY VENG"         
## [64] "SVAY RIENG"         "PREAH SIHANOUK"     "KAMPOT"            
## [67] "TAKEO"              "HO CHI MINH CITY"   "BA RIA-VUNG TAU"

setdiff(tentinh$PROVINCE_NAME2, tmax_2015$province)

##  [1] "HA NOI"            "VINH PHUC"         "BAC NINH"         
##  [4] "HAI PHONG"         "HUNG YEN"          "THAI BINH"        
##  [7] "HA NAM"            "DA NANG"           "BA RIA - VUNG TAU"
## [10] "HO CHI MINH"       "VINH LONG"         "CAN THO"          
## [13] "SOC TRANG"

tmax_2015 <- tmax_2015 %>% mutate(province =  case_when(province== "HANOI"  ~ "HA NOI" ,
                                                        province== "HANOI"  ~ "HA NOI",
                                                        province== "HAIPHONG"  ~ "HAI PHONG",
                                                        province== "BA RIA-VUNG TAU" ~  "BA RIA - VUNG TAU",
                                                        province== "HO CHI MINH CITY" ~ "HO CHI MINH" ,
                                                        province== "HO CHI MINH CITY" ~ "HO CHI MINH" ,
                                                        TRUE ~  province))
                                  
# No VINH PHUC, BAC NINH,   "HUNG YEN" ,  "THAI BINH", "HA NAM", "DA NANG",  "VINH LONG", "CAN THO" ,  "SOC TRANG"   

  
  
#filter only Vietnamese provinces

tmax_2015 <- tmax_2015 %>% filter(province %in% tentinh$PROVINCE_NAME2)
unique(tmax_2015$province)

##  [1] "HA GIANG"          "LAO CAI"           "CAO BANG"         
##  [4] "DIEN BIEN"         "LAI CHAU"          "TUYEN QUANG"      
##  [7] "BAC KAN"           "LANG SON"          "SON LA"           
## [10] "YEN BAI"           "THAI NGUYEN"       "PHU THO"          
## [13] "HA NOI"            "BAC GIANG"         "QUANG NINH"       
## [16] "HOA BINH"          "HAI DUONG"         "HAI PHONG"        
## [19] "THANH HOA"         "NINH BINH"         "NAM DINH"         
## [22] "NGHE AN"           "HA TINH"           "QUANG BINH"       
## [25] "QUANG TRI"         "THUA THIEN HUE"    "QUANG NAM"        
## [28] "QUANG NGAI"        "KON TUM"           "GIA LAI"          
## [31] "BINH DINH"         "DAK LAK"           "PHU YEN"          
## [34] "DAK NONG"          "KHANH HOA"         "TAY NINH"         
## [37] "BINH PHUOC"        "LAM DONG"          "NINH THUAN"       
## [40] "BINH DUONG"        "DONG NAI"          "BINH THUAN"       
## [43] "AN GIANG"          "LONG AN"           "HO CHI MINH"      
## [46] "BA RIA - VUNG TAU" "KIEN GIANG"        "DONG THAP"        
## [49] "BEN TRE"           "TIEN GIANG"        "HAU GIANG"        
## [52] "TRA VINH"          "CA MAU"            "BAC LIEU"

dim(tmax_2015)

## [1] 119 369

names(tmax_2015)[1] <- "longitude" 

2. Store all tmax for each province in list Ltmax2015

#---------Work with 1 provicine
require(reshape2)

unique(tmax_2015$province)

##  [1] "HA GIANG"          "LAO CAI"           "CAO BANG"         
##  [4] "DIEN BIEN"         "LAI CHAU"          "TUYEN QUANG"      
##  [7] "BAC KAN"           "LANG SON"          "SON LA"           
## [10] "YEN BAI"           "THAI NGUYEN"       "PHU THO"          
## [13] "HA NOI"            "BAC GIANG"         "QUANG NINH"       
## [16] "HOA BINH"          "HAI DUONG"         "HAI PHONG"        
## [19] "THANH HOA"         "NINH BINH"         "NAM DINH"         
## [22] "NGHE AN"           "HA TINH"           "QUANG BINH"       
## [25] "QUANG TRI"         "THUA THIEN HUE"    "QUANG NAM"        
## [28] "QUANG NGAI"        "KON TUM"           "GIA LAI"          
## [31] "BINH DINH"         "DAK LAK"           "PHU YEN"          
## [34] "DAK NONG"          "KHANH HOA"         "TAY NINH"         
## [37] "BINH PHUOC"        "LAM DONG"          "NINH THUAN"       
## [40] "BINH DUONG"        "DONG NAI"          "BINH THUAN"       
## [43] "AN GIANG"          "LONG AN"           "HO CHI MINH"      
## [46] "BA RIA - VUNG TAU" "KIEN GIANG"        "DONG THAP"        
## [49] "BEN TRE"           "TIEN GIANG"        "HAU GIANG"        
## [52] "TRA VINH"          "CA MAU"            "BAC LIEU"

#Now, store tmax in each province in a list, List tmax: Ltmax

funtmax <- function(PRO)
{
  TempPRO <- tmax_2015 %>% filter(province == PRO)
  
  TempPRO  <- reshape2::melt( TempPRO , 
                    id.vars = c("longitude",  "latitude" ,
                                "district",   "province"))
  TempPRO  <-  TempPRO  %>% rename(tmax = value)
  TempPRO  <-  TempPRO  %>% group_by(variable) %>%
    summarise(tmax = max(tmax))
  return( TempPRO )
}

#test funtmax for Hanoi

HN1 <- funtmax ("HA NOI")

#------APPLY funtmax()

Ltmax2015 <- list()
for (P in unique(tmax_2015$province))
{
  Ltmax2015[[P]] <- funtmax (P)
}

3. Store all compositional spline output in list Ltmax2015CS

#---------------Now, applying compositionalSpline() function

order = 4
der = 2
alpha = 0.3


funCS <- function(PRO)
{
  histempt <- hist( Ltmax2015[[PRO]]$tmax,
                 #xlab = "Tmax",
                # main = "Estimated histogram of tmax - Hanoi", 
                 #col = "blue", lwd = 2, 
                 breaks = 30,
                 plot = FALSE) # It does not work with 50 bins
  #t: histHN$breaks
  #density: histHN$density
  
  t_step = diff(histempt$breaks[1:2])
  
  
  #-----Using fcenLRF() by Huong
  
  knots <-  histempt$breaks[-length( histempt$breaks)]+0.5
  clrhisttempt <- fcenLRF(knots,   histempt$density ) #knots and  estimate
  
  #----using function in robCompositions
  t <- NULL
  clrf <- NULL
  knots<- NULL
  
  clrf <-  clrhisttempt$estimate
  t  <- clrhisttempt$knots
  knots <-  seq(min(Ltmax2015[[PRO]]$tmax), 
                max(Ltmax2015[[PRO]]$tmax), length = 10)
  w <- rep( 1/length(t), length(t))
  
  return(list(clrf=clrf, t = t, knots = knots, w = w))
  #compositionalSpline(t, clrf, knots, w, order,
   #                   der, alpha, 
    #                  spline.plot = TRUE, 
     #                 basis.plot = FALSE) #DONE
}

#---Now, apply for all provinces

Ltmax2015CS <- list()
for (P in unique(tmax_2015$province))
{
  Ltmax2015CS[[P]] <- funCS (P)
}

4. Now, plot all the spline of provinces, using compositionalSplineF

Plot for one province

P <- "DONG NAI"
order = 4
der = 2
alpha = 0.5
compositionalSplineF(Ltmax2015CS[[P]]$t, 
                    Ltmax2015CS[[P]]$clrf,
                    Ltmax2015CS[[P]]$knots,
                    Ltmax2015CS[[P]]$w, order,
                    der, alpha, 
                    spline.plot = TRUE, 
                    basis.plot = FALSE,
                    namemain = P) #DONE

## $J
##           [,1]
## [1,] 0.4575099
## 
## $ZB_coef
##              [,1]
##  [1,] -0.70757681
##  [2,] -1.97903576
##  [3,] -3.45542347
##  [4,] -4.53624184
##  [5,] -4.68004850
##  [6,] -4.00624818
##  [7,] -2.79506093
##  [8,] -1.42424492
##  [9,] -0.31368592
## [10,]  0.06411245
## [11,]  0.06292768
## 
## $CV
## [1] 0.8240508
## 
## $GCV
## [1] 0.7766156

Now, plot for all provinces, using the same smoothing parameters.

order = 4
der = 2
alpha = 0.5

for (P in unique(tmax_2015$province))
{
  compositionalSplineF(Ltmax2015CS[[P]]$t, 
                    Ltmax2015CS[[P]]$clrf,
                    Ltmax2015CS[[P]]$knots,
                    Ltmax2015CS[[P]]$w, order,
                    der, alpha, 
                    spline.plot = TRUE, 
                    basis.plot = FALSE,
                    namemain = P)
}