## This script caculates the mortality due to PM2.5 in China
rm(list=ls())
setwd('~/Documents/Research/3Chn/NCC/Final/Li&Zhang_2018/Data_and_Scripts')

# China mask
library(ncdf4)
file = nc_open('Province_mask/mask_prv_05x0667.nc')
mask = ncvar_get(file, varid='mask')
nc_close(file)
ind.ct = which(mask>0)

##== Gridded pop
file = nc_open('Population/SEDAC_2010_05x0667.nc')
pop.grid  = ncvar_get(file, varid='pop')
nc_close(file)
# ct total pop
pop.ct = sum(pop.grid[ind.ct],na.rm=T)

##== PM concentration
file = nc_open('GEOS-Chem/GC_PM25.nc')
conc0 = ncvar_get(file, varid='conc_nst')
nc_close(file)

# Annual average conc
conc1=apply(conc0[,,,,2:13],1:4,mean,na.rm=TRUE)

# RH correction
conc = array(NA,c(dim(conc1)[1:3],7))
#sna
conc[,,,1:3] = conc1[,,,1:3]*1.33
#bc
conc[,,,4] = conc1[,,,4]
#oc
conc[,,,5] = (conc1[,,,5]*1.16 + conc1[,,,6])*1.8
#sst
conc[,,,6] = conc1[,,,7]*1.86 + conc1[,,,8]
#dust
conc[,,,7] = conc1[,,,9]

# Annual averge total/anth. PM
library(abind)
cct1 = apply(conc[,,,1:7],c(1,2,3),sum,na.rm=TRUE)
cct2 = apply(conc[,,,1:5],c(1,2,3),sum,na.rm=TRUE)
cct  = abind(cct1,cct2,rev.along=0)
# include 2010 cct twice to be consistent with cases
cct = abind(cct[,,1,],cct[,,1,],cct[,,2:dim(cct)[3],],rev.along=2)
cas = c('2010','2010c2030p','bau','cint3','cint4','cint5','bau_nst','cint3_nst','cint4_nst','cint5_nst')
dimnames(cct)[[3]] = cas

##== pop 30+ and 5- from UN in 2010 and 2030
pop0 = read.csv('Population/WPP2015_DB03_Population_By_Age_Quinquennial.csv')
year = c('2010','2030')

for (i in 1:2){
    pop1 = pop0[pop0$Variant=='Medium' & pop0$Sex=='Both' & pop0$Time==year[i],]
    ipp = match(c("Pop_0_4","Pop_30_34","Pop_75_79","Pop_100"),colnames(pop1))
    ict = match(c('China'),pop1$Location)
    pop2 = as.matrix(pop1[ict,ipp[2]:ipp[3]])
    pop3 = apply(pop1[ict,(ipp[3]+1):ipp[4]],1,sum,na.rm=TRUE)
    # 30+ by age group
    pop.1 = c(pop2,pop3)
    # 30+
    pop.3 = apply(pop1[ict,ipp[2]:ipp[4]],1,sum,na.rm=TRUE)
    # 5-
    pop.5 = pop1[ict,ipp[1]]

    if (i==1) pop.2010 = list(pop.1,pop.1,pop.3,pop.3,pop.5)
    else pop.2030 = list(pop.1,pop.1,pop.3,pop.3,pop.5)
}

##== Incidence
load(file = "Health/GBD/y0.Rdata")

##== GBD model
prm0 = read.csv('Health/GBD/GBD_model/IHME_CRCurve_parameters.csv')

dss = names(y0.ch)
nds = length(dss)
prm  = list()

for (ids in 1:nds){
    prm1 = prm0[prm0$cause==dss[ids],]
    if (ids > 2){
       prm[[ids]] = as.matrix(prm1[,5:8])
    }else{
       age0 = unique(prm1$age)
       age  = age0[3:length(age0)]
       prm3 = array(NA,c(1000,4,length(age)))
       for (iag in 1:length(age)){
           prm2 = prm1[prm1$age==age[iag],]
           prm3[,,iag] = as.matrix(prm2[,5:8])
       }
       prm[[ids]] = prm3
    }
}

##== Mort using 2010/2030 pop in thousands of deaths
ngd = length(ind.ct)
ncs = length(cas)
npm = 2 #[total,anth.]

mt  = array(NA,c(ngd,ncs,npm,nds,1000))

