'''
Extract given list of sensors within given time period
Replace 4 character station ID with a integer no
Conver Lat & log to y & x based on given reference point
Chage time stamp to a relative one
'''

import math

def haversine(origin, destination):
    [lat1, lon1] = origin
    [lat2, lon2] = destination
    radius = 6371.009 # km

    dLat = math.radians(lat2 - lat1)
    dLon = math.radians(lon2 - lon1)
    a = math.sin(dLat/2) * math.sin(dLat/2) + math.cos(math.radians(lat1)) \
        * math.cos(math.radians(lat2)) * math.sin(dLon/2) * math.sin(dLon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = radius * c

    return d    


fdIn = open('sensorsAll.txt', 'r')
fdOut = open('subsetSensors.txt', 'w')

idx2Extract = [1, 0, 24, 22, 20, 6, 8, 10, 12, 14]
#idx2Extract = [1, 0, 24, 22, 20, 6, 8, 10, 12, 14, 26, 28, 30, 32, 34, 36, 38, 40, 42, 64, 65]
attNames = 'SID\tTime\tX\tY\tZ\tRH\tP\tT\tDD\tFF\n'
#attNames = 'SID\tTime\tX\tY\tZ\tRH\tP\tT\tDD\tFF\tNS1\tNS2\tNS3\tNS4\tNS5\tNS6\tNS7\tNS8\tPRESWEA\tSKYCOV\n'
sTime = 1340946000 
eTime = sTime + 86400 * 3
offset = 300
ll = [32.02, -102.98] #Lower left corner of sensor filed. Use Latitudes & longitudes
sID = {}
count = 0

fdOut.write(attNames)

for line in fdIn:
    #remove extra space, replace missing data, & split
    line = line[:-1] #Remove new line
    l = line.replace('    ', ' ').replace('-99999.000000', '-1').split(',')[:-1]
    if l[6] == '-1' and l[8] == '-1' and l[10] == '-1' and l[12] == '-1' and l[14] == '-1': continue #Skip missing ones
    t = int(float(l[0])) #time
    if t < sTime: continue
    if t > eTime: break
    t = t - sTime + offset #Cobnvert to relative time with given offset
    l[0] = str(t)

    if l[1] not in sID: #New numeric sensor ID
        sID[l[1]] = str(count)
        count += 1
    l[1] = sID[l[1]] 

    lat = haversine(ll, [float(l[22]), ll[1]])* 1000
    lon = haversine(ll, [ll[0], float(l[24])])* 1000
    if lat > 1000000.0 or lon > 1000000.0 : print line
    l[24] = str(round(lon, 4)) #Relative X location in meters
    l[22] = str(round(lat, 4)) #Y location
    tmpStr = ''
    for i in idx2Extract:
        if i != idx2Extract[-1]:
            tmpStr += l[i] + '\t'
        else:
            tmpStr += l[i] + '\n'
    fdOut.write(tmpStr)
fdIn.close()
fdOut.close()
