Ice targets Stat

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
from pylab import *

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

''' Create prj file '''
def create_prj(shp_filename):
    # 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)

''' Redefine type of an ice object based on length '''
def define_type(ll):
    if ll<5.:
        return 'Growler'
    if ll>=5. and ll<=15:
        return 'Bergy Bit'
    if ll>15. and ll<=60:
        return 'Small'
    if ll>60. and ll<=120:
        return 'Medium'
    if ll>120. and ll<=200:
        return 'Large'
    if ll>200. and ll<=300:
        return 'Very Large'
    if ll>300.:
        return 'Ice Island'

'''Statistics'''

def av_azimuth(direction):
    # Averaging azimuth
    try:
        d = np.array(direction)
        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)))
        return d_mean
    except:
        print 'cannot calculate azimuth'
        return 9999.

def statistics(width, length, height, speed, direction):
    print '\nTotal: %s' % total_number

    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 Higth: %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)
    #print u_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' % (total_number, all_targets_ch)

    # Statistics on ships
    u_ship = set(ship)
    print '*****************************'
    print '*** Statistics on vessels ***'
    print '*****************************'
    for iship in u_ship:
        i = 0
        for iiship in ship:
            if iiship == iship:
                i = i + 1
        print '%s = %d' % (iship, i)

def make_pie_plot(par, ititle, ilabels, ifracs, iexplode):
    plt.clf()
    # make a square figure and axes
    figure(1, figsize=(6,6))
    ax = axes([0.1, 0.1, 0.8, 0.8])

    # The slices will be ordered and plotted counter-clockwise.
    labels = ilabels
    fracs = ifracs
    explode = iexplode

    pie(fracs, explode=explode, labels=labels,
                autopct='%1.1f%%', shadow=True)
                # The default startangle is 0, which would start
                # the Frogs slice on the x-axis.  With startangle=90,
                # everything is rotated counter-clockwise by 90 degrees,
                # so the plotting starts on the positive y-axis.

    title(ititle, bbox={'facecolor':'0.8', 'pad':5})

    savefig('pics\\%s.png' % par)

############################
### Main Program
############################
source = 'vessel_targets'
ff = sys.argv[1]
width = []
length = []
height = []
speed = []
direction = []
type = []
ship = []
temp_width = []
temp_length = []
temp_hight = []
temp_speed = []
temp_direction = []
temp_type = []
all_target_statistics = []

###########################
#### Shapefile creation
###########################
shp_filename = 'shp\\%s' % source # ff.split('.')[0] #'updated_targets.shp'
# create the shapefile
w = sf.Writer(sf.POLYLINE)

w.field('name','C','90')
w.field('type','C','90')
w.field('start','C','90')
w.field('end','C','90')
w.field('number_obs','C','90')
w.field('av_width','C','90')
w.field('av_length','C','90')
w.field('av_height','C','90')
w.field('av_speed','C','90')
w.field('azimuth','C','90')
############################

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

unique_names = set(names)

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

lll = []

# Total number of tracks
total_number = 0

# Drift rose file
ff_rose = open('txt\\speed_dir.txt','w')
# File within numbers of observations
num_obs_ll = []
ff_obs = open('txt\\numbers_obs.csv','w')
ff_obs.write('ID,Observations,\n')
# Size statistics
ff_sizes = open('txt\\vessel_target_av_size.csv','w')
ff_sizes.write('ID,Type,Start Date,End Date,Num of obs,Width,Length,Height,Speed,Direction,\n')

