# This is a script to read and process LOCA-obs data
# Observations: 1/8 x 1/8 degree, 1950-2013
#===================================================================================================================
library(ncdf4)
library(caTools)
library(abind)
#===================================================================================================================

vars = c('pr','pet') 
dir = '/groups/anenberggrp/PA/dust/LOCA/obs/'

xlim=c(-125,-90); ylim=c(24,45) # Extract SW+Mex
years=c(1950:2013)
ts = seq(as.Date("1950/01/01"), as.Date("2013/12/01"), by = "month")
tmin = format(min(ts), format='%Y%m'); tmax = format(max(ts), '%Y%m')

for (j in 1:length(vars)){
  print(vars[j])
  files=(Sys.glob(paste(dir,vars[j],'/*.nc', sep='')))
  #df=array(NA, dim=c(229, 464, length(years)*12))

  raw.dat = nc_open(files[1])
    lon = ncvar_get(raw.dat, "centroidLons")
    lat = ncvar_get(raw.dat, "centroidLats")
    ind1=which(lon>=xlim[1] & lon<=xlim[2])
    ind2=which(lat>=ylim[1] & lat<=ylim[2])
    sub.lon = lon[ind1]; sub.lat=lat[ind2]

    df = ncvar_get(raw.dat, start=c(ind2[1],ind1[1],1), count=c(length(ind2),length(ind1),12))
    # lat x lon x months

  for (i in 2:length(files)){
    print(i)
    fname=files[i]
    raw.dat = nc_open(fname)
    x = ncvar_get(raw.dat, start=c(ind2[1],ind1[1],1), count=c(length(ind2),length(ind1),12))
    df = abind(df, x, along=3) 
  }
  print(dim(df))
  
  new.lat = rev(sub.lat)
  df = aperm(df, c(2,1,3)) # change to lon x lat x months
  df = df[,dim(df)[2]:1,] # reverse lat to be from S to N

  if (j==2){ # need to convert daily to monthly
    m1=c(1,3,5,7,8,10,12); m2=c(4,6,9,11); m3=2
    temp = array(df, dim=c(length(sub.lon), length(new.lat),12,dim(df)[3]/12))
    temp[,,m1,]=temp[,,m1,]*31; temp[,,m2,]=temp[,,m2,]*30; temp[,,m3,]=temp[,,m3,]*28
    df=array(temp, dim=dim(df))  
  }

  ss = list('data'=df, 'lon'=sub.lon, 'lat'=new.lat, 'time'=ts)
  print(dim(ss$data))

    #fillvalue = 1e+20 # missing value
    #x[x == fillvalue] = NA # set missing values to zero 

  fout = paste(vars[j],'_obs_',tmin,'-',tmax,'.Rdata', sep='') # filename
  save(ss, file=paste('./data/LOCA.obs/',fout,sep=''))
  print(fout)
}

#===================================================================================================================
#===================================================================================================================