#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Wed Aug 23 15:15:49 2017

@author: mschull
"""

import os
import subprocess
import h5py
import numpy as np
import csv
import pandas as pd
import glob
import datetime
import gzip
import shutil
from osgeo import gdal,osr
import argparse
import urllib2, base64
from bs4 import BeautifulSoup
import requests
import urllib
from pydap.client import open_url
from pydap.cas import urs
from joblib import Parallel, delayed
import time as timer
import math
import ephem


ncdcURL = 'https://nomads.ncdc.noaa.gov/modeldata/cfsv2_analysis_pgbh/'
ncdcfluxURL = 'https://nomads.ncdc.noaa.gov/modeldata/cfsv2_analysis_flxf/'

realtimeURL = 'http://nomads.ncep.noaa.gov/pub/data/nccf/com/cfs/prod/'

def folders(base):
    data_path = os.path.abspath(os.path.join(base,os.pardir,'VIIRS_DATA'))
    if not os.path.exists(data_path):
        os.makedirs(data_path) 
    static_path = os.path.join(base,"STATIC")
    if not os.path.exists(static_path):
        os.makedirs(static_path) 
    CFSR_path = os.path.join(static_path,"CFSR")   
    if not os.path.exists(CFSR_path):
        os.makedirs(CFSR_path)
    out = {'data_path':data_path,'static_path':static_path,'CFSR_path':CFSR_path}
    return out

base = os.getcwd()
Folders = folders(base)
data_path = Folders['data_path']
static_path = Folders['static_path']
CFSR_path = Folders['CFSR_path']

def warp(args):
    """with a def you can easily change your subprocess call"""
    # command construction with binary and options
    options = ['gdalwarp']
    options.extend(args)
    # call gdalwarp 
    subprocess.check_call(options)

def writeArray2Tiff(data,res,UL,inProjection,outfile,outFormat):

    xres = res[0]
    yres = res[1]

    ysize = data.shape[0]
    xsize = data.shape[1]

    ulx = UL[0] #- (xres / 2.)
    uly = UL[1]# - (yres / 2.)
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(outfile, xsize, ysize, 1, outFormat)
    #ds = driver.Create(outfile, xsize, ysize, 1, gdal.GDT_Int16)
    
    srs = osr.SpatialReference()
    
    if isinstance(inProjection, basestring):        
        srs.ImportFromProj4(inProjection)
    else:
        srs.ImportFromEPSG(inProjection)
        
    ds.SetProjection(srs.ExportToWkt())
    
    gt = [ulx, xres, 0, uly, 0, -yres ]
    ds.SetGeoTransform(gt)
    
    ds.GetRasterBand(1).WriteArray(data)
    #ds = None
    ds.FlushCache() 

def is_odd(num):
   return num % 2 != 0

def getGrabTime(time):    
    return int(((time/600)+1)*600)


def getGrabTimeInv(grab_time,doy):
    if is_odd(grab_time):
        hr = grab_time-3
        forecastHR = 3
    else:
        hr = grab_time
        forecastHR = 0
    if hr == 24:
        hr = 0
        doy+=1
    return hr, forecastHR,doy 

def listFD(url, ext=''):
    page = requests.get(url).text
    soup = BeautifulSoup(page, 'html.parser')
    return [url + '/' + node.get('href') for node in soup.find_all('a') if node.get('href').endswith(ext)]

def tile2latlon(tile):
    row = tile/24
    col = tile-(row*24)
    # find lower left corner
    lat= (75.-row*15.)-15.
    lon=(col*15.-180.)-15. 
    return [lat,lon]

def get_VIIRS_bounds(fn):
    f = h5py.File(fn, 'r')
    east = []
    west = []
    north = []
    south = []
    for i in range(4):
        east = np.append(east,f['Data_Products']['VIIRS-I5-SDR']['VIIRS-I5-SDR_Gran_%d' % i ].attrs['East_Bounding_Coordinate'][0][0])
        west = np.append(west,f['Data_Products']['VIIRS-I5-SDR']['VIIRS-I5-SDR_Gran_%d' % i ].attrs['West_Bounding_Coordinate'][0][0])
        north = np.append(north,f['Data_Products']['VIIRS-I5-SDR']['VIIRS-I5-SDR_Gran_%d' % i ].attrs['North_Bounding_Coordinate'][0][0])
        south = np.append(south,f['Data_Products']['VIIRS-I5-SDR']['VIIRS-I5-SDR_Gran_%d' % i ].attrs['South_Bounding_Coordinate'][0][0])
    east = east.min()
    west = west.max()
    north = north.max()
    south = south.min()
    date = f['Data_Products']['VIIRS-I5-SDR']['VIIRS-I5-SDR_Gran_0'].attrs['Beginning_Date'][0][0]
    N_Day_Night_Flag = f['Data_Products']['VIIRS-I5-SDR']['VIIRS-I5-SDR_Gran_0'].attrs['N_Day_Night_Flag'][0][0]
    bounds = {'filename':fn,'east':east,'west':west,'north':north,'south':south,
              'N_Day_Night_Flag':N_Day_Night_Flag,'date':date}
    database = os.path.join(data_path,'I5_database.csv')
    if not os.path.exists(database):
        with open(database, 'wb') as f:  # Just use 'w' mode in 3.x
            w = csv.DictWriter(f, bounds.keys())
            w.writeheader()
            w.writerow(bounds)
    else:
        print("writing")
        with open(database, 'ab') as f:  # Just use 'w' mode in 3.x
            w = csv.DictWriter(f, bounds.keys())
            w.writerow(bounds)
            
    return        

def downloadSubscriptionSDR(year=None,doy=None,url=None): 
    if year==None:
        dd = datetime.date.today()+datetime.timedelta(days=-1)
        year = dd.year
    if doy==None:
        dd = datetime.date.today()+datetime.timedelta(days=-1)
        month = dd.month
        day = dd.day
    else:
        dd=datetime.datetime(year,1,1)+datetime.timedelta(days=doy)
        month = dd.month
        day = dd.day
    filePath = os.path.join(data_path,"%d" % year,"%02d" % month)
    #download I5 data
    ext = 'h5'
    if url==None: # Use subscription
        url = 'https://download.class.ngdc.noaa.gov/download/sub/hain/85113/'

    for fn in listFD(url, ext):
        print fn
        fileName = str(fn.split('/')[-1])  
        if (fileName.split("_")[2]=='d%d%02d%02d' % (year,month,day)):
#            filePath = os.path.join(outPath,'%s' % fileName.split('_')[0])
            if not os.path.exists(filePath):
                os.makedirs(filePath)
        
            outName=os.path.join(filePath,fileName)
            
            if not os.path.isfile(outName):
                print "downloading:  %s" % fileName
                #wget.download(url+fileName,out=outName)
                urllib.urlretrieve(url+fileName, outName)

    # download cloud data
    if url==None:
        url = 'https://download.class.ngdc.noaa.gov/download/sub/hain/85123/'
        
        for fn in listFD(url, ext):
            print fn
            fileName = str(fn.split('/')[-1])  
            if (fileName.split("_")[2]=='d%d%02d%02d' % (year,month,day)):
    #            filePath = os.path.join(outPath,'%s' % fileName.split('_')[0])
                if not os.path.exists(filePath):
                    os.makedirs(filePath)
            
                outName=os.path.join(filePath,fileName)
                
                if not os.path.isfile(outName):
                    print "downloading:  %s" % fileName
                    #wget.download(url+fileName,out=outName)
                    urllib.urlretrieve(url+fileName, outName)
                
                
def getInsolation(earthLoginUser,earthLoginPass,tile,year=None,doy=None):
    if year==None:
        dd = datetime.date.today()+datetime.timedelta(days=-1)
        year = dd.year
    if doy==None:
        doy = (datetime.date.today()-datetime.date(year,1,1)).days
        dd = datetime.date.today()+datetime.timedelta(days=-1)
        month = dd.month
        day = dd.day
    llLat,llLon = tile2latlon(tile)
    ulx = llLon
    uly = llLat+15.
    lrx = llLon+15.
    lry = llLat
    MERRA2_ulLat = 90.0
    MERRA2_ulLon = -180.0
    MERRA2LatRes = 0.5
    MERRA2LonRes = 0.625
    nrow = 3750
    ncol = 3750
    date = '%d%03d' % (year,doy)
    inProj4 = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
    #========get MERRA2 Insolation data at overpass time====================== 
    inRes = [MERRA2LonRes,MERRA2LatRes]
    inUL = [MERRA2_ulLon,MERRA2_ulLat]

    if year <1992:
        fileType = 100
    elif year >1991 and year < 2001:
        fileType=200
    elif year > 2000 and year<2011:
        fileType = 300
    else: