inStations = r'I:\Documents\Klusjes\GIDS Interpolation Air Temperature\GSOD//NEARNLSTATIONS_MERGED_8_monday.csv'
folderGSOD = r'I:\Documents\Klusjes\GIDS Interpolation Air Temperature\GSOD2009\selection_nl_2'
RasterDEM = r'I:\Documents\Klusjes\GIDS Interpolation Air Temperature\GSOD//AALBERS_25KM_DEM.tif'
outFolder = r'I:\Documents\Klusjes\GIDS Interpolation Air Temperature\GSOD2009\outTifs3//'
prefix = 'WDSP'
date_all = [20120723, 20120724, 20120725, 20120726, 20120727, 20120728, 20120729, 20120730]

from scipy import stats
from scipy.spatial.distance import cdist
from scipy.spatial import cKDTree as KDTree
import statsmodels.api as sm
import numpy as np
import pandas as pd
import gdal
import os
from __future__ import division

def GSODfiles(inGSODFolder):
    st_wmo = [os.path.join(root, name)
               for root, dirs, files in os.walk(inGSODFolder)
                 for name in files
                 if name.endswith('.op')]
    return st_wmo

def inRaster(fileDEM):
    raster = gdal.Open(fileDEM, gdal.GA_ReadOnly)
    band = raster.GetRasterBand(1)
    dem = band.ReadAsArray()
    extent = raster.GetGeoTransform()
    return raster, dem, extent

def saveRaster(path, array, datatype=6, formatraster="GTiff"):
    # Set Driver
    format_ = formatraster #save as format
    driver = gdal.GetDriverByName( format_ )
    driver.Register()
    
    # Set Metadata for Raster output
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    bands = raster.RasterCount
    datatype = 6#band.DataType
    
    # Set Projection for Raster
    outDataset = driver.Create(outFilename, cols, rows, bands, datatype)
    geoTransform = raster.GetGeoTransform()
    outDataset.SetGeoTransform(geoTransform)
    proj = raster.GetProjection()
    outDataset.SetProjection(proj)
    
    # Write output to band 1 of new Raster
    outBand = outDataset.GetRasterBand(1)
    outBand.WriteArray(array) #save input array
    #outBand.WriteArray(dem)
    
    # Close and finalise newly created Raster
    #F_M01 = None
    outBand = None
    proj = None
    geoTransform = None
    outDataset = None
    driver = None
    datatype = None
    bands = None
    rows = None
    cols = None
    driver = None
    array = None

def IDWKDtree(x,y):    
    x_y = [(x,y)] # for kd_tree that starts counting at bottom left with 0,0
    #print x_y
    x_ = (x*extent[1]+extent[0]+extent[1]/2) # longitude aalbers projection (meters)
    y_ = (y*extent[5]+extent[3]+extent[5]/2) # latitude aalbers projection (meters)
    long_lat = np.array([[x_,y_]])
    
    #print long_lat
    
    dem_1 = dem[y,x] # elevation x_y coordinate
    dist_tree, ix_tree = tree.query(x_y, k=8, eps=0, p=1) # returns distance and index
    dist_tree = dist_tree.flatten()
    ix_tree  = ix_tree.flatten()
    df_selection = df.ix[ix_tree.ravel()] # seleclt 8 neighbours
    c = df_selection.ix[:,14]
    
    if dist_tree[0] < 1e-10:
        wz = c.ix[ix_tree[0]]
    else:
        w = 1 / dist_tree**2
    w /= np.sum(w)
    wz = np.dot(w, c.ix[ix_tree])
    
    #df_selection = df.ix[ix_tree.ravel()]

    return wz

st_wmo = GSODfiles(folderGSOD)
raster, dem, extent = inRaster(RasterDEM)
for date in date_all:

    df = pd.read_csv(inStations)
    #df_select = None
    df_select = pd.DataFrame()
    for st in st_wmo:
        with open(st) as f:
            """
            0,6     = STNN---  : USAF number
            14,22   = YEARMODA : YearMonthDay
            24,30   = TEMP     : Mean temperature in Fahreinheit
            35,41   = DEWP     : Mean dew point in Fahreinheit
            57,63   = STP      : Mean station pressure in millibars
            78,83   = WDSP     : Mean wind speed in knots
            118,123 = PRCP     : Total precipitation in inches

            """
            colspecs_data = [(0,6), (14,22), (24,30), (35,41), (57,63), (78,83), (118,123)]
            df_date = pd.read_fwf(f, colspecs=colspecs_data)
            for moda in df_date.YEARMODA:
                #print moda
                if moda == date:
                    df_select = df_select.append(df_date[df_date.YEARMODA==date])
                    break

    # Set Missing values
    df_select.replace({'STP':{9999.9:np.nan}}, inplace=True)
    df_select.replace({'TEMP':{9999.9:np.nan}}, inplace=True)
    df_select.replace({'DEWP':{9999.9:np.nan}}, inplace=True)
    df_select.replace({'WDSP':{999.9:np.nan}}, inplace=True)
    df_select.replace({'PRCP':{99.99:0}}, inplace=True)

    #print df_select.head()
    # Merge values with stations
    df = pd.merge(df, df_select, how='inner', left_on='USAF', right_on='STN---')
    #print df.head()

    ##STP
    #df = df[pd.notnull(df['STP'])]

    ##TEMP
    #df = df[pd.notnull(df['TEMP'])]
    #print df.head()

    ##DEWP
    #df = df[pd.notnull(df['DEWP'])]
    #df.reset_index(drop=True, inplace=True)
    #print df.head()

    ##WDSP
    df = df[pd.notnull(df['WDSP'])]
    df.reset_index(drop=True, inplace=True)

    ##PRCP
    #df = df[pd.notnull(df['PRCP'])]

    Longscaled = df.ix[:,8]
    #print 'longscaled\n', Longscaled
    Latscaled = df.ix[:,7]
    #print '\nlatscaled\n', Latscaled
    tree = KDTree(zip(Longscaled,Latscaled), leafsize=11)

    tp = np.zeros([dem.shape[1],dem.shape[0]])

    #x = 0
    #y = 0
    #GIDS(x,y)
    for x in range(0,20,1):
        for y in range(0,19,1):
            tp[x][y] = IDWKDtree(x,y)
    tp = tp.T

    #save output as raster
    outFilename = outFolder+prefix+str(date)+'.tif'
    saveRaster(outFilename, tp)
    outFilename = None
    tp = None

df.shape