import os import sys import shutil from datetime import datetime # Author: Alistair Knox # Based on Emilie Klein's Predict_weeklyDispl_2ndStep.com, recorded speedup 360x+ # # This version sets a t0 at an arbitrary point in the data. Epochs before will have a # negative difference, and after will have the same difference we normally calculate. # GRID_FOLDER = 'displGrid/' INPUT_FOLDER = 'displGrid/DataST/ST/' OUTPUT_FOLDER = 'displGrid/DataST/TablesDisp/' # Additional output folder for data predicted by SECTOR PRED_INPUT_FOLDER = 'displGrid/DataST/ST_predicted/' PRED_OUTPUT_FOLDER = 'displGrid/DataST/TablesDisp_predicted/' # Dates must be set to an exact epoch in downsampled file, otherwise 0.0192 increment will not match # 4 digit decimal places necessary. #START_DATE = 1999.9968 START_DATE = 2010.0000 T0_DATE = 2010.0000 TMAX = 2970 def main(): station_list, station_coords = get_stations() iter_stations(INPUT_FOLDER, OUTPUT_FOLDER, station_list, station_coords) #iter_stations(PRED_INPUT_FOLDER, PRED_OUTPUT_FOLDER, station_list, station_coords) def iter_stations(input_folder, output_folder, station_list, station_coords): """Check the existence of files within the input folder and append data to files in the output folder. Args: input_folder (str): Input folder of station.dat files. output_folder (str): Output folder of station.dat files. station_list (list): List of active stations processing data. station_coords (dict): Station : (Lat, lon) dictionary. """ print("**** Iterating stations for input folder: " + input_folder + " ****") # clear output directory if os.path.isdir(output_folder): shutil.rmtree(output_folder) os.makedirs(output_folder) # iterate our stations for idx, station in enumerate(station_list): station_file = input_folder + station + '.dat' station_file_non_acc = input_folder + 'nonAcc_' + station + '.dat' if os.path.isfile(station_file): return_value = append_data(station, station_file, output_folder, station_coords, date_file_prefix='accumulated_displacements_') if return_value: print("### Predict2 done for station " + station + " (" + str(idx+1) + "/" + str(len(station_list)) + "), results stored ###") # non-accumulating data depends on accumulating data if os.path.isfile(station_file_non_acc): append_data(station, station_file_non_acc, output_folder, station_coords, date_file_prefix='non_accumulating_') else: print("Skipped non-accumulating data for station " + station + " - not found in " + input_folder) else: print("*** Skipped station " + station + " - could not find t0: " + str(T0_DATE) + " ***") else: print("--- Skipped station " + station + " - not found in " + input_folder + " ---") def append_data(station, station_file, output_folder, coords, date_file_prefix): """Append the station data to each date file (1 file per station -> 1 file per epoch) Args: station (str): 4-char station code. station_file (str): Data filename containing either accumulating or non-accumulating values. output_folder (str): Filename for destination containing data organized by epoch. coords (tuple): (lon, lat) date_file_prefix (str): Prefix for output filename. Returns: bool: Whether the operation was successful. Fails if station data ends before the preset t0. """ global START_DATE, T0_DATE, TMAX date = START_DATE # Find t0 in the file to set dn0, de0, du0 # then read again from first epoch to write negative diffs until t0 t0_reached = False with open(station_file, 'r') as sdata: # first block to run - get dn0, de0, du0 initialized if not t0_reached: for line in sdata: find_t0_line_items = [float(i) for i in line.split()] # first block to run if T0_DATE == find_t0_line_items[0]: # initial epoch data dn0 = float("{:7.2f}".format(find_t0_line_items[1])) de0 = float("{:7.2f}".format(find_t0_line_items[2])) du0 = float("{:7.2f}".format(find_t0_line_items[3])) t0_reached = True # read file again if t0_reached: sdata.seek(0) for line in sdata: line_items = [float(i) for i in line.split()] # increment date # if file is missing epochs, keep incrementing to catch up while date < line_items[0] and date <= TMAX: date += 0.0192 date = float("{:.4f}".format(date)) if date == line_items[0]: # append this station's data and differences to each date file date_file = output_folder + date_file_prefix + "{:.4f}".format(date) + '_mm.dat' with open(date_file, 'a') as df: df.write(station+' ') df.write(coords[station][0] + ' ' + coords[station][1] + ' ') dn, de, du = line_items[1], line_items[2], line_items[3] sn, se, su = line_items[4], line_items[5], line_items[6] diffn, diffe, diffu = dn - dn0, de - de0, du - du0 # print(station_file,line_items[0],dn,de,du,sn,se,su,diffn,diffe,diffu) df.write("{:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f} {:5.2f}\n".format(diffn, diffe, diffu, sn, se, su)) return t0_reached def get_stations(): """Get the active stations and coordinates. Returns: (tuple): (List of stations, Map of station : coordinates) """ TSLen = get_active_station_file() if not TSLen: print("Couldn't find the active station list") sys.exit() station_list = [] with open(GRID_FOLDER + TSLen, 'r') as ts: station_coords = {} for line in ts: line_items = line.split() station_list.append(line_items[0]) # station_coords[line_items[0]] = "{:.5f}".format(float(line_items[1])), "{:.5f}".format(float(line_items[2])) station_coords[line_items[0]] = "{:.9f}".format(float(line_items[1])), "{:.9f}".format(float(line_items[2])) return station_list, station_coords def get_active_station_file(): """Get the latest file with stations used for gridding. Returns: str: Filename of 3-column station file with Site lat lon """ for f in sorted(os.listdir(GRID_FOLDER), reverse=True): if f.startswith('Station_Locations_'): print('got active: ', f) return f if __name__ == '__main__': s = datetime.now() main() print("Duration: " + str(datetime.now() - s)[:-4])