source: misc/tools/event_selection/find_faullt/find_fault.py @ 7378

""" Created by: Jonathan Griffin
    Date: 29 July 2009
    Purpose: This script searches output from Ross Wilson's event selector to find events that are within a geographic area.
    Maximum and Minimum longitude and latitude values are given, and a search distance within which the max/min lon/lat of
    the quake must fall within to be selected.
"""


search_dist = 0.75 # radius in degrees to search 

# List coordinates to search for
MaxLon = 120.2605
MinLon = 113.7402
MaxLat = -9.4653
MinLat = -11.2732

# Specify input Folder
inDir = r'/nas/gemd/georisk_models/inundation/data/western_australia/broome_tsunami_scenario_2009/anuga/boundaries/Results_Australia_1629_0.01_15.00/'
inFile = 'fault.xy'
#inFile = 'fault.txt'
inFilePath = file(inDir + inFile)

# Output file to list I.D. of possible quakes
outFile = 'possible_quakes.txt'
outFilePath = file(inDir + outFile, 'w')

#Establish variable and lists
quake = 0
possible_quakes = []
quake_list = []
count = 0
# Read data from fault.xy 
for line in inFilePath.readlines():
    #Skip header row
    if "Lon" in line:
        continue
    #Split data and read into variables
    data = line.split(',')
    Lon = float(data[0])
    Lat = float(data[1])
    Quake_ID = int(data[2])
    subfault_ID = data[3].strip('\r\n')
    subfault_ID = int(subfault_ID)

    # if new quake_ID, find latitude/longitude range and determine if
    # within search distance
    if Quake_ID != quake:
        #print 'quake list', quake_list
        if len(quake_list) != 0:
            lonlist = []
            latlist = []
            for i in range(len(quake_list)):
                lonlist.append(quake_list[i][0])
                latlist.append(quake_list[i][1])
                #print 'lonlist', lonlist
            maxlon = max(lonlist) # find max/min lat/lon here
            minlon = min(lonlist)
            maxlat = max(latlist)
            minlat = min(latlist)
            if ((abs(maxlon-MaxLon) < search_dist) & (abs(minlon-MinLon) < search_dist) & (abs(maxlat-MaxLat) < search_dist) & (abs(minlat-MinLat) < search_dist)):
                possible_quakes.append(Quake_ID)
        quake = Quake_ID
        quake_list = [] #reset quake_list for next earthquake
        quake_list.append([Lon,Lat,Quake_ID,subfault_ID])
        count += 1
    else:   #append data to list of subfault coordinates for that quake
        quake_list.append([Lon,Lat,Quake_ID,subfault_ID])

# Print out all earthqaukes that fall within search area.
print 'possible quakes', possible_quakes, 'out of ',count, ' earthquakes'

#write out possible earthquakes to text file
outFilePath.write('Possible earthquakes\n')
for quakes in possible_quakes:
    print quakes
    outFilePath.write(str(quakes) + '\n')


inFilePath.close()
outFilePath.close()