targets2shp

import xlrd
import datetime
import sys
import shapefile as sf
import numpy as np
from math import pi, cos, sin, radians, degrees, atan2
import math

def nautical(kilometers=0, meters=0, miles=0, feet=0):
    nm = 0.
    if feet:
        miles += feet / ft(1.)
    if miles:
        kilometers += km(miles=miles)
    if meters:
        kilometers += meters / 1000.
    nm += kilometers / 1.852
    return nm

''' Distance and azimuth '''
def dist(llat1,llong1,llat2,llong2):
    rad = 6372795

    lat1 = llat1*math.pi/180.
    lat2 = llat2*math.pi/180.
    long1 = llong1*math.pi/180.
    long2 = llong2*math.pi/180.

    cl1 = math.cos(lat1)
    cl2 = math.cos(lat2)
    sl1 = math.sin(lat1)
    sl2 = math.sin(lat2)
    delta = long2 - long1
    cdelta = math.cos(delta)
    sdelta = math.sin(delta)

    y = math.sqrt(math.pow(cl2*sdelta,2)+math.pow(cl1*sl2-sl1*cl2*cdelta,2))
    x = sl1*sl2+cl1*cl2*cdelta
    ad = math.atan2(y,x)
    dist = ad*rad

    x = (cl1*sl2) - (sl1*cl2*cdelta)
    y = sdelta*cl2
    z = math.degrees(math.atan(-y/x))

    if (x < 0):
        z = z+180.

    z2 = (z+180.) % 360. - 180.
    z2 = - math.radians(z2)
    anglerad2 = z2 - ((2*math.pi)*math.floor((z2/(2*math.pi))) )
    angledeg = (anglerad2*180.)/math.pi

    #print 'Distance >> %.0f %s' % (dist/1000., ' [km]')
    #print 'Initial bearing >> ', angledeg, '[degrees]'
    return dist, angledeg

ff = "updated_targets.xls" #sys.argv[1]
width = []
length = []
height = []
speed = []
direction = []
type = []

###########################
#### Shapefile creation
###########################
shp_filename = 'shp\\%s' % ff.split('.')[0] #'updated_targets.shp'
# create the shapefile
w = sf.Writer(sf.POINT)
w.field('name','C','90')
w.field('type','C','90')
w.field('date','C','90')
w.field('direction','C','90')
w.field('speed','C','90')
w.field('height','C','90')
w.field('width','C','90')
w.field('length','C','90')
w.field('mass','C','90')
w.field('sources','C','90')
w.field('note','C','200')
############################

# Model file
f_model_flag =0
try:
    f_model = open('T:\\KARA\\ICEBERG_FORECAST\\target.txt','w')
    date_model = open('T:\\KARA\\ICEBERG_FORECAST\\data.txt','w')
    f_model_flag == 1
except:
    print "Cannot mount disk T !"

iid = 0

rb = xlrd.open_workbook(ff,formatting_info=True)
sheet = rb.sheet_by_index(0)

names = []

for rownum in range(sheet.nrows):
    row = sheet.row_values(rownum)
    try:
        names.append(row[0])
    except:
        pass
    #for c_el in row:
    #    print c_el

unique_names = set(names)

ilist = []
# lat/lons lists
llat = []
llon = []
names = []

ff_last_copd_positions = open('last_copd_positions.txt','w')

