from osgeo import gdal
gdal.UseExceptions()

import struct

objs = {}

def wakeup(filename):

    objs['ds'] = gdal.Open(filename)
    objs['gt'] = objs['ds'].GetGeoTransform()
    objs['band'] = objs['ds'].GetRasterBand(1)
    return objs


def read(x, y, objs):

    # stolen from http://stackoverflow.com/questions/24537450/python-struct-error-unpack-requires-a-string-argument-of-length-2
    gt = objs['gt']
    ds = objs['ds']
    band = objs['band']

    px = int((x - gt[0]) / gt[1])
    py = int((y - gt[3]) / gt[5])

    val = band.ReadRaster(px,py,1,1,buf_type=gdal.GDT_Float32)
    z = struct.unpack('f' , val)
    return z[0]

def filter(ins,outs):
    HAG = ins['HAG']
    X = ins['X']
    Y = ins['Y']
    Z = ins['Z']

    objs = wakeup('@CMAKE_SOURCE_DIR@/test/data/autzen/autzen-surface.tif.min.tif')
    print (objs)
    for i in range(len(X)):
        x = X[i]
        y = Y[i]
        z = Z[i]
        surface_z = read(x, y, objs)
#        print (x, y, z, surface_z)
        hag = z - surface_z
        if (surface_z  != -9999):
            HAG[i] = hag
    outs['Z'] = HAG
    return True

if __name__=='__main__':
    objs = wakeup('@CMAKE_SOURCE_DIR@/test/data/autzen/autzen-surface.tif.min.tif')
    print (read(636436,850412, objs))