#!/usr/bin/python

# Reformat the Global CMT Project file in NDK format to a GEON IDV ascii "csv" point data file.

# First get the ASCII file of Global CMT Project data, 
# from the CMT catalog at http://www.globalcmt.org/CMTfiles.html
# The CMT catalog is available in the ASCII "NDK" format, called 'jan76_dec10.ndk'
# and has all GMT CMT solutions for  1976-2010; more recent data may be found there, by month.
 
# This program creates a csv ASCII file for IDV use of selected values per event, 
# for focal mechanisms and other data.

# on Linux, run this program with a command line like this: 
 
# GMT_ndk_file_to_csv_with_subsetting.py  howMany mbMin mbMax MsMin MsMax depthMin depthMax latMin latMax lonMin lonMax
#                                   0        1      2    3     4     5      6        7       8       9     10     11       

# GMT_ndk_file_to_csv_with_subsetting.py 500 0.0 9.0 0.0 9.0 30 800 -50.0 10.0 -90 -55
 
# depthMax is largest POSITIVE; depthMin also POSITIVE kms below surface 
#   (depth value are converted later to negative values for IDV use)
# The valid input longitude range is -180.0 (west) to +180.0 (east); same as data file uses.
 
# created Sep. 11-12 2008
#  SKW, UNAVCO 
# revisions 
# Oct 27 2008 SKW
# Aug 16 2012 SKW

import os
import sys
import string
import math

# main program
def main():
    extract()
# end of main 


