#---------------------------
# Packages

library(reshape2)
library(futile.matrix)
library(maptools)
library(mapproj)
library(gtools)
library(ncdf4)
#---------------------------


#---------------------------
# Convert g/km2/yr to kg/m2/s

g_to_kg = 1e-3
km2_to_m2 = 1e-6 # Convert /km2 to /m2
sec_per_year = 365*24*60*60
#---------------------------



#---------------------------
# Emissions

setwd('/Users/angot/Documents/Projects/GMA 2018/emissions_inventory/GMA2018_emissions_v2.0/__POWERGEN_orig_speciation/layers_g_km2/')

#-------
# Hg0

Hg0_0_50 = read.table("allSectors__Hg0_level0_50_kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
# We now need to convert the data frames into a matrix 
# dim: 1440x720 for a 0.25x0.25 resolution
Hg0_0_50 = as.matrix(Hg0_0_50)
Hg0_0_50 = t(Hg0_0_50)

Hg0_50_150 = read.table("allSectors__Hg0_level50_150_kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
Hg0_50_150 = as.matrix(Hg0_50_150)
Hg0_50_150 = t(Hg0_50_150)

Hg0_150 = read.table("allSectors__Hg0_level150___kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
Hg0_150 = as.matrix(Hg0_150)
Hg0_150 = t(Hg0_150)

# Convert g/km2/yr to kg/m2/s
Hg0_0_50 = Hg0_0_50 * g_to_kg * km2_to_m2 / sec_per_year
Hg0_50_150 = Hg0_50_150 * g_to_kg * km2_to_m2 / sec_per_year
Hg0_150 = Hg0_150 * g_to_kg * km2_to_m2 / sec_per_year

Hg0 = Hg0_0_50 + Hg0_50_150 + Hg0_150
#---------

#---------
# Hg2

Hg2_0_50 = read.table("allSectors__Hg2_level0_50_kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
# We now need to convert the data frames into a matrix 
# dim: 1440x720 for a 0.25x0.25 resolution
Hg2_0_50 = as.matrix(Hg2_0_50)
Hg2_0_50 = t(Hg2_0_50)

Hg2_50_150 = read.table("allSectors__Hg2_level50_150_kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
Hg2_50_150 = as.matrix(Hg2_50_150)
Hg2_50_150 = t(Hg2_50_150)

Hg2_150 = read.table("allSectors__Hg2_level150___kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
Hg2_150 = as.matrix(Hg2_150)
Hg2_150 = t(Hg2_150)

# Convert g/km2/yr to kg/m2/s
Hg2_0_50 = Hg2_0_50 * g_to_kg * km2_to_m2 / sec_per_year
Hg2_50_150 = Hg2_50_150 * g_to_kg * km2_to_m2 / sec_per_year
Hg2_150 = Hg2_150 * g_to_kg * km2_to_m2 / sec_per_year

Hg2 = Hg2_0_50 + Hg2_50_150 + Hg2_150
#---------

#---------
# HgP

HgP_0_50 = read.table("allSectors__HgP_level0_50_kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
# We now need to convert the data frames into a matrix 
# dim: 1440x720 for a 0.25x0.25 resolution
HgP_0_50 = as.matrix(HgP_0_50)
HgP_0_50 = t(HgP_0_50)

HgP_50_150 = read.table("allSectors__HgP_level50_150_kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
HgP_50_150 = as.matrix(HgP_50_150)
HgP_50_150 = t(HgP_50_150)

HgP_150 = read.table("allSectors__HgP_level150___kg_cell_g_km2.txt",header = FALSE, skip = 6, nrows = 720, dec = ".", na.strings = "-9999", colClasses = "numeric")
HgP_150 = as.matrix(HgP_150)
HgP_150 = t(HgP_150)

# Convert g/km2/yr to kg/m2/s
HgP_0_50 = HgP_0_50 * g_to_kg * km2_to_m2 / sec_per_year
HgP_50_150 = HgP_50_150 * g_to_kg * km2_to_m2 / sec_per_year
HgP_150 = HgP_150 * g_to_kg * km2_to_m2 / sec_per_year

HgP = HgP_0_50 + HgP_50_150 + HgP_150
#---------
#---------------------------


#---------------------------
# lon/lat format
# In the source file (.txt), the lower left corner corresponds to the point (-180,-90)
# meaning that columns: n = 1440, longitudes from -180 to + 180 (left to right)
# and rows: n = 720, latitudes from +90 to -90 (top to bottom)

# In GEOS-Chem we need the transpose of this matrix (done above)
# We now have columns: n = 720, latitudes from +90 to -90 (left to right)
# and rows: n = 1440, longitudes from -180 to 180 (top to bottom)

# We need to reorganize the order of the rows/columns

# We give names to the rows/columns according to the current order:
rows.i = seq(-179.875,179.875,0.25)
rows.i = as.character(round(rows.i,2))
cols.i = seq(89.875,-89.875,-0.25)
cols.i = as.character(round(cols.i,2))
dimnames(Hg0) = list(rows.i, cols.i)
dimnames(Hg2) = list(rows.i, cols.i)
dimnames(HgP) = list(rows.i, cols.i)

# We want the final order or the rows/columns to be:
lon_1 = seq(0.125,179.875,0.25)
lon_2 = seq(-179.875,-0.125,0.25)
rows.f = as.character(round(c(lon_1,lon_2),2))
cols.f = round(seq(-89.875,89.875,0.25),2)
cols.f = as.character(cols.f)

# We change the order of the rows/columns:
row.order = rows.f
Hg0 = Hg0[rows.f,]
Hg2 = Hg2[rows.f,]
HgP = HgP[rows.f,]
col.order = cols.f
Hg0 = Hg0[,col.order]
Hg2 = Hg2[,col.order]
HgP = HgP[,col.order]
#---------------------------


#---------------------------
# Creation .nc file

#---------
# Hg0

# Define dimensions of the new file
londim = ncdim_def("lon",units = "degrees_east",vals = seq(0.125,359.875,0.25))
latdim = ncdim_def("lat",units = "degrees_north",vals = seq(-89.875,89.875,0.25))
# Define variables of the new file
emi_hg_0 = ncvar_def("emi_hg_0",units = "kg m-2 s-1",dim = list(londim,latdim))
# Create file with associated dimensions
filename_new = "/Users/angot/Documents/Projects/GMA 2018/emissions_inventory/nc_files/GMA_emissions_POWERGEN_Hg0.0.25x0.25.2015.nc"
ncnew = nc_create(filename_new, emi_hg_0)
# Put variable in file
ncvar_put(nc = ncnew, varid = emi_hg_0, vals = Hg0)
nc_close(ncnew)
#---------

#---------
# Hg2 

# Define dimensions of the new file
londim = ncdim_def("lon",units = "degrees_east",vals = seq(0.125,359.875,0.25))
latdim = ncdim_def("lat",units = "degrees_north",vals = seq(-89.875,89.875,0.25))
# Define variables of the new file
emi_hg_2 = ncvar_def("emi_hg_2",units = "kg m-2 s-1",dim = list(londim,latdim))
# Create file with associated dimensions
filename_new = "/Users/angot/Documents/Projects/GMA 2018/emissions_inventory/nc_files/GMA_emissions_POWERGEN_Hg2.0.25x0.25.2015.nc"
ncnew = nc_create(filename_new, emi_hg_2)
# Put variable in file
ncvar_put(nc = ncnew, varid = emi_hg_2, vals = Hg2)
nc_close(ncnew)
#---------

#---------
# HgP

# Define dimensions of the new file
londim = ncdim_def("lon",units = "degrees_east",vals = seq(0.125,359.875,0.25))
latdim = ncdim_def("lat",units = "degrees_north",vals = seq(-89.875,89.875,0.25))
# Define variables of the new file
emi_hg_p = ncvar_def("emi_hg_p",units = "kg m-2 s-1",dim = list(londim,latdim))
# Create file with associated dimensions
filename_new = "/Users/angot/Documents/Projects/GMA 2018/emissions_inventory/nc_files/GMA_emissions_POWERGEN_HgP.0.25x0.25.2015.nc"
ncnew = nc_create(filename_new, emi_hg_p)
# Put variable in file
ncvar_put(nc = ncnew, varid = emi_hg_p, vals = HgP)
nc_close(ncnew)
#---------
#---------------------------