diff --git a/cz.py b/cz.py index db8aa13..9934e47 100644 --- a/cz.py +++ b/cz.py @@ -4,7 +4,7 @@ warnings.filterwarnings('ignore') import argparse, json from common import NoIndent, MyEncoder import os, sys -sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))+'/pyat') +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + '/pyat') def _argparse(): parser = argparse.ArgumentParser(description='Average Noise Spectral Level') @@ -18,13 +18,13 @@ def main(): if argfile[-5:] != '.json': print('输入的参数文件名有误') return - with open(argfile,'r') as jsonFile: + with open(argfile, 'r') as jsonFile: paramData = json.load(jsonFile) if ("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or \ - ("freq" not in paramData) or ("sourceDepth" not in paramData) or ("receiveDepth" not in paramData) or\ - ("lat" not in paramData) or ("lon" not in paramData) or\ - ("timeType" not in paramData) or\ - ("outParam" not in paramData): + ("freq" not in paramData) or ("sourceDepth" not in paramData) or ("receiveDepth" not in paramData) or \ + ("lat" not in paramData) or ("lon" not in paramData) or \ + ("timeType" not in paramData) or \ + ("outParam" not in paramData): print(f'{argfile}中的参数有误，请核实！') return ETOPO1Path = paramData["ETOPO1Path"] @@ -42,7 +42,7 @@ def main(): outParam = paramData["outParam"] if "outPath" not in outParam: outParam["outPath"] = "" - + from DataBase import DataBase from czcore import call_cz import numpy as np @@ -50,7 +50,7 @@ def main(): # 定义时间段 timeIndex = timeType # 创建类 - database = DataBase(ETOPO1Path, STPath,dataSource="") + database = DataBase(ETOPO1Path, STPath, dataSource="") # 设置时间段 database.set_timePeriod(timeIndex) @@ -60,11 +60,11 @@ def main(): # czdata.TL是损失、czdata.bathym是地形数据 # czdata.Rays是声线数据 # czdata.CzDetail是会聚区数据，包括会聚区宽度和距离 - res = call_cz(database, freq, lat, lon, azi, sourceDepth, receiveDepth) + res = call_cz(database, freq, lat, lon, azi, sourceDepth) # call_cz(25, 18.0, 114.0, 10, 3500) # write data to file - #输出json文件 + # 输出json文件 outData = {} outData["param"] = paramData if len(res) == 0: @@ -72,10 +72,10 @@ def main(): outData["data"] = [] else: outData["result"] = "有会聚区" - outDataData, outBathymData,outTlData,outRayData = [],[],[],[] + outDataData, outBathymData, outTlData, outRayData = [], [], [], [] json.encoder.FLOAT_REPR = lambda x: format(x, '.3f') for index in range(len(res)): - bathymData,tlData,rayData = {},{},{} + bathymData, tlData, rayData = {}, {}, {} # 保存地形数据 res[index].bathym.r = np.around(res[index].bathym.r, decimals=0) res[index].bathym.d = np.around(res[index].bathym.d, decimals=0) @@ -84,18 +84,20 @@ def main(): "azi": f'{round(res[index].azi)}', "count": len(res[index].bathym.d), "descript": "距离、深度", - "data": NoIndent(list(zip(res[index].bathym.r,res[index].bathym.d))) + "data": NoIndent(list(zip(res[index].bathym.r, res[index].bathym.d))) } # 保存传播损失数据 tlpoint = [] res[index].depths = np.around(res[index].depths, decimals=0) res[index].ranges = np.around(res[index].ranges, decimals=0) res[index].TL = np.around(res[index].TL, decimals=2) - tlpoint = list(zip(np.tile(res[index].depths.reshape(len(res[index].depths),1),len(res[index].ranges)).flatten(), np.tile(res[index].ranges,(len(res[index].depths),1)).flatten(), res[index].TL.flatten())) + tlpoint = list( + zip(np.tile(res[index].depths.reshape(len(res[index].depths), 1), len(res[index].ranges)).flatten(), + np.tile(res[index].ranges, (len(res[index].depths), 1)).flatten(), res[index].TL.flatten())) tlData = { "param": paramData, "azi": f'{round(res[index].azi)}', - "count": len(res[index].ranges)*len(res[index].depths), + "count": len(res[index].ranges) * len(res[index].depths), "descript": "距离、深度、损失", "data": NoIndent(tlpoint) } @@ -106,7 +108,7 @@ def main(): ray.xy[1] = np.around(ray.xy[1], decimals=2) tempdata = { "pointCount": len(ray.xy[0]), - "data": NoIndent(list(zip(ray.xy[0],ray.xy[1]))) + "data": NoIndent(list(zip(ray.xy[0], ray.xy[1]))) } raypoint.append(tempdata) rayData = { @@ -124,9 +126,9 @@ def main(): directczdetaillist = [] for directcz in res[index].CzDetail.directCz: directczdetail = { - "minRange": f'{round(directcz.minRange,1)}', - "maxRange": f'{round(directcz.maxRange,1)}', - "width": f'{round(directcz.width,1)}', + "minRange": f'{round(directcz.minRange, 1)}', + "maxRange": f'{round(directcz.maxRange, 1)}', + "width": f'{round(directcz.width, 1)}', } directczdetaillist.append(directczdetail) czDeatilData["directCz"] = directczdetaillist @@ -136,9 +138,9 @@ def main(): reflectczdetaillist = [] for reflectcz in res[index].CzDetail.reflectCz: reflectczdetail = { - "minRange": f'{round(reflectcz.minRange,1)}', - "maxRange": f'{round(reflectcz.maxRange,1)}', - "width": f'{round(reflectcz.width,1)}', + "minRange": f'{round(reflectcz.minRange, 1)}', + "maxRange": f'{round(reflectcz.maxRange, 1)}', + "width": f'{round(reflectcz.width, 1)}', } reflectczdetaillist.append(reflectczdetail) czDeatilData["reflectCz"] = reflectczdetaillist @@ -146,19 +148,26 @@ def main(): outBathymData.append(bathymData) outTlData.append(tlData) outRayData.append(rayData) - with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"czOut_bathym.json", 'w') as jsonFile: + with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "czOut_bathym.json", + 'w') as jsonFile: json_data = json.dumps(outBathymData, cls=MyEncoder, indent=2, ensure_ascii=False) jsonFile.write(json_data) - with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"czOut_tl.json", 'w') as jsonFile: + with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "czOut_tl.json", + 'w') as jsonFile: json_data = json.dumps(outTlData, cls=MyEncoder, indent=2, ensure_ascii=False) jsonFile.write(json_data) - with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"czOut_ray.json", 'w') as jsonFile: + with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "czOut_ray.json", + 'w') as jsonFile: json_data = json.dumps(outRayData, cls=MyEncoder, indent=2, ensure_ascii=False) jsonFile.write(json_data) outData["data"] = outDataData - - with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"czOut.json", 'w') as jsonFile: + + with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "czOut.json", 'w') as jsonFile: json.dump(outData, jsonFile, indent=4, ensure_ascii=False) - + +# import time if __name__ == '__main__': - main() \ No newline at end of file + # start = time.perf_counter() + main() + # end = time.perf_counter() + # print('time = ',(end-start)) \ No newline at end of file diff --git a/pyat/bellhop.py b/pyat/bellhop.py index d2628ee..d957253 100644 --- a/pyat/bellhop.py +++ b/pyat/bellhop.py @@ -7,6 +7,7 @@ from multiprocessing import Pool from scipy.stats import norm import math + def write_ssp(sspfile, ssp, bathm, NZmax): if sspfile[-3:] != 'ssp': sspfile += '.ssp' @@ -35,9 +36,7 @@ def write_bathy(btyfile, bathm): f.write(' ' + str(bathm.r[row]) + ' ' + str(bathm.d[row]) + '\r\n') return -def call_Bellhop_p(filename): - system(AtBinPath + "/bellhop " + filename) - + def calZmax(ssp): d = np.array([], dtype=int) for it in range(len(ssp)): @@ -47,15 +46,18 @@ def calZmax(ssp): Zmax = ssp[ssp_idx].z[-1] return NZmax, Zmax, ssp_idx -def call_Bellhop(freq, sd, rd, bathm, ssp, sed, base): +def call_Bellhop_p(filename): + system(AtBinPath + "/bellhop " + filename) + +def call_Bellhop(freq, sd, bathm, ssp, sed, base): # Source and receiving position filename = 'envB' Rmax = bathm.r[-1] # unit: km - dr = (1500 / freq) / 1000 # wavelength / 10, unit: km + dr = 0.1 # wavelength / 10, unit: km ran = np.linspace(0, Rmax, round(Rmax / dr)) # Receiving range (unit: km) NZmax, Zmax, ssp_idx = calZmax(ssp) RD = np.linspace(0, Zmax, round(Zmax / 10) + 1) - pos = Pos(Source(sd), Dom(ran, rd)) + pos = Pos(Source(sd), Dom(ran, RD)) # Range of phase velocity cInt.Low = 1400 cInt.High = 15000 @@ -100,18 +102,22 @@ def call_Bellhop(freq, sd, rd, bathm, ssp, sed, base): bdy = Bndry(top, bottom) # Beam params run_type = 'C' # incoherently sum beams, see AT docs for more info + NAlphaRange = 12 - Alpha = alphadiv(NAlphaRange,Rmax) # min and max launch angle -20 degrees to 20 degrees - totalBeams = beamsnumber(freq, Rmax, max(bathm.r)) + Alpha = alphadiv(NAlphaRange,Rmax) # min and max launch angle -20 degrees to 20 degrees + totalBeams = beamsnumber(freq,Rmax,max(bathm.r)) box = Box(Zmax, max(bathm.r)) # bound the region you let the beams go, depth in meters and range in km deltas = 0 # length step of ray trace, 0 means automatically choose + + # Write *.env file + ialpha = 0 Filenames = [] - for iAlphaRange in range(len(Alpha) - 1): + for iAlphaRange in range(len(Alpha)-1): alpha = np.array([Alpha[iAlphaRange], Alpha[iAlphaRange + 1]]) - nbeams = int(totalBeams * (alpha[1] - alpha[0]) / 180) + nbeams = int(totalBeams * (alpha[1]-alpha[0])/180) beam = Beam(RunType=run_type, Nbeams=nbeams, alpha=alpha, box=box, deltas=deltas) # package - filenameI = filename + str(iAlphaRange) + filenameI = filename+ str(iAlphaRange) write_env(filenameI + '.env', 'BELLHOP', 'Pekeris profile', freq, sspB, bdy, pos, beam, cInt, Rmax) write_ssp(filenameI, ssp, bathm, NZmax) write_bathy(filenameI, bathm) @@ -122,7 +128,7 @@ def call_Bellhop(freq, sd, rd, bathm, ssp, sed, base): # Read sound field [x, x, x, x, Pos1, pressure] = read_shd(filename + '0.shd') - for iAlphaRange in range(NAlphaRange - 1): + for iAlphaRange in range(NAlphaRange-1): [x, x, x, x, Pos1, pressure1] = read_shd(filename + str(iAlphaRange) + '.shd') pressure = pressure + pressure1 TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss @@ -212,61 +218,63 @@ def find_cvgcRays(rays_total): Rays.append(ray) return Rays -def alphadiv(NalphaRange, Rmax): - angle = [-90, 90] - if Rmax > 100 and Rmax <= 1000: - sigma = 25 - elif Rmax > 1000: +def alphadiv(NalphaRange,Rmax): + angle = [-90,90] + if Rmax >100 and Rmax<=1000: + sigma = 20 + elif Rmax >=1000: sigma = 15 - else: - sigma = 40 - for i in range(1, NalphaRange): - angle.append(norm.ppf(i / NalphaRange, loc=0, scale=sigma)) + else : + sigma = 25 + for i in range(1,NalphaRange): + angle.append(norm.ppf(i/NalphaRange,loc=0, scale=sigma)) angle.sort() return angle -def beamsnumber(freq, Rmax, depth): - if Rmax > 700: - nbeams = max(int(0.01 * Rmax * 1000 * freq / 1500), 300) - elif Rmax > 300 and Rmax <= 700: - nbeams = max(int(0.05 * Rmax * 1000 * freq / 1500), 300) - elif Rmax > 100 and Rmax <= 300: - nbeams = max(int(0.1 * Rmax * 1000 * freq / 1500), 300) - else: - nbeams = max(int(0.3 * Rmax * 1000 * freq / 1500), 300) - d_theta_recommended = math.atan(depth / (10.0 * Rmax * 1000)) - nbeams = max(int(math.pi / d_theta_recommended), nbeams) + +def beamsnumber(freq,Rmax,depth): + if Rmax >700 : + nbeams = max(int(0.01 * Rmax*1000 * freq /1500),300) + elif Rmax >300 and Rmax<=700: + nbeams = max(int(0.05 * Rmax*1000 * freq / 1500), 300) + elif Rmax >100 and Rmax<=300: + nbeams = max(int(0.1 * Rmax*1000 * freq / 1500), 300) + else: + nbeams = max(int(0.3 * Rmax*1000 * freq / 1500), 300) + + d_theta_recommended =math.atan(depth/(10.0 * Rmax*1000)) + nbeams = min(int(math.pi/d_theta_recommended),nbeams,3000) return nbeams # if __name__=="__main__": -''' ********************************************** ''' -''' ******************** Test ******************** ''' -# freq = 25 # Frequency -# sd = np.array([10.0]) # Source depth -# rd = np.linspace(0, 4000, 401) # Receiving depth -# # 定义坐标 -# lat = 19.00 -# lon = 116.00 -# # 定义时间段 -# timeIndex = db.timeDict['Annual'] -# # 创建类 -# data = db.DataBase() -# # 设置时间段 -# data.set_timePeriod(timeIndex) + ''' ********************************************** ''' + ''' ******************** Test ******************** ''' + # freq = 25 # Frequency + # sd = np.array([10.0]) # Source depth + # rd = np.linspace(0, 4000, 401) # Receiving depth + # # 定义坐标 + # lat = 19.00 + # lon = 116.00 + # # 定义时间段 + # timeIndex = db.timeDict['Annual'] + # # 创建类 + # data = db.DataBase() + # # 设置时间段 + # data.set_timePeriod(timeIndex) -# [bathm, ssp, sed, base] = GetEnv(data, freq, Coord(lon, lat), R=100, azi=90) -# [Pos, TL] = call_Bellhop(freq, sd, rd, bathm, ssp, sed, base) + # [bathm, ssp, sed, base] = GetEnv(data, freq, Coord(lon, lat), R=100, azi=90) + # [Pos, TL] = call_Bellhop(freq, sd, rd, bathm, ssp, sed, base) -# #[x, x, x, x, Pos, pressure] = read_shd('envB.shd') -# #TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss -# # Plot -# plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130,cmap='jet_r', shading='auto') -# plt.xlabel("Range/m") -# plt.ylabel("Depth/m") -# plt.gca().invert_yaxis() -# plt.colorbar() -# plt.plot(bathm.r*1000, bathm.d, color='k') -# plt.show() \ No newline at end of file + # #[x, x, x, x, Pos, pressure] = read_shd('envB.shd') + # #TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss + # # Plot + # plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130,cmap='jet_r', shading='auto') + # plt.xlabel("Range/m") + # plt.ylabel("Depth/m") + # plt.gca().invert_yaxis() + # plt.colorbar() + # plt.plot(bathm.r*1000, bathm.d, color='k') + # plt.show() \ No newline at end of file diff --git a/pyat/czcore.py b/pyat/czcore.py index 4ad25bd..34aa5a1 100644 --- a/pyat/czcore.py +++ b/pyat/czcore.py @@ -52,7 +52,7 @@ def call_Bellhop(freq, sd, bathm, ssp, sed, base,azi): # Source and receiving position filename = 'cz'+str(azi) Rmax = bathm.r[-1] # unit: km - dr = (1500 / freq) / 1000 # wavelength / 10, unit: km + dr = 0.1 # wavelength / 10, unit: km ran = np.linspace(0, Rmax, round(Rmax / dr)) # Receiving range (unit: km) Zmax, ssp_idx = calZmax(ssp, sed) """RD 已修改 8.26""" @@ -100,7 +100,7 @@ def call_Bellhop(freq, sd, bathm, ssp, sed, base,azi): bdy = Bndry(top, bottom) # Beam params run_type = 'I' # incoherently sum beams, see AT docs for more info - nbeams = 3000 + nbeams = 1000 alpha = np.linspace(-90, 90, 2) # min and max launch angle -20 degrees to 20 degrees box = Box(Zmax, 100) # bound the region you let the beams go, depth in meters and range in km deltas = 0 # length step of ray trace, 0 means automatically choose @@ -122,7 +122,7 @@ def call_Bellhop_Rays(freq, sd, bathm, ssp, sed, base,azi): # Source and receiving position filename = 'cz_ray'+str(azi) Rmax = bathm.r[-1] # unit: km - dr = (1500 / freq) / 1000 # wavelength / 10, unit: km + dr = 0.1 # wavelength / 10, unit: km ran = np.linspace(0, Rmax, round(Rmax / dr)) # Receiving range (unit: km) Zmax, ssp_idx = calZmax(ssp, sed) """RD 已修改 8.26""" @@ -297,7 +297,7 @@ def call_fangwei(freq, sd, bathm, ssp, sed, base, azi, isPlot=False): czdetail.reflectCz = [] return Pos.r.range, Pos.r.depth, TL[0, 0, :, :], bathm, Rays, czdetail -def child(database, freq, lat, lon, sd, rd, i): +def child(database, freq, lat, lon, sd, isPlot, i): descripts = ["正北方向", "东北方向", "正东方向", "东南方向", "正南方向", "西南方向", "正西方向", "西北方向"] theta = np.arange(0, 359, 45) [bathm, ssp, sed, base] = GetEnv(database, 25, Coord(lon, lat), R=100, azi=theta[i]) @@ -306,7 +306,7 @@ def child(database, freq, lat, lon, sd, rd, i): ddd=[] return None else: - ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, rd, bathm, ssp, sed, base, azi=theta[i],isPlot=True) + ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, bathm, ssp, sed, base, azi=theta[i],isPlot=isPlot) czdetail.descript = descripts[i] return CzData(theta[i], ranges, depths, TL, bathym, Rays, czdetail) @@ -325,7 +325,7 @@ def call_cz(database, freq, lat, lon, azi, sourceDepth=10, isPlot=False): if np.mean(bathm.d) < 2000: return None else: - ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, rd, bathm, ssp, sed, base, azi, + ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, bathm, ssp, sed, base, azi, isPlot=isPlot) czdetail.descript = "{:.1f}°方向".format(azi) results.append(CzData(0, ranges, depths, TL, bathym, Rays, czdetail)) @@ -340,7 +340,7 @@ def call_cz(database, freq, lat, lon, azi, sourceDepth=10, isPlot=False): pool = Pool(processes=8) applyResults = [] for pi in range(8): - applyResults.append(pool.apply_async(partial(child, database, freq, lat, lon, sd, rd), (pi,))) + applyResults.append(pool.apply_async(partial(child, database, freq, lat, lon, sd, isPlot), (pi,))) # pi=Process(target=child_0,args=(pi,)) # pi.start() pool.close() @@ -348,3 +348,12 @@ def call_cz(database, freq, lat, lon, azi, sourceDepth=10, isPlot=False): for applyResult in applyResults: results.append(applyResult.get()) return results +# timeIndex = '20190101' +# # 设置时间段 +# +# STPath = '/home/gan/下载/code' +# ETOPO1Path = '/home/gan/下载/code/20220920/ETOPO1' +# database = DataBase(ETOPO1Path, STPath) +# database.set_timePeriod(timeIndex) +# res = call_cz(database, 300, 13, 115, 0, 60, True) +# a = 1 diff --git a/pyat/krakenc.py b/pyat/krakenc.py index d1ee796..a833a61 100755 --- a/pyat/krakenc.py +++ b/pyat/krakenc.py @@ -1,8 +1,8 @@ from os import system import numpy as np -from pyat.readwrite import write_env, write_fieldflp, read_shd -from pyat.env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry -from pyat.config import AtBinPath +from readwrite import write_env, write_fieldflp, read_shd +from env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry +from config import AtBinPath def calZmax(ssp, sed): d = np.array([]) @@ -29,6 +29,7 @@ def call_krakenc(freq, sd, bathm, ssp, sed, base): RD = np.linspace(0, Zmax, round(Zmax) + 1) pos = Pos(Source(sd), Dom(ran, RD)) + # Range of phase velocity cInt.Low = 1400 cInt.High = 15000 @@ -108,6 +109,7 @@ def call_krakenc(freq, sd, bathm, ssp, sed, base): return Pos1, TL + # if __name__=="__main__": ''' ********************************************************************************* ''' ''' ***************************** Input Parameters ********************************** ''' diff --git a/pyat/ram.py b/pyat/ram.py index 1c992fb..459f2c6 100755 --- a/pyat/ram.py +++ b/pyat/ram.py @@ -1,7 +1,8 @@ -from pyat.ramRW import PosRAM, PROF, SSPRAM, write_env_RAM, readTL_RAM +from ramRW import PosRAM, PROF, SSPRAM, write_env_RAM, readTL_RAM from os import system import numpy as np -from pyat.config import AtBinPath +from config import AtBinPath + def calZmax(ssp, sed): d = np.array([]) for it in range(len(ssp)): @@ -11,6 +12,7 @@ def calZmax(ssp, sed): d = np.append(d, ssp[it].z[-1] + sed[it].z[-1]) Zmax = max(d) return Zmax + def call_ram(freq, sd, bathm, ssp, sed): ''' ********************************************************************************* ''' ''' ************************** Sound Field Calculation ****************************** ''' @@ -50,6 +52,7 @@ def call_ram(freq, sd, bathm, ssp, sed): # Read tl.grid return readTL_RAM() + # if __name__=="__main__": # ''' ********************************************************************************* ''' # ''' ***************************** Input Parameters ********************************** ''' diff --git a/pyat/tlcore.py b/pyat/tlcore.py index a03d494..7717eac 100644 --- a/pyat/tlcore.py +++ b/pyat/tlcore.py @@ -1,9 +1,9 @@ import numpy as np -from pyat.AcousticPara import GetEnv, Coord -from pyat.DataBase import DataBase -from pyat.bellhop import call_Bellhop -from pyat.krakenc import call_krakenc -from pyat.ram import call_ram +from AcousticPara import GetEnv, Coord +from DataBase import DataBase +from bellhop import call_Bellhop +from krakenc import call_krakenc +from ram import call_ram import matplotlib.pyplot as plt def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, @@ -79,11 +79,13 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, -# # 创建类 -# database = DataBase() + +# timeIndex = '20190101' # # 设置时间段 -# database.set_timePeriod('06') -# # 设置时间段 -# #database.set_timePeriod(timeIndex) +# +# STPath = '/home/gan/下载/code' +# ETOPO1Path = '/home/gan/下载/code/20220920/ETOPO1' +# database = DataBase(ETOPO1Path, STPath) +# database.set_timePeriod(timeIndex) # range, depth, tl, bathym = call_tl(database=database,freq=200, sourceDepth=60, sourceLat=19, sourceLon=127, -# receiveLat=17.5, receiveLon=127.1, modelType='ram', isPlot=True) +# receiveLat=17.5, receiveLon=129, modelType='bellhop', isPlot=True) diff --git a/tl.py b/tl.py index 29630d2..65e6216 100644 --- a/tl.py +++ b/tl.py @@ -1,31 +1,36 @@ # -*- coding: utf-8 -*- import warnings + warnings.filterwarnings('ignore') import argparse, json from common import NoIndent, MyEncoder import os, sys -sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))+'/pyat') +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + '/pyat') + + +# import time def _argparse(): parser = argparse.ArgumentParser(description='tl client') parser.add_argument('-a', '--argfile', action='store', dest='argfile', required=True, help='name of arguments file') return parser.parse_args() + def main(): parser = _argparse() argfile = parser.argfile if argfile[-5:] != '.json': print('输入的参数文件名有误') return - with open(argfile,'r') as jsonFile: + with open(argfile, 'r') as jsonFile: paramData = json.load(jsonFile) - if (("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or\ - "freq" not in paramData) or ("sourceDepth" not in paramData) or\ - ("receiveDepth" not in paramData) or ("sLat" not in paramData) or\ - ("sLon" not in paramData) or ("rLat" not in paramData) or\ - ("rLon" not in paramData) or ("timeType" not in paramData) or\ - ("modelType" not in paramData) or ("png" not in paramData) or\ - ("outParam" not in paramData): + if (("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or \ + "freq" not in paramData) or ("sourceDepth" not in paramData) or \ + ("receiveDepth" not in paramData) or ("sLat" not in paramData) or \ + ("sLon" not in paramData) or ("rLat" not in paramData) or \ + ("rLon" not in paramData) or ("timeType" not in paramData) or \ + ("modelType" not in paramData) or ("png" not in paramData) or \ + ("outParam" not in paramData): print(f'{argfile}中的参数有误，请核实！') return ETOPO1Path = paramData["ETOPO1Path"] @@ -41,11 +46,11 @@ def main(): modelType = paramData["modelType"] png = paramData["png"] outParam = paramData["outParam"] - + from DataBase import DataBase from tlcore import call_tl import numpy as np - + # 定义时间段 timeIndex = timeType # 创建类 @@ -53,10 +58,11 @@ def main(): # 设置时间段 database.set_timePeriod(timeIndex) - [ranges, depths, TL, bathym] = call_tl(database=database, freq=freq, sourceDepth=sourceDepth, sourceLat=sLat, sourceLon=sLon, receiveDepth=receiveDepth, - receiveLat=rLat, receiveLon=rLon, modelType=modelType, isPlot=png) + [ranges, depths, TL, bathym] = call_tl(database=database, freq=freq, sourceDepth=sourceDepth, sourceLat=sLat, + sourceLon=sLon, + receiveLat=rLat, receiveLon=rLon, modelType=modelType, isPlot=png) - # write data to file + # write data to file # 保存地形数据 bathym.r = np.around(bathym.r, decimals=0) bathym.d = np.around(bathym.d, decimals=0) @@ -66,7 +72,8 @@ def main(): "descript": "距离、深度", "data": NoIndent(list(zip(bathym.r, bathym.d))) } - with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"tlOut_bathym.json", 'w') as jsonFile: + with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "tlOut_bathym.json", + 'w') as jsonFile: json_data = json.dumps(bathymData, cls=MyEncoder, indent=2, ensure_ascii=False) jsonFile.write(json_data) @@ -75,16 +82,21 @@ def main(): depths = np.around(depths, decimals=0) ranges = np.around(ranges, decimals=0) TL = np.around(TL, decimals=2) - tlpoint = list(zip(np.tile(depths.reshape(len(depths),1),len(ranges)).flatten(), np.tile(ranges,(len(depths),1)).flatten(), TL.flatten())) + tlpoint = list(zip(np.tile(depths.reshape(len(depths), 1), len(ranges)).flatten(), + np.tile(ranges, (len(depths), 1)).flatten(), TL.flatten())) tlData = { "param": paramData, - "count": len(ranges)*len(depths), + "count": len(ranges) * len(depths), "descript": "距离、深度、损失", "data": NoIndent(tlpoint) } - with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"tlOut_tl.json", 'w') as jsonFile: + with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "tlOut_tl.json", 'w') as jsonFile: json_data = json.dumps(tlData, cls=MyEncoder, indent=2, ensure_ascii=False) jsonFile.write(json_data) - + + if __name__ == '__main__': - main() \ No newline at end of file + # start = time.perf_counter() + main() + # end = time.perf_counter() + # print('time = ',(end-start)) \ No newline at end of file