def extract ():

        # get and repair command line argument input values
        howMany = 0;
        value=int(sys.argv[1])
        if (value>0):
           howMany=value
           print "  how many is "+`howMany`
        mbMin = 0;
        value=float(sys.argv[2])
        if (value>0):
           mbMin=value
           print "   mbMin is "+`mbMin`
        mbMax = 0;
        value=float(sys.argv[3])
        if (value>0):
           mbMax=value
           print "   mbMax is "+`mbMax`
        if (mbMax < mbMin):
           hold=mbMin
           mbMin=mbMax
           mbMax=hold
        print "   mb range is "+`mbMin`+" to "+`mbMax`
        MsMin = 0;
        value=float(sys.argv[4])
        if (value>0):
           MsMin=value
           #print "   MsMin is "+`MsMin`
        MsMax = 0;
        value=float(sys.argv[5])
        if (value>0):
           MsMax=value
           #print "   MsMax is "+`MsMax`
        if (MsMax < MsMin):
           hold=MsMin
           MsMin=MsMax
           MsMax=hold
        print "   Ms range is "+`MsMin`+" to "+`MsMax`
        depthMin =  float(sys.argv[6]) #positive-down coordinates for depth ; same as data in input file
        depthMax =  float(sys.argv[7]) #positive-down coordinates for depth ; same as data in input file
        latMin = float(sys.argv[8]) 
        latMax = float(sys.argv[9]) 
        lonMin = float(sys.argv[10]) 
        lonMax = float(sys.argv[11])
        print "   lat min max;  long min max "+`latMin`+" "+`latMax`+";    "+`lonMin`+"  "+`lonMax`

        # open the input file 
        #          jan76_dec10.ndk
        filename= 'jan76_dec10.ndk'
        # the file handler is
        file1 = open (filename)

        # read and count how many lines in file, five per quake 
        allLines = file1.readlines()
        linecount = len(allLines)
        ct = linecount
        #print " count of lines  in input file is "+`ct`
        nevs = ct / 5.0
        nloop  = int(nevs)
        print "  there are "+`nloop`+ " earthquakes in the input file "+filename

        file1.seek(0) # rewind to beginning

        # name and open the output file
        mb1= (`mbMin`)[0:4]
        mb2= (`mbMax`)[0:4]
        ms1= (`MsMin`)[0:4]
        ms2= (`MsMax`)[0:4]
        d1=(sys.argv[6])
        d2=(sys.argv[7])
        lt1=(sys.argv[8])
        lt2=(sys.argv[9])
        lo1=(sys.argv[10])
        lo2=(sys.argv[11])
        #ofname = 'GMT_focal_mechanisms_'+`howMany`+"_"+mb1+"_"+mb2+"_"+ms1+"_"+ms2+"_"+d1+"_"+d2+"_"+lt1+"_"+lt2+"_"+lo1+"_"+lo2+'.csv'
        ofname = 'GMT_focal_mechanisms_'+mb1+"_"+mb2+"_"+ms1+"_"+ms2+"_"+d1+"_"+d2+"_"+lt1+"_"+lt2+"_"+lo1+"_"+lo2+'.csv'
        #ofname = 'GMT_focal_mechanisms_SouthAmerica.csv'
        csvfile = open (ofname, 'w')
        print "  output filename is "+ofname


        # write the two header lines used and required by the IDV for an IDV "Text Point Data File"
        # which describe the data charater in this case

        # the IDV recognizes and must have special variable names for focal mechanisms:
        # float strike  // degrees CW from North   
        # float dip     // degrees downwards from horizon level
        # float rake    // degrees down from horizon, in plane of mechanism, of motion

        csvfile.write("(index) -> (Time,Latitude,Longitude,Altitude,mb,Ms,magnitude,strike,dip,rake)")
        csvfile.write("\n")

        csvfile.write("Time[fmt=\"yyyy/MM/dd HH:mm:ss\"],Latitude[unit=\"deg\"],Longitude[unit=\"deg\"],Altitude[unit=\"km\"],mb[unit=\"null\"],Ms[unit=\"null\"],magnitude[unit=\"null\"],strike[unit=\"null\"],dip[unit=\"null\"],rake[unit=\"null\"]")
        csvfile.write("\n")

        goodcount = 0
        counthowMany=0;


        # each event is a five line set like this
        #0    5    10   15   20   25   30   35   40   45   50   55
        #0    5    .    .    .    .    .    .    .    .    .    .
        #PDE  2005/12/31 12:14:02.2 -28.99  -71.52  25.0 5.1 0.0 NEAR COAST OF CENTRAL CH
        #C200512311214A   B: 12   13  40 S: 40   63  50 M:  0    0   0 CMT: 1 TRIHD:  0.8
        #CENTROID:      4.9 0.5 -29.21 0.03  -71.79 0.05  24.7  2.3 FREE S-20060404111308
        #23  0.071 0.192  1.300 0.130 -1.370 0.151 -2.250 0.411 -3.240 0.441  0.009 0.098
        #V10   3.956 44 144   0.390 18  36  -4.346 41 290   4.151 311 18    5 217 88  108
        # 0     1    2   3    4     5   6    7     8   9    10    11  12   13  14 15   16


        # PDE  2005/12/28 12:21:58.7 -19.86  -68.92 103.6 5.0 0.0 CHILE-BOLIVIA BORDER REG

        # read each event 
        for i in range(nloop) :
    
             mb = "0.0";
             MS = "0.0";
             magnitude= "0.0";

             # read line 1
             line = file1.readline()
             # print line  
             split = string.split(line) # split on whitespace
             date  = line[5:15]  # split[1]
             time1 = line [16:24] # split[2]
             lat   = line [26:33] #split[3]
             dv = string.atof(lat)
             #print " lat is  "+`dv`
             # check for  range, if wanted; skip to next file if fails this test:
             if ( (latMin!=0 or latMax!=0) and (dv<latMin or dv>latMax) ) :
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 continue # go try next quake
             #else :
             #    print " lat is good : "+`dv`
             #print " latitude "+lat+" is ok "
             long  = line[34:41] #split[4]
             dv = string.atof(long)
             # check for  range, if wanted; skip to next file if fails this test:
             if ( (lonMin!=0 or lonMax!=0) and (dv<lonMin or dv>lonMax) ) :
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 continue # go try next quake
             #print " longitude "+long+" is ok "

             deep  = line[42:47] # split[5]
             dv = string.atof(deep)
             # check for depth range if wanted; skip to next file if fails this test:
             if ((depthMin!=0 or depthMax!=0) and (dv==0 or (dv<depthMin or dv>depthMax)) ) :
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 continue # go try next quake
             #print "    depth "+deep+" is ok "
             dv  =  -1.0 * string.atof(deep)  # IDV uses negative downwards depths
             deep = ""+`dv`

             mb = line [48:51]
             MS = line [52:55]

             #print `goodcount` + "   " +line
             #print "         line has mb = "+mb +"   Ms = "+MS

             # check magnitude ranges
             mv = string.atof(mb)
             if (mbMax>0 and (mv<mbMin or mv>mbMax)) :
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 continue # go try next quake
             #print "       mb "+mb+" is in range"
             mv = string.atof(MS)
             if (MsMax>0 and (mv<MsMin or mv>MsMax) ):
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 line = file1.readline()
                 continue # go try next quake
             #print "          ms "+MS+" is in range"

             timestr = date+" "+time1 

             # read line 2
             line = file1.readline()
             # read line 3
             line = file1.readline()
             # read line 4
             line = file1.readline()

             # read line 5
             line = file1.readline()
             split = string.split(line) # split on whitespace
             # 1st modal plane values:
             #strike  = split[11]
             #dip  = split[12]
             #rake  = split[13]
             # 2nd modal plane values:
             strike  = split[14]
             dip  =    split[15]
             rake  =   split[16]

             mbv = string.atof(mb)
             MSv = string.atof(MS)
             # set generic "magnitude" variable vlaue used by IDV to color dots in some cases
             if (mbv > 0.0) :
                magnitude = mb
             if (MSv>0.0 and MSv> mbv) :
                magnitude = MS 

             # if (i<10): 
             #   print "   date is "+date+"  time1 is "+time1+"  time all is "+timestr
             #   print timestr +"," +lat+","+long+","+ deep+","+mb+","+MS+","+magnitude+","+strike+","+dip+","+rake 

             # flip-o
             hold=dip
             dip =rake
             rake=hold
             csvfile.write(timestr+"," +lat+","+long+","+ deep+","+mb+","+MS+","+magnitude+","+strike+","+dip+","+rake)

             csvfile.write("\n")

             goodcount +=1
             counthowMany += 1;

             #print "      "+timestr +"," +lat+","+long+","+ deep+","+mb+","+MS+","+magnitude+","+strike+","+dip+","+rake 
             #print "   found "+`counthowMany`+" matching quakes "
             if (howMany>0 and counthowMany >= howMany):
                #print "  made file with "+ `howMany` +" earthquakes "
                break
 
             # end for loop on each event  

        csvfile.close()

        print "     processed "+ `goodcount` +" earthquakes"
        print "   found "+`counthowMany`+" matching quakes "
        print "  output file is "+ofname

# END OF MODULE extract


#call main
main()

# ALL DONE