for iname in unique_names:
    print iname
    for rownum in range(sheet.nrows):
        row = sheet.row_values(rownum)
        if row[0]==iname:
            ilist.append({'name':row[0], 'type':row[2], 'date':row[4], \
                      'lat':row[5], 'lon':row[6], 'direction':row[7], \
                      'speed':row[8], 'height':row[9], 'width':row[10], \
                      'length':row[11], 'mass':row[12], 'sources':row[13], \
                      'note':row[14]})
    try:
        sorted_list = sorted(ilist, key=lambda x: datetime.datetime.strptime(x['date'],\
                                                                     '%m/%d/%Y %I:%M:%S %p'))
    except:
        print 'Error! iname: %s' % iname
    ilast = sorted_list[-1]

    # add lats/lons to lists
    llat.append(float(ilast['lat']))
    llon.append(float(ilast['lon']))
    names.append(int(ilast['name'].split('-')[-1]))

    # for statistics
    if ilast['height']<>'':
        height.append(float(ilast['height']))
    if ilast['length']<>'':
        length.append(float(ilast['length']))
    if ilast['width']<>'':
        width.append(float(ilast['width']))
    if ilast['direction']<>'':
        direction.append(float(ilast['direction']))
    if ilast['speed']<>'':
        speed.append(float(ilast['speed']))
    if ilast['type']<>'':
        type.append(ilast['type'])

    # to TXT file
    ff_last_copd_positions.write(('%s;%s;%s;%s;%s;%s;%s;%s;%s;%s;%s;%s;%s;\n') % (ilast['name'],ilast['type'],\
                                                 ilast['date'],ilast['lat'], \
                                                 ilast['lon'], ilast['direction'], \
                                                 ilast['speed'], ilast['height'], \
                                                 ilast['width'], ilast['length'], \
                                                 ilast['mass'], ilast['sources'], \
                                                 ilast['note']))
    # to SHAPE file
    w.point(float(ilast['lon']), float(ilast['lat']))
    w.record(ilast['name'],ilast['type'],ilast['date'],ilast['direction'],ilast['speed'],ilast['height'], \
             ilast['width'],ilast['length'],ilast['mass'],ilast['sources'],ilast['note'])

    # to model file
    if ilast['width']<>'' and ilast['length']<>'' and f_model_flag==1:
        if float(ilast['width'])>=10 and float(ilast['length'])>=10:
            f_model.write('%d %s %s %s %s 3 15\n' % (iid,ilast['lat'],ilast['lon'],ilast['width'],ilast['length']))
        else:
            # if less size than 10m
            f_model.write('%d %s %s 20 20 3 15\n' % (iid,ilast['lat'],ilast['lon']))
    else:
        if f_model_flag==1:
            f_model.write('%d %s %s 30 30 3 15\n' % (iid, ilast['lat'],ilast['lon']))
    # model date
    if f_model_flag==1:
        mod_date = datetime.datetime.strptime(ilast['date'],'%m/%d/%Y %I:%M:%S %p')
        date_model.write('%s %s %s %s\n' % (mod_date.year,mod_date.month,mod_date.day,mod_date.hour))

    iid = iid + 1
    # make an empty list
    ilist = []

try:
    f_model.close()
    date_model.close()
except:
    print 'Error!'

# create the PRJ file
prjfn = shp_filename
prj = open("%s.prj" % prjfn, "w")
epsg = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
prj.write(epsg)
prj.close()
w.save(shp_filename)

ff_last_copd_positions.close()

######################
# Statistics
######################

print '\nTotal: %s' % iid

w = np.array(width)
l = np.array(length)
h = np.array(height)
s = np.array(speed)
d = np.array(direction)

# width/length/heigth mean values
w_mean = w.mean()
l_mean = l.mean()
h_mean = h.mean()
s_mean = s.mean()

# Averaging azimuth
try:
    rads = [radians(i) for i in d]
    ave = lambda x: sum(x)/len(x)
    rotate = lambda x: x+2*(pi) if x<0 else x
    # Calculate the averages
    sa = ave([sin(i) for i in rads]) # Average of the sines
    ca = ave([cos(i) for i in rads]) # Average of the cosines
    # Take the arctan of the averages, rotate to positive values, and then
    # convert back to degrees.
    d_mean = degrees(rotate(atan2(sa,ca)))
except:
    print 'cannot calculate azimuth'

try:
    print '\nMean values\n Width: %5.2f m \n Length: %5.2f m \n Heigth: %5.2f m\
\n Speed: %5.2f kt\n Direction: %5.0f degrees\n' % \
      (w_mean, l_mean, h_mean, s_mean, d_mean)
except:
    pass

u_type = set(type)
# all targets number
all_targets_ch = 0

for itype in u_type:
    i = 0
    for iitype in type:
        if iitype == itype:
            i = i + 1
    print '%s = %d' % (itype, i)
    all_targets_ch = all_targets_ch + i

# checking total targets
print '\nTotal\All types: %s\%d\n' % (iid, all_targets_ch)

# closest targets
ilat = np.array(llat)
ilon = np.array(llon)
inames = np.array(names)

# West Alpha point
lat0 = 74.4561
lon0 = 63.8400

mdist = 10000.

fdist = open('txt//distances.txt','w')

for i in range(len(inames)):
    d, a = dist(lat0, lon0, ilat[i], ilon[i])
    dist_nm = nautical(meters = d)
    #print 'Nomer: %d %5.2f' % (i,dist_nm)
    if dist_nm < mdist:
        mdist = dist_nm
        aa = a
        print 'Min distance: %5.1f   %5.1f' % (mdist,a)
        idx = i
    # distance to a file
    fdist.write('%5d %5.1f %5.0f\n' % (names[i],dist_nm,a))

fdist.close()

print '\nClosest to (%s %s) is (%s %s) (#%s)' % (lat0, lon0, ilat[idx], ilon[idx], int(inames[idx]))
print 'Distance: %7.1f [nm]/%.0f' % (mdist,aa)