dimnames(mt)[[2]] = cas
dimnames(mt)[[3]] = c('total','anth')
dimnames(mt)[[4]] = dss

for (ipm in 1:npm){
for (ics in 1:ncs){
    print(ics)
    cct0 = cct[,,ics,ipm]
    if (ics==1) pop0 = pop.2010 else pop0 = pop.2030

    for (ids in 1:nds){
        prm0 = prm[[ids]]
        pop.un = pop0[[ids]]
        y0   = y0.ch[[ids]]

        for (igd in 1:ngd){
            pm = cct0[ind.ct[igd]]
            if (ids > 2){             
               RR = 1+prm0[,1]*(1-exp(-prm0[,2]*(pm-prm0[,4])^prm0[,3]))
               # RR.itv = quantile(RR, c(.025, .50, .975), na.rm=TRUE)  
               mt[igd,ics,ipm,ids,] = y0*((RR-1)/RR)*pop.grid[ind.ct[igd]]/pop.ct*pop.un
            } else{
               mt0 = array(NA,c(11,1000)) #(age group)
               RR = 1+prm0[,1,]*(1-exp(-prm0[,2,]*(pm-prm0[,4,])^prm0[,3,]))
               for (iag in 1:11){
                   # RR.itv = quantile(RR, c(.025, .50, .975), na.rm=TRUE)  
                   mt0[iag,] = y0[iag]*((RR[,iag]-1)/RR[,iag])*pop.grid[ind.ct[igd]]/pop.ct*pop.un[iag]
               }
               mt[igd,ics,ipm,ids,] = apply(mt0,2,sum,na.rm=TRUE) 
            }       
        }   
    }
}
}

# Median index
mt.ch = apply(mt,c(2:5),sum,na.rm=TRUE)
whichmedian <- function(x) which.min(abs(x - median(x)))
ind.med = apply(mt.ch,c(1:3),whichmedian)

# Median value by grid cell
mt.med = array(NA,c(ngd,ncs,npm,nds))
for (ics in 1:ncs){
	for (ipm in 1:npm){
		for (ids in 1:nds){
			mt.med[,ics,ipm,ids]=mt[,ics,ipm,ids,ind.med[ics,ipm,ids]]
		}
	}
}

mt.ch.tt = apply(mt.med,2:4,sum,na.rm=TRUE)
dimnames(mt.ch.tt)[[1]] = cas
dimnames(mt.ch.tt)[[2]] = c('total','anth')
dimnames(mt.ch.tt)[[3]] = dss

print(round(mt.ch.tt['2010','total',],0))
   # IHD STROKE   COPD     LC   ALRI 
   # 232    617    238    112     11
      
  
##== Calculate prv mortality and save in csv files
mt.med.tt = apply(mt.med,1:3,sum,na.rm=TRUE) #ngd,ncs,npm

library(xlsx)
prv.nms = read.xlsx("Province_mask/mask_prv_nms.xlsx",1)
npv = dim(prv.nms)[1]
mt.prv = array(NA,c(npv,ncs))

# Re-order rows based on r2 
data = read.xlsx("../Scripts_for_Figures/plot.xlsx",1)
name1 = data[,'r2']

for (i in 1:npm){
    for (j in 1:npv){
        ind.prv = match(which(mask==prv.nms$number[j]),ind.ct)
        # total pop
        # mt.prv[j,1] = round(sum(pop.grid[ind.ct[ind.prv]],na.rm=T),0)
        # mortality
        mt.prv[j,1:ncs] = round(apply(mt.med.tt[ind.prv,,i]*1000,2,sum,na.rm=T),0)
    }
    rownames(mt.prv) = prv.nms$short_name
    colnames(mt.prv) = cas
    
    # Remove Tibet and add Chinese total
    ind = which(prv.nms$long_name=='Tibet')
    mt.prv = mt.prv[-ind,]
    mt.chn = apply(mt.prv,2,sum,na.rm=T)
    mt.prv = rbind(mt.prv,mt.chn)
    rownames(mt.prv)[31] = 'China'

    # Reorder rows
    name2 = rownames(mt.prv)
    ind = match(name1,name2)
    mt.prv = mt.prv[ind,]
    
    if (i==1) write.xlsx(mt.prv, file = "Health/GBD/mortality.xlsx", sheetName="total") else write.xlsx(mt.prv, file = "Health/GBD/mortality.xlsx", sheetName="anth",append=TRUE)
}