#!/bin/bash
#***************************************************
#     CALC YEARLY ANNUAL AVERAGES FROM DAILY DATA
#***************************************************

#** Check that there was an output directory supplied
if [[ $# -eq 0 ]] ; 
then
  echo 'Usage: ./climatology.sh output_directory'
  exit 1
fi


OUTDIR=$1


if [ -d $OUTDIR ] ;
then 
  echo "Creating folder for temporary files"
  TMPDIR=${OUTDIR}/tmp
  if [ ! -d $TMPDIR ] ;
  then
    mkdir ${TMPDIR}
  fi
fi
#Clear log file 
> log.txt

URL="http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/cpc_us_precip/RT/precip.V1.0"

function createPrecip() {
for year in {2007..2014}
do
  echo $year
  #Clip the input to (-125.4,44.0),(-115.0,50.0)
  tmpSubset=${TMPDIR}/subset.${year}.nc
  echo "Subsetting precipitation data to study area to ${tmpSubset}"
  ncks -O -d lon,234.0,244.0 -d lat,44.0,50.0 ${URL}.${year}.nc ${tmpSubset}
  echo "" 
  echo "Finished subsetting the data, setting min/max and null values"
  echo ""
  echo `ncdump -h -v precip ${tmpSubset}`
  
  echo ""
  echo "Finished subsetting ${tmpSubset}, removing Null values..."
  tmpNaRm=${TMPDIR}/narm.${year}.nc
  ncap2 -O -s 'where(precip < 0.0 || precip > 10000 ) precip=0;' ${tmpSubset} ${tmpNaRm}


  echo ""
  echo "Finished subsetting ${tmpNaRm}, summarizing the daily precipitation data..."
  echo ""
  #Summarize the daily precipitation data for the last 7 years
  tmpAnnual=${TMPDIR}/annual.${year}.nc
  ncra -O --op_typ=ttl -d time,,,1 -v precip ${tmpNaRm} ${tmpAnnual}

  echo ""
  echo "Regridding ${tmpAnnual} to 1km resolution"
  echo ""
  #Regrid the annual data to 1km
  tmpRegrid=${TMPDIR}/regrid.${year}.nc
  export FILE=${tmpAnnual}
  export  OUTFILE="${tmpRegrid}"
  export OUTDIR="${TMPDIR}"
  ncl regrid.ncl
  echo ""
  echo "Finished regridding ${tmpRegrid}"
  echo ""
done
avgAnnual=${TMPDIR}/avg_annual.nc
echo "Averaging annuals precipitation to ${avgAnnual}"
nces -O -n 8,2,1 ${TMPDIR}/regrid.2007.nc $avgAnnual


}

createPrecip > log.txt

################################################
##     This R script will take open a 
##     netcdf-4 file and output it to a geotiff
##     
#################################################

library(raster)
library(ncdf4)
library(rgdal)

setwd('/home/ncasler/bootcamp/nco/tmp')


r <- raster('avg_annual.nc')

oldExt <- extent(r)


#Subtract 360 from the x coordinates
# will make coordinates lat/lon friendly
xMin <- oldExt[1] - 360
xMax <- oldExt[2] - 360
yMin <- oldExt[3]
yMax <- oldExt[4]

geoExt <- extent(xMin, xMax, yMin, yMax)
r <- setExtent(r, geoExt, keepres=TRUE)
newProj <- '+proj=lcc +lat_1=47.33333333333334 +lat_2=45.83333333333334 +lat_0=45.33333333333334 +lon_0=-120.5 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +to_meter=0.3048006096012192 +no_defs' 
xMinPrj <- 548554.658
xMaxPrj <- 2558498.351
yMinPrj <- 77271.774
yMaxPrj <- 1366159.434

xRes <- 1000
yRes <- 1000

newRaster <- raster(xmn=xMinPrj, xmx=xMaxPrj, ymn=yMinPrj, ymx=yMaxPrj, crs=newProj, res=1000)

newR <- projectRaster(r, newRaster)

print("Finished projecting raster")

summary(newR)
fileout <- "testRout.tif"
print(paste0("Writing new Raster to ", fileout))

writeRaster(newR, filename=fileout, format="GTiff", options=c("COMPRESS=ZIP"), overwrite=TRUE)

print("Finished outputting raster. Exiting...")

;***********************************************
; This script will read in netcdf precipitation files
;***********************************************
; Load the necessary scripts
;***********************************************
load "$NCARG_NCARG/nclscripts/csm/gsn_code.ncl"
load "$NCARG_NCARG/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_NCARG/nclscripts/csm/contributed.ncl"
load "$NCARG_NCARG/nclscripts/csm/shea_util.ncl"
;***********************************************
begin
;***********************************************
; Read in the netCDF file

  

  fn = getenv("FILE")
  print("Input file is " + fn)
  fileout = getenv("OUTFILE")
  print("Output file is " + fileout)


  fi = addfile(fn, "r")

  setfileoption("nc", "Format", "NetCDF4Classic")

;***********************************************
; Read in original dimensions
;***********************************************
  time = fi->time
;   Remove the time dimension from the precip variable(dimension order is time,lon,lat)
  precip = fi->precip(0,:,:)
  print("Precip max = "+ max(precip))
  print("Precip min = "+ min(precip))
 
 ; Read in the latitude and longitude dimensions
  lat = fi->lat
  lon = fi->lon
  printVarSummary(precip)
  
;**********************************************
; Create new dimension definitions
;**********************************************
;Set new extent
;*********************************************

  newXMin = 360 - 125.0463254981687271
  newYMin = 45.4641026583566514
  newXMax = 360 -  116.6757915137716282
  newYMax = 49.0780358101814755
;**********************************************
; Create new dimensions
;**********************************************
  newlat = fspan(newYMin,newYMax,1312)
  newlon = fspan(newXMin,newXMax,2088)
  newlat@units = "degrees_north"
  newlat@long_name = "Latitude"
  newlat@actual_range = (/newXMin, newXMax/)
  newlat@standard_name = "latitude"
  newlat@axis = "Y"

  newlon@units = "degrees_east"
  newlon@long_name = "Longitude"
  newlon@actual_range = (/newYMin, newYMax/)
  newlon@standard_name = "longitude"
  newlon@axis = "X"


  newprecip = linint2(fi&lon,fi&lat,precip,False,newlon,newlat,0)
;****************************************************
; Change all missing values to 0
;****************************************************
  printVarSummary(newprecip)
  print("NewPrecip Max: "+max(newprecip))
  print("NewPrecip Min: "+min(newprecip))
  newprecip!0 = "lat"
  newprecip!1 = "lon"
  newprecip&lat = newlat
  newprecip&lon = newlon
  newprecip@valid_range = (/min(newprecip), max(newprecip)/)
  newprecip@long_name = "Average Annual Precipitation 2007-2014"
  newprecip@units = "mm"
  newprecip@var_desc = "Precipitation"
  newprecip@level_desc = "Surface"
  newprecip@statistic = "Average Annual Accumulation"

;*************************************************
; Settings for output file
;************************************************

  nlat = dimsizes(newlat)           ; get dimension sizes
  nlon = dimsizes(newlon)
  
  print("File out is: "+fileout)
  system("/bin/rm -f " + fileout) ; remove if exists
  fout = addfile(fileout, "c")   ; open output file
;************************************************
; explicitly declare file definition mode. Improve efficiency
;************************************************
  setfileoption(fout,"DefineMode",True)

;************************************************
; create global attributes of the output file
;************************************************
  fAtt                = True        ;assign file attributes
  fAtt@title          = "Interpolated Annual Precipitation"
  fAtt@source_file    = fn
  fAtt@Conventions    = "COARDS"
  fAtt@description    = "Interpolated annual Precipitation"
  fAtt@creation_date  = systemfunc  ("date")
  fileattdef( fout, fAtt )          ; copy file attributes

;*************************************************
; predefine the coordinate variables and their dimensionality
; Note: to get an UNLIMITED record dimension, we set the dimensionality
; to -1 (or the actual size) and set the dimension name to True.
;*************************************************
  dimNames  = (/"lat", "lon"/)
  dimSizes  = (/ nlat, nlon/)
  dimUnlim  = (/ False, False/)
  filedimdef(fout,dimNames,dimSizes,dimUnlim)

;*************************************************
; Predefine the dimensionality of the variables to be written out
;*************************************************
; Here we are using NCL functions to facilitate defining
; each variable's dimension name(s) and type.
; The following could be replaced with explicit, user-defined dimension
; names different from those associated with the variable in memory
; Say, PS(time,lat,lon) in the NCL script. They could be redefined for the file via: 
; filevardef(fout, "PS", typeof(PS), (/"TIME", "latitude", "longitude"/))
;*************************************************

  filevardef(fout, "lat", typeof(lat), getvardims(lat))
  filevardef(fout, "lon", typeof(lon), getvardims(lon))
  filevardef(fout, "precip", typeof(newprecip), getvardims(newprecip))

;*************************************************
; Copy attributes associated with each variable to the file
; All attributes associated with each variable will be copied.
;*************************************************
  filevarattdef(fout,"precip",newprecip)      ; copy newprecip attributes
  filevarattdef(fout,"lat", newlat)           ; copy lat attributes
  filevarattdef(fout,"lon", newlon)           ; copy lon attributes

;*************************************************
; explicitly exit file definition mode. **NOT REQUIRED**
;*************************************************
  setfileoption(fout,"DefineMode", False)

;*************************************************
; output only the data values since the dimensionality and such have
; been predefined. The "(/","/)" syntax tells NCL to only output the
; data values to the predefined locations on the file
;**************************************************
  fout->lat     = (/newlat/)
  fout->lon     = (/newlon/)
  fout->precip  = (/newprecip/)
  print(fout)
end