for iname in unique_names:
    for rownum in range(sheet.nrows):
        row = sheet.row_values(rownum)
        if row[0]==iname and row[13]=="Vessel":
            # List withing all data for a current target
            if (float(row[6])>50. and float(row[6])<90. and float(row[5])>69. and float(row[5])<78.):
                ilist.append({'name':row[0], 'type':row[2], 'date':row[4], \
                      'lat':row[5], 'lon':row[6], 'direction':row[7], \
                      'speed':row[8], 'hight':row[9], 'width':row[10], \
                      'length':row[11], 'mass':row[12], 'sources':row[13], \
                      'note':row[14], 'entered':row[15]})
    try:
        if len(ilist)>0:
            #print '%s %s' % (iname, len(ilist))
            # List within number of observations
            num_obs_ll.append([iname, len(ilist)])
            # Total number of tracks
            total_number = total_number + 1
            #print iname
            sorted_list = sorted(ilist, key=lambda x: datetime.datetime.strptime(x['date'],\
                                                                    '%m/%d/%Y %I:%M:%S %p'))

            # Start and End date of measaurments
            start_date = sorted_list[0]['date']
            end_date = sorted_list[-1]['date']
            #print start_date, end_date

            # Looking through elements of a list and retrieve parameters
            for key in sorted_list:
                # lat lon
                lll.append([float(key['lon']),float(key['lat'])])
                # Sizes and types
                try:
                    # Lists for overall and current statistics
                    if key['width']<>'':
                        width.append(float(key['width']))
                        temp_width.append(float(key['width']))
                    if key['length']<>'':
                        length.append(float(key['length']))
                        temp_length.append(float(key['length']))
                    if key['hight']<>'':
                        height.append(float(key['hight']))
                        temp_hight.append(float(key['hight']))
                    if key['speed']<>'':
                        if float(key['speed'])<3.:
                            speed.append(float(key['speed']))
                            temp_speed.append(float(key['speed']))
                            # Speed and direction to a file
                            ff_rose.write('%s %s\n' % (float(key['speed']),float(key['direction'])))
                    if key['direction']<>'':
                        direction.append(float(key['direction']))
                        temp_direction.append(float(key['direction']))

                    if key['entered']<>'':
                        ship.append(key['entered'])

                     # Redefine type
                    if key['length']<>'':
                        ll = float(key['length'])
                        new_type = define_type(ll)
                        type.append(new_type)
                        temp_type.append(new_type)
                    else:
                        type.append(key['type'])
                        temp_type.append(key['type'])
                except:
                    pass


            # Filtering wrong coordinates by distance
            # Firs values
            lon_prev = float(lll[0][0])
            lat_prev = float(lll[0][1])
            #print "%s lon: %s lat: %s" % (iname, lon_prev, lat_prev)
            ch = 0
            for ii in lll:
                # Check the distance between coordinates
                #print ii
                llat2  = float(ii[1])
                llong2 = float(ii[0])
                try:
                    d, a = dist(lat_prev, lon_prev, llat2,llong2)
                    dist_nm = nautical(meters = d)
                    # If distance > 30 then delete an element
                    if dist_nm>30:
                        #print "%.2f [nm]" % dist_nm
                        del lll[ch]
                except:

                    dist_nm = 0.
                ch = ch + 1

            # Current statistics on each track
            w_mean = np.array(temp_width).mean()
            l_mean = np.array(temp_length).mean()
            h_mean = np.array(temp_hight).mean()
            s_mean = np.array(temp_speed).mean()
            d_mean = av_azimuth(temp_direction)

            if len(temp_type)==0:
                temp_type='Unknown'
            else:
                temp_type = temp_type[0]

            if np.isnan(w_mean)==True:
                w_mean = 9999.
            if np.isnan(l_mean)==True:
                l_mean = 9999.
            if np.isnan(h_mean)==True:
                h_mean = 9999.
            if np.isnan(s_mean)==True:
                s_mean = 9999.
            if np.isnan(d_mean)==True:
                d_mean = 9999.

            # Add al data to a dictonary
            all_target_statistics.append({'name':iname, 'type':temp_type, 'date1':start_date, \
                                    'date2':end_date, 'num_obs': len(sorted_list), \
                                    'width':w_mean, 'length':l_mean, 'height':h_mean, \
                                    'speed':s_mean, 'direction':d_mean})

            temp_width = []
            temp_length = []
            temp_hight = []
            temp_speed = []
            temp_direction = []


    except:
        print 'Error! iname: %s' % iname
    try:
        ilast = sorted_list[-1]
    except:
        pass

    #############################
    ### Shapefile creation
    #############################
    try:
        if len(lll)>0:
            w.line(parts=[lll])
            w.record('%s' % iname,'%s' % temp_type,'%s' % start_date, '%s' % end_date, str(len(sorted_list)),\
                     str('%5.1f' % w_mean),str('%5.1f' % l_mean), str('%5.1f' % h_mean), str('%5.1f' % s_mean), str('%5.1f' % d_mean))

    except:
        pass
    temp_type = []
    lll = []
    # make an empty list
    ilist = []

try:
    all_targets_sorted_list = sorted(all_target_statistics, key=lambda x: (x['name']))
    for i in range(0,len(all_targets_sorted_list)):
        ff_sizes.write('%s,%s,%s,%s,%d,%5.1f,%5.1f,%5.1f,%5.1f,%5.0f\n' % (all_targets_sorted_list[i]['name'],all_targets_sorted_list[i]['type'],all_targets_sorted_list[i]['date1'],all_targets_sorted_list[i]['date2'],\
                                                                                       all_targets_sorted_list[i]['num_obs'],all_targets_sorted_list[i]['width'],all_targets_sorted_list[i]['length'],\
                                                                                       all_targets_sorted_list[i]['height'],all_targets_sorted_list[i]['speed'],all_targets_sorted_list[i]['direction']))
except:
    pass


# Get statistics
try:
    statistics(width, length, height, speed, direction)
except:
    "Statistics error"

create_prj(shp_filename)

# Plots
# Size Distribution
ilabels = 'Growler', 'Bergy Bit', 'Small'
ifracs = [27, 49, 91]
iexplode=(0, 0.05, 0)
par = 'Types'
ititle = 'Size Distribution of Ice Targets'
make_pie_plot(par, ititle, ilabels, ifracs, iexplode)

# Close txt files
ff_rose.close()
sorted_num_obs_ll = sorted(num_obs_ll, key=lambda x: x[1], reverse=True)
for iid,inum in sorted_num_obs_ll:
    ff_obs.write('%s,%s\n' % (iid,inum))

ff_obs.close()
ff_sizes.close()