From 2564b0166073790d668bff7ecb996452756e0bd1 Mon Sep 17 00:00:00 2001 From: Hao Date: Mon, 8 Jul 2024 07:37:55 +0000 Subject: [PATCH 01/10] add exporter of velest module --- pyrocko | 1 + src/qseek/exporters/velest.py | 269 +++++++++++++++++++++++++++++++++- 2 files changed, 264 insertions(+), 6 deletions(-) create mode 160000 pyrocko diff --git a/pyrocko b/pyrocko new file mode 160000 index 00000000..7148797f --- /dev/null +++ b/pyrocko @@ -0,0 +1 @@ +Subproject commit 7148797f3bf57e37d944c2274c3a7a71e83cbc41 diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index 7286650a..a16d303c 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -2,34 +2,291 @@ import logging from pathlib import Path - +from io import TextIOWrapper import rich from rich.prompt import FloatPrompt from qseek.exporters.base import Exporter from qseek.search import Search -logger = logging.getLogger(__name__) +from qseek.models.station import Station, Stations +from qseek.models.detection import EventDetection +from qseek.utils import PhaseDescription +logger = logging.getLogger(__name__) class Velest(Exporter): """Crate a VELEST project folder for 1D velocity model estimation.""" - min_pick_semblance: float = 0.3 + min_pick_semblance: float = 0.2 + pp:float = 0.3 + ps:float = 0.3 + tdiff:float = 2.5 n_picks: dict[str, int] = {} n_events: int = 0 async def export(self, rundir: Path, outdir: Path) -> Path: rich.print("Exporting qseek search to VELEST project folder") - min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.3) - + min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.2) + pp = FloatPrompt.ask("Minimum pick probability for P phase", default=0.3) + ps = FloatPrompt.ask("Minimum pick probability for S phase", default=0.3) + tdiff = FloatPrompt.ask("Maximum difference between model and observed arrival", default=2.5) self.min_pick_semblance = min_pick_semblance + self.pp = pp + self.ps = ps + self.tdiff = tdiff + self.n_picks= {"P":0,"S":0} outdir.mkdir() search = Search.load_rundir(rundir) + phases = search.image_functions.get_phases() + for phase in phases: + if 'P' in phase: + phaseP=phase + if 'S' in phase: + phaseS=phase + if 'phaseP' not in locals() and 'phaseS' not in locals(): + print("Flail to find both P and S phase name." ) + return outdir + catalog = search.catalog # noqa + #export station file + stations=search.stations.stations + station_file=outdir/"stations_velest.sta" + self.export_station(stations=stations,filename=station_file) + + #export phase file + phase_file = outdir / "phase_velest.pha" + fp=open(phase_file,'w') + neq=0 + for event in catalog: + if event.semblance > min_pick_semblance: + countp, counts=self.export_phase(event,fp,phaseP,phaseS,tdiff,pp,ps) + self.n_picks['P']+=countp + self.n_picks['S']+=counts + neq+=1 + self.n_events=neq + #export control file + params = { + "reflat": search.octree.location.lat, # Reference Latitude + "reflon": -search.octree.location.lon, # Reference Longitude (should be negative!) + "neqs": neq, # Number of earthquakes + "distmax": 200, # Maximum distance from the station + "zmin": -0.2, # Minimum depth + "lowvelocity": 0, # Allow low velocity zones (0 or 1) + "vthet": 1, # Damping parameter for the velocity + "stathet": 0.1, # Damping parameter for the station + "iuseelev": 0, # Use elevation (0 or 1) + "iusestacorr": 0 # Use station corrections (0 or 1) + } + control_file=outdir / "velest.cmn" + self.make_control_file(control_file,params, 1,'model.mod','phase_velest.pha','main.out','log.out','final.cnv','stacor.dat') + + #export velocity model file + dep=search.ray_tracers.root[0].earthmodel.layered_model.profile("z") + vp=search.ray_tracers.root[0].earthmodel.layered_model.profile("vp") + vs=search.ray_tracers.root[0].earthmodel.layered_model.profile("vs") + dep_velest=[] + vp_velest=[] + vs_velest=[] + for i, d in enumerate(dep): + if float(d)/1000 not in dep_velest: + dep_velest.append(float(d)/1000) + vp_velest.append(float(vp[i])/1000) + vs_velest.append(float(vs[i])/1000) + velmod_file=outdir / "model.mod" + self.make_velmod_file(velmod_file,vp_velest,vs_velest,dep_velest ) + export_info = outdir / "export_info.json" export_info.write_text(self.model_dump_json(indent=2)) - + print("Done!") return outdir + + @staticmethod + def export_phase( + event: EventDetection, + fp:TextIOWrapper, + phaseP:PhaseDescription, + phaseS:PhaseDescription, + tdiff:float, + pp:float, + ps:float, + ) -> int: + year=int(str(event.time.year)[-2::]) + # export phase catalog into txt file (for velest, format: ised=1) + # tdiff: time diff tobs-tmodel threshold P,S + # pp,ps : minimun probability threshold P,S + # there is quality weigth (0-4, 0 is best) in velest, here use probability to define it + month=event.time.month + day=event.time.day + hour=event.time.hour + min=event.time.minute + sec=event.time.second+event.time.microsecond/1000000 + lat=event.effective_lat + lon=event.effective_lon + dep=event.depth/1000 + if event.magnitude !=None: + mag=event.magnitude.average + else: + mag=0.0 + if lat<0: + vsn='S' + lat=abs(lat) + else: + vsn='N' + if lon<0: + vew='W' + lon=abs(lon) + else: + vew='E' + fp.write("%2d%2d%2d %2d%2d %5.2f %7.4f%s %8.4f%s %7.2f %5.2f\n"%(year,month,day,hour,min,sec,lat,vsn,lon,vew,dep,mag)) + countP=0 + countS=0 + for receiver in event.receivers.receivers: + station=receiver.station + if receiver.phase_arrivals[phaseP].observed !=None and receiver.phase_arrivals[phaseP].observed.detection_value >=pp: + tpick=receiver.phase_arrivals[phaseP].observed.time-event.time + tpick=tpick.total_seconds() + dt=receiver.phase_arrivals[phaseP].traveltime_delay.total_seconds() + if abs(dt)>tdiff: + continue + if receiver.phase_arrivals[phaseP].observed.detection_value <0.4: + iwt=3 + elif receiver.phase_arrivals[phaseP].observed.detection_value <0.6: + iwt=2 + elif receiver.phase_arrivals[phaseP].observed.detection_value <0.8: + iwt=1 + else: + iwt=0 + phase='P' + fp.write(" %-6s %-1s %1d %6.2f\n"%(station,phase,iwt,tpick)) + countP+=1 + if receiver.phase_arrivals[phaseS].observed !=None and receiver.phase_arrivals[phaseS].observed.detection_value >=ps: + tpick=receiver.phase_arrivals[phaseS].observed.time-event.time + tpick=tpick.total_seconds() + dt=receiver.phase_arrivals[phaseS].traveltime_delay.total_seconds() + if abs(dt)>tdiff: + continue + if receiver.phase_arrivals[phaseS].observed.detection_value <0.4: + iwt=3 + elif receiver.phase_arrivals[phaseS].observed.detection_value <0.6: + iwt=2 + elif receiver.phase_arrivals[phaseS].observed.detection_value <0.8: + iwt=1 + else: + iwt=0 + phase='S' + fp.write(" %-6s %-1s %1d %6.2f\n"%(station,phase,iwt,tpick)) + countS+=1 + if countP == 0 and countS == 0 : + print("No good phases obesered for event%s, lower pick probability or increase pick confidence"%event.time) + fp.write("\n") + return countP , countS + + @staticmethod + def export_station( + stations: list[Station], + filename:Path + )-> None: + fpout=open(filename,'w') + fpout.write("(a6,f7.4,a1,1x,f8.4,a1,1x,i4,1x,i1,1x,i3,1x,f5.2,2x,f5.2)\n") + p2=1 + p3=1 + v1=0.00 + v2=0.00 + for station in stations: + lat=station.lat + lon=station.lon + sta=station.station + elev=station.elevation + if lat<0: + vsn='S' + lat=abs(lat) + else: + vsn='N' + if lon<0: + vew='W' + lon=abs(lon) + else: + vew='E' + fpout.write("%-6s%7.4f%1s %8.4f%s %4d %1d %3d %5.2f %5.2f\n"%(sta,lat,vsn,lon,vew,elev,p2,p3,v1,v2)) + p3+=1 + fpout.close() + + @staticmethod + def make_control_file(filename, params,isingle,modname,phasefile,mainoutfile,outcheckfile,finalcnv,stacorfile): + + if isingle == 1: + ittmax = 99 + invertratio = 0 + else: + ittmax = 9 + invertratio = 3 + fp=open(filename,'w') + fp.write("velest parameters are below, please modify according to their documents\n") + #*** olat olon icoordsystem zshift itrial ztrial ised + fp.write("%s %s 0 0.0 0 0.00 1\n"%(params["reflat"],params["reflon"])) + #*** neqs nshot rotate + fp.write("%s 0 0.0\n"%params["neqs"]) + #*** isingle iresolcalc + fp.write("%s 0\n"%isingle) + #*** dmax itopo zmin veladj zadj lowveloclay + fp.write("%s 0 %s 0.20 5.00 %s\n"%(params["distmax"],params["zmin"],params["lowvelocity"] )) + #*** nsp swtfac vpvs nmod + fp.write( "2 0.75 1.650 1\n") + #*** othet xythet zthet vthet stathet + fp.write( "0.01 0.01 0.01 %s %s\n"%(params["vthet"],params["stathet"])) + #*** nsinv nshcor nshfix iuseelev iusestacorr + fp.write( "1 0 0 %s %s\n"%(params["iuseelev"],params["iusestacorr"])) + #*** iturbo icnvout istaout ismpout + fp.write( "1 1 2 0\n") + #*** irayout idrvout ialeout idspout irflout irfrout iresout + fp.write( "0 0 0 0 0 0 0\n") + #*** delmin ittmax invertratio + fp.write( "0.001 %s %s\n"%(ittmax,invertratio)) + #*** Modelfile: + fp.write( "%s\n"%modname) + #*** Stationfile: + fp.write( "stations_velest.sta\n") + #*** Seismofile: + fp.write( " \n") + #*** File with region names: + fp.write( "regionsnamen.dat\n") + #*** File with region coordinates: + fp.write( "regionskoord.dat\n") + #*** File #1 with topo data: + fp.write( " \n") + #*** File #2 with topo data: + fp.write( " \n") + #*** File with Earthquake data (phase): + fp.write( "%s\n"%phasefile) + #*** File with Shot data: + fp.write( " \n") + #*** Main print output file: + fp.write( "%s\n"%mainoutfile) + #*** File with single event locations: + fp.write( "%s\n"%outcheckfile) + #*** File with final hypocenters in *.cnv format: + fp.write( "%s\n"%finalcnv) + #*** File with new station corrections: + fp.write( "%s\n"%stacorfile) + fp.close() + + @staticmethod + def make_velmod_file(modname,vp,vs,dep ): + nlayer=len(dep) + vdamp = 1.0 + fp=open(modname,'w') + fp.write("initial 1D-model for velest\n") + # the second line - indicate the number of layers for Vp + fp.write("%3d vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n"%nlayer) + #vp model + for i,v in enumerate(vp): + fp.write("%5.2f %7.2f %7.3f\n"%(v,dep[i],vdamp)) + #vs model + fp.write("%3d\n"%nlayer) + for i,v in enumerate(vs): + fp.write("%5.2f %7.2f %7.3f\n"%(v,dep[i],vdamp)) + fp.close() + \ No newline at end of file From eaceff8fcf99915e2dab28e94c4bbbe98e0afe7e Mon Sep 17 00:00:00 2001 From: Hao Date: Mon, 8 Jul 2024 07:45:38 +0000 Subject: [PATCH 02/10] add exporter of velest module --- src/qseek/exporters/velest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index a16d303c..77c26e2e 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -45,7 +45,7 @@ async def export(self, rundir: Path, outdir: Path) -> Path: phaseP=phase if 'S' in phase: phaseS=phase - if 'phaseP' not in locals() and 'phaseS' not in locals(): + if 'phaseP' not in locals() or 'phaseS' not in locals(): print("Flail to find both P and S phase name." ) return outdir From d6965cf06fd0bb580095a632082e3a91e29f2706 Mon Sep 17 00:00:00 2001 From: Hao Date: Mon, 8 Jul 2024 07:37:55 +0000 Subject: [PATCH 03/10] add exporter of velest module --- pyrocko | 1 + src/qseek/exporters/velest.py | 269 +++++++++++++++++++++++++++++++++- 2 files changed, 264 insertions(+), 6 deletions(-) create mode 160000 pyrocko diff --git a/pyrocko b/pyrocko new file mode 160000 index 00000000..7148797f --- /dev/null +++ b/pyrocko @@ -0,0 +1 @@ +Subproject commit 7148797f3bf57e37d944c2274c3a7a71e83cbc41 diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index 7286650a..a16d303c 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -2,34 +2,291 @@ import logging from pathlib import Path - +from io import TextIOWrapper import rich from rich.prompt import FloatPrompt from qseek.exporters.base import Exporter from qseek.search import Search -logger = logging.getLogger(__name__) +from qseek.models.station import Station, Stations +from qseek.models.detection import EventDetection +from qseek.utils import PhaseDescription +logger = logging.getLogger(__name__) class Velest(Exporter): """Crate a VELEST project folder for 1D velocity model estimation.""" - min_pick_semblance: float = 0.3 + min_pick_semblance: float = 0.2 + pp:float = 0.3 + ps:float = 0.3 + tdiff:float = 2.5 n_picks: dict[str, int] = {} n_events: int = 0 async def export(self, rundir: Path, outdir: Path) -> Path: rich.print("Exporting qseek search to VELEST project folder") - min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.3) - + min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.2) + pp = FloatPrompt.ask("Minimum pick probability for P phase", default=0.3) + ps = FloatPrompt.ask("Minimum pick probability for S phase", default=0.3) + tdiff = FloatPrompt.ask("Maximum difference between model and observed arrival", default=2.5) self.min_pick_semblance = min_pick_semblance + self.pp = pp + self.ps = ps + self.tdiff = tdiff + self.n_picks= {"P":0,"S":0} outdir.mkdir() search = Search.load_rundir(rundir) + phases = search.image_functions.get_phases() + for phase in phases: + if 'P' in phase: + phaseP=phase + if 'S' in phase: + phaseS=phase + if 'phaseP' not in locals() and 'phaseS' not in locals(): + print("Flail to find both P and S phase name." ) + return outdir + catalog = search.catalog # noqa + #export station file + stations=search.stations.stations + station_file=outdir/"stations_velest.sta" + self.export_station(stations=stations,filename=station_file) + + #export phase file + phase_file = outdir / "phase_velest.pha" + fp=open(phase_file,'w') + neq=0 + for event in catalog: + if event.semblance > min_pick_semblance: + countp, counts=self.export_phase(event,fp,phaseP,phaseS,tdiff,pp,ps) + self.n_picks['P']+=countp + self.n_picks['S']+=counts + neq+=1 + self.n_events=neq + #export control file + params = { + "reflat": search.octree.location.lat, # Reference Latitude + "reflon": -search.octree.location.lon, # Reference Longitude (should be negative!) + "neqs": neq, # Number of earthquakes + "distmax": 200, # Maximum distance from the station + "zmin": -0.2, # Minimum depth + "lowvelocity": 0, # Allow low velocity zones (0 or 1) + "vthet": 1, # Damping parameter for the velocity + "stathet": 0.1, # Damping parameter for the station + "iuseelev": 0, # Use elevation (0 or 1) + "iusestacorr": 0 # Use station corrections (0 or 1) + } + control_file=outdir / "velest.cmn" + self.make_control_file(control_file,params, 1,'model.mod','phase_velest.pha','main.out','log.out','final.cnv','stacor.dat') + + #export velocity model file + dep=search.ray_tracers.root[0].earthmodel.layered_model.profile("z") + vp=search.ray_tracers.root[0].earthmodel.layered_model.profile("vp") + vs=search.ray_tracers.root[0].earthmodel.layered_model.profile("vs") + dep_velest=[] + vp_velest=[] + vs_velest=[] + for i, d in enumerate(dep): + if float(d)/1000 not in dep_velest: + dep_velest.append(float(d)/1000) + vp_velest.append(float(vp[i])/1000) + vs_velest.append(float(vs[i])/1000) + velmod_file=outdir / "model.mod" + self.make_velmod_file(velmod_file,vp_velest,vs_velest,dep_velest ) + export_info = outdir / "export_info.json" export_info.write_text(self.model_dump_json(indent=2)) - + print("Done!") return outdir + + @staticmethod + def export_phase( + event: EventDetection, + fp:TextIOWrapper, + phaseP:PhaseDescription, + phaseS:PhaseDescription, + tdiff:float, + pp:float, + ps:float, + ) -> int: + year=int(str(event.time.year)[-2::]) + # export phase catalog into txt file (for velest, format: ised=1) + # tdiff: time diff tobs-tmodel threshold P,S + # pp,ps : minimun probability threshold P,S + # there is quality weigth (0-4, 0 is best) in velest, here use probability to define it + month=event.time.month + day=event.time.day + hour=event.time.hour + min=event.time.minute + sec=event.time.second+event.time.microsecond/1000000 + lat=event.effective_lat + lon=event.effective_lon + dep=event.depth/1000 + if event.magnitude !=None: + mag=event.magnitude.average + else: + mag=0.0 + if lat<0: + vsn='S' + lat=abs(lat) + else: + vsn='N' + if lon<0: + vew='W' + lon=abs(lon) + else: + vew='E' + fp.write("%2d%2d%2d %2d%2d %5.2f %7.4f%s %8.4f%s %7.2f %5.2f\n"%(year,month,day,hour,min,sec,lat,vsn,lon,vew,dep,mag)) + countP=0 + countS=0 + for receiver in event.receivers.receivers: + station=receiver.station + if receiver.phase_arrivals[phaseP].observed !=None and receiver.phase_arrivals[phaseP].observed.detection_value >=pp: + tpick=receiver.phase_arrivals[phaseP].observed.time-event.time + tpick=tpick.total_seconds() + dt=receiver.phase_arrivals[phaseP].traveltime_delay.total_seconds() + if abs(dt)>tdiff: + continue + if receiver.phase_arrivals[phaseP].observed.detection_value <0.4: + iwt=3 + elif receiver.phase_arrivals[phaseP].observed.detection_value <0.6: + iwt=2 + elif receiver.phase_arrivals[phaseP].observed.detection_value <0.8: + iwt=1 + else: + iwt=0 + phase='P' + fp.write(" %-6s %-1s %1d %6.2f\n"%(station,phase,iwt,tpick)) + countP+=1 + if receiver.phase_arrivals[phaseS].observed !=None and receiver.phase_arrivals[phaseS].observed.detection_value >=ps: + tpick=receiver.phase_arrivals[phaseS].observed.time-event.time + tpick=tpick.total_seconds() + dt=receiver.phase_arrivals[phaseS].traveltime_delay.total_seconds() + if abs(dt)>tdiff: + continue + if receiver.phase_arrivals[phaseS].observed.detection_value <0.4: + iwt=3 + elif receiver.phase_arrivals[phaseS].observed.detection_value <0.6: + iwt=2 + elif receiver.phase_arrivals[phaseS].observed.detection_value <0.8: + iwt=1 + else: + iwt=0 + phase='S' + fp.write(" %-6s %-1s %1d %6.2f\n"%(station,phase,iwt,tpick)) + countS+=1 + if countP == 0 and countS == 0 : + print("No good phases obesered for event%s, lower pick probability or increase pick confidence"%event.time) + fp.write("\n") + return countP , countS + + @staticmethod + def export_station( + stations: list[Station], + filename:Path + )-> None: + fpout=open(filename,'w') + fpout.write("(a6,f7.4,a1,1x,f8.4,a1,1x,i4,1x,i1,1x,i3,1x,f5.2,2x,f5.2)\n") + p2=1 + p3=1 + v1=0.00 + v2=0.00 + for station in stations: + lat=station.lat + lon=station.lon + sta=station.station + elev=station.elevation + if lat<0: + vsn='S' + lat=abs(lat) + else: + vsn='N' + if lon<0: + vew='W' + lon=abs(lon) + else: + vew='E' + fpout.write("%-6s%7.4f%1s %8.4f%s %4d %1d %3d %5.2f %5.2f\n"%(sta,lat,vsn,lon,vew,elev,p2,p3,v1,v2)) + p3+=1 + fpout.close() + + @staticmethod + def make_control_file(filename, params,isingle,modname,phasefile,mainoutfile,outcheckfile,finalcnv,stacorfile): + + if isingle == 1: + ittmax = 99 + invertratio = 0 + else: + ittmax = 9 + invertratio = 3 + fp=open(filename,'w') + fp.write("velest parameters are below, please modify according to their documents\n") + #*** olat olon icoordsystem zshift itrial ztrial ised + fp.write("%s %s 0 0.0 0 0.00 1\n"%(params["reflat"],params["reflon"])) + #*** neqs nshot rotate + fp.write("%s 0 0.0\n"%params["neqs"]) + #*** isingle iresolcalc + fp.write("%s 0\n"%isingle) + #*** dmax itopo zmin veladj zadj lowveloclay + fp.write("%s 0 %s 0.20 5.00 %s\n"%(params["distmax"],params["zmin"],params["lowvelocity"] )) + #*** nsp swtfac vpvs nmod + fp.write( "2 0.75 1.650 1\n") + #*** othet xythet zthet vthet stathet + fp.write( "0.01 0.01 0.01 %s %s\n"%(params["vthet"],params["stathet"])) + #*** nsinv nshcor nshfix iuseelev iusestacorr + fp.write( "1 0 0 %s %s\n"%(params["iuseelev"],params["iusestacorr"])) + #*** iturbo icnvout istaout ismpout + fp.write( "1 1 2 0\n") + #*** irayout idrvout ialeout idspout irflout irfrout iresout + fp.write( "0 0 0 0 0 0 0\n") + #*** delmin ittmax invertratio + fp.write( "0.001 %s %s\n"%(ittmax,invertratio)) + #*** Modelfile: + fp.write( "%s\n"%modname) + #*** Stationfile: + fp.write( "stations_velest.sta\n") + #*** Seismofile: + fp.write( " \n") + #*** File with region names: + fp.write( "regionsnamen.dat\n") + #*** File with region coordinates: + fp.write( "regionskoord.dat\n") + #*** File #1 with topo data: + fp.write( " \n") + #*** File #2 with topo data: + fp.write( " \n") + #*** File with Earthquake data (phase): + fp.write( "%s\n"%phasefile) + #*** File with Shot data: + fp.write( " \n") + #*** Main print output file: + fp.write( "%s\n"%mainoutfile) + #*** File with single event locations: + fp.write( "%s\n"%outcheckfile) + #*** File with final hypocenters in *.cnv format: + fp.write( "%s\n"%finalcnv) + #*** File with new station corrections: + fp.write( "%s\n"%stacorfile) + fp.close() + + @staticmethod + def make_velmod_file(modname,vp,vs,dep ): + nlayer=len(dep) + vdamp = 1.0 + fp=open(modname,'w') + fp.write("initial 1D-model for velest\n") + # the second line - indicate the number of layers for Vp + fp.write("%3d vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n"%nlayer) + #vp model + for i,v in enumerate(vp): + fp.write("%5.2f %7.2f %7.3f\n"%(v,dep[i],vdamp)) + #vs model + fp.write("%3d\n"%nlayer) + for i,v in enumerate(vs): + fp.write("%5.2f %7.2f %7.3f\n"%(v,dep[i],vdamp)) + fp.close() + \ No newline at end of file From 3d63cddcb224b5831299ec889a0e1ee5cc1b0fbc Mon Sep 17 00:00:00 2001 From: Hao Date: Mon, 8 Jul 2024 07:45:38 +0000 Subject: [PATCH 04/10] add exporter of velest module --- src/qseek/exporters/velest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index a16d303c..77c26e2e 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -45,7 +45,7 @@ async def export(self, rundir: Path, outdir: Path) -> Path: phaseP=phase if 'S' in phase: phaseS=phase - if 'phaseP' not in locals() and 'phaseS' not in locals(): + if 'phaseP' not in locals() or 'phaseS' not in locals(): print("Flail to find both P and S phase name." ) return outdir From c4a2a64e0b215cb5e59238804568b9f1e4b82f94 Mon Sep 17 00:00:00 2001 From: Hao Date: Thu, 11 Jul 2024 09:13:31 +0000 Subject: [PATCH 05/10] velest --- src/qseek/exporters/velest.py | 499 +++++++++++++++++----------------- 1 file changed, 254 insertions(+), 245 deletions(-) diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index 77c26e2e..da3d23aa 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -2,291 +2,300 @@ import logging from pathlib import Path -from io import TextIOWrapper +from typing import NamedTuple + import rich from rich.prompt import FloatPrompt from qseek.exporters.base import Exporter +from qseek.models.detection import EventDetection, PhaseDetection, Receiver +from qseek.models.station import Station from qseek.search import Search -from qseek.models.station import Station, Stations -from qseek.models.detection import EventDetection -from qseek.utils import PhaseDescription - logger = logging.getLogger(__name__) +CONTROL_FILE_TPL = """velest parameters are below, please modify according to their documents +{ref_lat} {ref_lon} 0 0.0 0 0.00 1 +{n_earthquakes} 0 0.0 +{isingle:1d} 0 +{max_distance_station} 0 {min_depth} 0.20 5.00 {use_station_correction:1d} +2 0.75 {vp_vs_ratio} 1 +0.01 0.01 0.01 {velocity_damping} {station_correction_damping} +1 0 0 {use_elevation:1d} {use_station_correction:1d} +1 1 2 0 +0 0 0 0 0 0 0 +0.001 {iteration_number} {invertratio} +{model_file} +stations_velest.sta + +regionsnamen.dat +regionskoord.dat + + +{phase_file} + +{mainout_file} +{outcheck_file} +{finalcnv_file} +{stacorrection_file} +""" + + +class VelestControlFile(NamedTuple): + ref_lat: float + ref_lon: float # should be negative for East + n_earthquakes: int + isingle: bool + max_distance_station: float + min_depth: float + allow_low_velocity: bool + vp_vs_ratio: float + velocity_damping: float # Damping parameter for the velocity + station_correction_damping: float # Damping parameter for the station + use_elevation: bool + use_station_correction: bool + iteration_number: int + invertratio: int + model_file: str + phase_file: str + mainout_file: str + outcheck_file: str + finalcnv_file: str + stacorrection_file: str + + def write_config_file(self, file: Path): + with file.open("w") as fp: + fp.write(CONTROL_FILE_TPL.format(**self._asdict())) + + class Velest(Exporter): """Crate a VELEST project folder for 1D velocity model estimation.""" min_pick_semblance: float = 0.2 - pp:float = 0.3 - ps:float = 0.3 - tdiff:float = 2.5 - n_picks: dict[str, int] = {} + min_receivers_number: int = 10 + min_p_phase_confidence: float = 0.3 + min_s_phase_confidence: float = 0.3 + max_traveltime_delay: float = 2.5 + n_picks_p: int = 0 + n_picks_s: int = 0 n_events: int = 0 async def export(self, rundir: Path, outdir: Path) -> Path: rich.print("Exporting qseek search to VELEST project folder") min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.2) - pp = FloatPrompt.ask("Minimum pick probability for P phase", default=0.3) - ps = FloatPrompt.ask("Minimum pick probability for S phase", default=0.3) - tdiff = FloatPrompt.ask("Maximum difference between model and observed arrival", default=2.5) + min_receivers_number = FloatPrompt.ask( + "Minimum number of receivers (P phase)", default=10 + ) + min_p_phase_confidence = FloatPrompt.ask( + "Minimum pick probability for P phase", default=0.3 + ) + min_s_phase_confidence = FloatPrompt.ask( + "Minimum pick probability for S phase", default=0.3 + ) + max_traveltime_delay = FloatPrompt.ask( + "Maximum difference between theoretical and observed arrival", default=2.5 + ) self.min_pick_semblance = min_pick_semblance - self.pp = pp - self.ps = ps - self.tdiff = tdiff - self.n_picks= {"P":0,"S":0} + self.min_receivers_number = min_receivers_number + self.min_p_phase_confidence = min_p_phase_confidence + self.min_s_phase_confidence = min_s_phase_confidence + self.max_traveltime_delay = max_traveltime_delay outdir.mkdir() search = Search.load_rundir(rundir) phases = search.image_functions.get_phases() for phase in phases: - if 'P' in phase: - phaseP=phase - if 'S' in phase: - phaseS=phase - if 'phaseP' not in locals() or 'phaseS' not in locals(): - print("Flail to find both P and S phase name." ) - return outdir - - catalog = search.catalog # noqa + if "P" in phase: + phase_p = phase + if "S" in phase: + phase_s = phase - #export station file - stations=search.stations.stations - station_file=outdir/"stations_velest.sta" - self.export_station(stations=stations,filename=station_file) + catalog = search.catalog - #export phase file + # export station file + stations = search.stations.stations + station_file = outdir / "stations_velest.sta" + self.export_station(stations=stations, filename=station_file) + + # export phase file phase_file = outdir / "phase_velest.pha" - fp=open(phase_file,'w') - neq=0 + n_earthquakes = 0 for event in catalog: - if event.semblance > min_pick_semblance: - countp, counts=self.export_phase(event,fp,phaseP,phaseS,tdiff,pp,ps) - self.n_picks['P']+=countp - self.n_picks['S']+=counts - neq+=1 - self.n_events=neq - #export control file - params = { - "reflat": search.octree.location.lat, # Reference Latitude - "reflon": -search.octree.location.lon, # Reference Longitude (should be negative!) - "neqs": neq, # Number of earthquakes - "distmax": 200, # Maximum distance from the station - "zmin": -0.2, # Minimum depth - "lowvelocity": 0, # Allow low velocity zones (0 or 1) - "vthet": 1, # Damping parameter for the velocity - "stathet": 0.1, # Damping parameter for the station - "iuseelev": 0, # Use elevation (0 or 1) - "iusestacorr": 0 # Use station corrections (0 or 1) - } - control_file=outdir / "velest.cmn" - self.make_control_file(control_file,params, 1,'model.mod','phase_velest.pha','main.out','log.out','final.cnv','stacor.dat') + if event.semblance < min_pick_semblance: + continue + if event.receivers.n_observations(phase_p) < min_receivers_number: + continue + + observed_arrivals: list[tuple[Receiver, PhaseDetection]] = [] + + for receiver in event.receivers: + for _phase, detection in receiver.phase_arrivals.items(): + if detection.observed is None: + continue + observed = detection.observed + if ( + detection.phase == phase_p + and observed.detection_value <= min_p_phase_confidence + ): + continue + if ( + detection.phase == phase_s + and observed.detection_value <= min_s_phase_confidence + ): + continue + if ( + detection.traveltime_delay.total_seconds() + > max_traveltime_delay + ): + continue + observed_arrivals.append((receiver, detection)) - #export velocity model file - dep=search.ray_tracers.root[0].earthmodel.layered_model.profile("z") - vp=search.ray_tracers.root[0].earthmodel.layered_model.profile("vp") - vs=search.ray_tracers.root[0].earthmodel.layered_model.profile("vs") - dep_velest=[] - vp_velest=[] - vs_velest=[] + countp, counts = self.export_phases_slim( + phase_file, event, observed_arrivals + ) + self.n_picks_p += countp + self.n_picks_s += counts + n_earthquakes += 1 + self.n_events = n_earthquakes + + # export control file + control_file = outdir / "velest.cmn" + control_file_parameters = VelestControlFile( + ref_lat=search.octree.location.lat, + ref_lon=-search.octree.location.lon, + n_earthquakes=n_earthquakes, + max_distance_station=200, + min_depth=-0.2, + allow_low_velocity=False, + velocity_damping=1.0, + station_correction_damping=0.1, + use_elevation=False, + use_station_correction=False, + model_file="model.mod", + phase_file="phase_velest.pha", + mainout_file="main.out", + outcheck_file="log.out", + finalcnv_file="final.cnv", + stacorrection_file="stacor.dat", + isingle=False, + vp_vs_ratio=1.65, + iteration_number=99, + invertratio=0, + ) + control_file_parameters.write_config_file(control_file) + # export velocity model file + dep = search.ray_tracers.root[0].earthmodel.layered_model.profile("z") + vp = search.ray_tracers.root[0].earthmodel.layered_model.profile("vp") + vs = search.ray_tracers.root[0].earthmodel.layered_model.profile("vs") + dep_velest = [] + vp_velest = [] + vs_velest = [] for i, d in enumerate(dep): - if float(d)/1000 not in dep_velest: - dep_velest.append(float(d)/1000) - vp_velest.append(float(vp[i])/1000) - vs_velest.append(float(vs[i])/1000) - velmod_file=outdir / "model.mod" - self.make_velmod_file(velmod_file,vp_velest,vs_velest,dep_velest ) + if float(d) / 1000 not in dep_velest: + dep_velest.append(float(d) / 1000) + vp_velest.append(float(vp[i]) / 1000) + vs_velest.append(float(vs[i]) / 1000) + velmod_file = outdir / "model.mod" + self.make_velmod_file(velmod_file, vp_velest, vs_velest, dep_velest) export_info = outdir / "export_info.json" export_info.write_text(self.model_dump_json(indent=2)) - print("Done!") return outdir - - @staticmethod - def export_phase( - event: EventDetection, - fp:TextIOWrapper, - phaseP:PhaseDescription, - phaseS:PhaseDescription, - tdiff:float, - pp:float, - ps:float, - ) -> int: - year=int(str(event.time.year)[-2::]) - # export phase catalog into txt file (for velest, format: ised=1) - # tdiff: time diff tobs-tmodel threshold P,S - # pp,ps : minimun probability threshold P,S - # there is quality weigth (0-4, 0 is best) in velest, here use probability to define it - month=event.time.month - day=event.time.day - hour=event.time.hour - min=event.time.minute - sec=event.time.second+event.time.microsecond/1000000 - lat=event.effective_lat - lon=event.effective_lon - dep=event.depth/1000 - if event.magnitude !=None: - mag=event.magnitude.average - else: - mag=0.0 - if lat<0: - vsn='S' - lat=abs(lat) + + def export_phases_slim( + self, + outfile: Path, + event: EventDetection, + observed_arrivals: list[tuple[Receiver, PhaseDetection]], + ): + mag = event.magnitude.average if event.magnitude is not None else 0.0 + lat = event.effective_lat + lon = event.effective_lon + if lat < 0: + vsn = "S" + lat = abs(lat) else: - vsn='N' - if lon<0: - vew='W' - lon=abs(lon) + vsn = "N" + if lon < 0: + vew = "W" + lon = abs(lon) else: - vew='E' - fp.write("%2d%2d%2d %2d%2d %5.2f %7.4f%s %8.4f%s %7.2f %5.2f\n"%(year,month,day,hour,min,sec,lat,vsn,lon,vew,dep,mag)) - countP=0 - countS=0 - for receiver in event.receivers.receivers: - station=receiver.station - if receiver.phase_arrivals[phaseP].observed !=None and receiver.phase_arrivals[phaseP].observed.detection_value >=pp: - tpick=receiver.phase_arrivals[phaseP].observed.time-event.time - tpick=tpick.total_seconds() - dt=receiver.phase_arrivals[phaseP].traveltime_delay.total_seconds() - if abs(dt)>tdiff: - continue - if receiver.phase_arrivals[phaseP].observed.detection_value <0.4: - iwt=3 - elif receiver.phase_arrivals[phaseP].observed.detection_value <0.6: - iwt=2 - elif receiver.phase_arrivals[phaseP].observed.detection_value <0.8: - iwt=1 + vew = "E" + with outfile.open("a") as file: + file.write( + f"{event.time.strftime('%y%m%d %H%M')} {event.time.second:2d}.{str(event.time.microsecond)[:2]} {lat:7.4f}{vsn:1s} {lon:8.4f}{vew:1s} {event.depth/1000:7.2f} {mag:5.2f}\n" + ) + count_p = 0 + count_s = 0 + for rec, dectection in observed_arrivals: + if dectection.observed.detection_value < 0.4: + quality_weight = 3 + elif dectection.observed.detection_value < 0.6: + quality_weight = 2 + elif dectection.observed.detection_value < 0.8: + quality_weight = 1 else: - iwt=0 - phase='P' - fp.write(" %-6s %-1s %1d %6.2f\n"%(station,phase,iwt,tpick)) - countP+=1 - if receiver.phase_arrivals[phaseS].observed !=None and receiver.phase_arrivals[phaseS].observed.detection_value >=ps: - tpick=receiver.phase_arrivals[phaseS].observed.time-event.time - tpick=tpick.total_seconds() - dt=receiver.phase_arrivals[phaseS].traveltime_delay.total_seconds() - if abs(dt)>tdiff: - continue - if receiver.phase_arrivals[phaseS].observed.detection_value <0.4: - iwt=3 - elif receiver.phase_arrivals[phaseS].observed.detection_value <0.6: - iwt=2 - elif receiver.phase_arrivals[phaseS].observed.detection_value <0.8: - iwt=1 + quality_weight = 0 + if dectection.phase.endswith("P"): + phase = "P" + count_p += 1 else: - iwt=0 - phase='S' - fp.write(" %-6s %-1s %1d %6.2f\n"%(station,phase,iwt,tpick)) - countS+=1 - if countP == 0 and countS == 0 : - print("No good phases obesered for event%s, lower pick probability or increase pick confidence"%event.time) - fp.write("\n") - return countP , countS + phase = "S" + count_s += 1 + traveltime = (dectection.observed.time - event.time).total_seconds() + file.write( + f" {rec.station:6s} {phase:1s} {quality_weight:1d} {traveltime:7.2f}\n" + ) + file.write("\n") + if count_p == 0 and count_s == 0: + logging.warning("Warning:No phases obesered for event{event.time}, removed") + with outfile.open("r") as file: + lines = file.readlines() + with outfile.open("w") as file: + file.writelines(lines[:-2]) - @staticmethod - def export_station( - stations: list[Station], - filename:Path - )-> None: - fpout=open(filename,'w') - fpout.write("(a6,f7.4,a1,1x,f8.4,a1,1x,i4,1x,i1,1x,i3,1x,f5.2,2x,f5.2)\n") - p2=1 - p3=1 - v1=0.00 - v2=0.00 - for station in stations: - lat=station.lat - lon=station.lon - sta=station.station - elev=station.elevation - if lat<0: - vsn='S' - lat=abs(lat) - else: - vsn='N' - if lon<0: - vew='W' - lon=abs(lon) - else: - vew='E' - fpout.write("%-6s%7.4f%1s %8.4f%s %4d %1d %3d %5.2f %5.2f\n"%(sta,lat,vsn,lon,vew,elev,p2,p3,v1,v2)) - p3+=1 - fpout.close() + return count_p, count_s @staticmethod - def make_control_file(filename, params,isingle,modname,phasefile,mainoutfile,outcheckfile,finalcnv,stacorfile): - - if isingle == 1: - ittmax = 99 - invertratio = 0 - else: - ittmax = 9 - invertratio = 3 - fp=open(filename,'w') - fp.write("velest parameters are below, please modify according to their documents\n") - #*** olat olon icoordsystem zshift itrial ztrial ised - fp.write("%s %s 0 0.0 0 0.00 1\n"%(params["reflat"],params["reflon"])) - #*** neqs nshot rotate - fp.write("%s 0 0.0\n"%params["neqs"]) - #*** isingle iresolcalc - fp.write("%s 0\n"%isingle) - #*** dmax itopo zmin veladj zadj lowveloclay - fp.write("%s 0 %s 0.20 5.00 %s\n"%(params["distmax"],params["zmin"],params["lowvelocity"] )) - #*** nsp swtfac vpvs nmod - fp.write( "2 0.75 1.650 1\n") - #*** othet xythet zthet vthet stathet - fp.write( "0.01 0.01 0.01 %s %s\n"%(params["vthet"],params["stathet"])) - #*** nsinv nshcor nshfix iuseelev iusestacorr - fp.write( "1 0 0 %s %s\n"%(params["iuseelev"],params["iusestacorr"])) - #*** iturbo icnvout istaout ismpout - fp.write( "1 1 2 0\n") - #*** irayout idrvout ialeout idspout irflout irfrout iresout - fp.write( "0 0 0 0 0 0 0\n") - #*** delmin ittmax invertratio - fp.write( "0.001 %s %s\n"%(ittmax,invertratio)) - #*** Modelfile: - fp.write( "%s\n"%modname) - #*** Stationfile: - fp.write( "stations_velest.sta\n") - #*** Seismofile: - fp.write( " \n") - #*** File with region names: - fp.write( "regionsnamen.dat\n") - #*** File with region coordinates: - fp.write( "regionskoord.dat\n") - #*** File #1 with topo data: - fp.write( " \n") - #*** File #2 with topo data: - fp.write( " \n") - #*** File with Earthquake data (phase): - fp.write( "%s\n"%phasefile) - #*** File with Shot data: - fp.write( " \n") - #*** Main print output file: - fp.write( "%s\n"%mainoutfile) - #*** File with single event locations: - fp.write( "%s\n"%outcheckfile) - #*** File with final hypocenters in *.cnv format: - fp.write( "%s\n"%finalcnv) - #*** File with new station corrections: - fp.write( "%s\n"%stacorfile) - fp.close() + def export_station(stations: list[Station], filename: Path) -> None: + with filename.open("w") as fpout: + fpout.write("(a6,f7.4,a1,1x,f8.4,a1,1x,i4,1x,i1,1x,i3,1x,f5.2,2x,f5.2)\n") + station_index = 1 + for station in stations: + lat = station.lat + lon = station.lon + sta = station.station + elev = station.elevation + if lat < 0: + vsn = "S" + lat = abs(lat) + else: + vsn = "N" + if lon < 0: + vew = "W" + lon = abs(lon) + else: + vew = "E" + fpout.write( + f"{sta:6s}{lat:7.4f}{vsn} {lon:8.4f}{vew} {int(elev):4d} 1 {index:3d} 0.00 0.00\n" + ) + station_index += 1 + fpout.write("\n") @staticmethod - def make_velmod_file(modname,vp,vs,dep ): - nlayer=len(dep) + def make_velmod_file(modname: Path, vp: list, vs: list, dep: list): + nlayer = len(dep) vdamp = 1.0 - fp=open(modname,'w') - fp.write("initial 1D-model for velest\n") - # the second line - indicate the number of layers for Vp - fp.write("%3d vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n"%nlayer) - #vp model - for i,v in enumerate(vp): - fp.write("%5.2f %7.2f %7.3f\n"%(v,dep[i],vdamp)) - #vs model - fp.write("%3d\n"%nlayer) - for i,v in enumerate(vs): - fp.write("%5.2f %7.2f %7.3f\n"%(v,dep[i],vdamp)) - fp.close() - \ No newline at end of file + with modname.open("w") as fp: + fp.write("initial 1D-model for velest\n") + # the second line - indicate the number of layers for Vp + fp.write( + f"{nlayer} vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n" + ) + # vp model + for i, v in enumerate(vp): + fp.write(f"{v:5.2f} {dep[i]:7.2f} {vdamp:7.3f}\n") + # vs model + fp.write("%3d\n" % nlayer) + for i, v in enumerate(vs): + fp.write(f"{v:5.2f} {dep[i]:7.2f} {vdamp:7.3f}\n") From 4fb8c883ed2b1cfeaf1739693f3faa6a31097409 Mon Sep 17 00:00:00 2001 From: Hao Date: Thu, 11 Jul 2024 09:29:54 +0000 Subject: [PATCH 06/10] add velest exporter --- src/qseek/exporters/velest.py | 301 ++++++++++++++++++++++++++++++++++ 1 file changed, 301 insertions(+) create mode 100644 src/qseek/exporters/velest.py diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py new file mode 100644 index 00000000..7492de43 --- /dev/null +++ b/src/qseek/exporters/velest.py @@ -0,0 +1,301 @@ +from __future__ import annotations + +import logging +from pathlib import Path +from typing import NamedTuple + +import rich +from rich.prompt import FloatPrompt + +from qseek.exporters.base import Exporter +from qseek.models.detection import EventDetection, PhaseDetection, Receiver +from qseek.models.station import Station +from qseek.search import Search + +logger = logging.getLogger(__name__) + +CONTROL_FILE_TPL = """velest parameters are below, please modify according to their documents +{ref_lat} {ref_lon} 0 0.0 0 0.00 1 +{n_earthquakes} 0 0.0 +{isingle:1d} 0 +{max_distance_station} 0 {min_depth} 0.20 5.00 {use_station_correction:1d} +2 0.75 {vp_vs_ratio} 1 +0.01 0.01 0.01 {velocity_damping} {station_correction_damping} +1 0 0 {use_elevation:1d} {use_station_correction:1d} +1 1 2 0 +0 0 0 0 0 0 0 +0.001 {iteration_number} {invertratio} +{model_file} +stations_velest.sta + +regionsnamen.dat +regionskoord.dat + + +{phase_file} + +{mainout_file} +{outcheck_file} +{finalcnv_file} +{stacorrection_file} +""" + + +class VelestControlFile(NamedTuple): + ref_lat: float + ref_lon: float # should be negative for East + n_earthquakes: int + isingle: bool + max_distance_station: float + min_depth: float + allow_low_velocity: bool + vp_vs_ratio: float + velocity_damping: float # Damping parameter for the velocity + station_correction_damping: float # Damping parameter for the station + use_elevation: bool + use_station_correction: bool + iteration_number: int + invertratio: int + model_file: str + phase_file: str + mainout_file: str + outcheck_file: str + finalcnv_file: str + stacorrection_file: str + + def write_config_file(self, file: Path): + with file.open("w") as fp: + fp.write(CONTROL_FILE_TPL.format(**self._asdict())) + + +class Velest(Exporter): + """Crate a VELEST project folder for 1D velocity model estimation.""" + + min_pick_semblance: float = 0.2 + min_receivers_number: int = 10 + min_p_phase_confidence: float = 0.3 + min_s_phase_confidence: float = 0.3 + max_traveltime_delay: float = 2.5 + n_picks_p: int = 0 + n_picks_s: int = 0 + n_events: int = 0 + + async def export(self, rundir: Path, outdir: Path) -> Path: + rich.print("Exporting qseek search to VELEST project folder") + min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.2) + min_receivers_number = FloatPrompt.ask( + "Minimum number of receivers (P phase)", default=10 + ) + min_p_phase_confidence = FloatPrompt.ask( + "Minimum pick probability for P phase", default=0.3 + ) + min_s_phase_confidence = FloatPrompt.ask( + "Minimum pick probability for S phase", default=0.3 + ) + max_traveltime_delay = FloatPrompt.ask( + "Maximum difference between theoretical and observed arrival", default=2.5 + ) + self.min_pick_semblance = min_pick_semblance + self.min_receivers_number = min_receivers_number + self.min_p_phase_confidence = min_p_phase_confidence + self.min_s_phase_confidence = min_s_phase_confidence + self.max_traveltime_delay = max_traveltime_delay + + outdir.mkdir() + search = Search.load_rundir(rundir) + phases = search.image_functions.get_phases() + for phase in phases: + if "P" in phase: + phase_p = phase + if "S" in phase: + phase_s = phase + + catalog = search.catalog + + # export station file + stations = search.stations.stations + station_file = outdir / "stations_velest.sta" + self.export_station(stations=stations, filename=station_file) + + # export phase file + phase_file = outdir / "phase_velest.pha" + n_earthquakes = 0 + for event in catalog: + if event.semblance < min_pick_semblance: + continue + if event.receivers.n_observations(phase_p) < min_receivers_number: + continue + + observed_arrivals: list[tuple[Receiver, PhaseDetection]] = [] + + for receiver in event.receivers: + for _phase, detection in receiver.phase_arrivals.items(): + if detection.observed is None: + continue + observed = detection.observed + if ( + detection.phase == phase_p + and observed.detection_value <= min_p_phase_confidence + ): + continue + if ( + detection.phase == phase_s + and observed.detection_value <= min_s_phase_confidence + ): + continue + if ( + detection.traveltime_delay.total_seconds() + > max_traveltime_delay + ): + continue + observed_arrivals.append((receiver, detection)) + + countp, counts = self.export_phases_slim( + phase_file, event, observed_arrivals + ) + self.n_picks_p += countp + self.n_picks_s += counts + n_earthquakes += 1 + self.n_events = n_earthquakes + + # export control file + control_file = outdir / "velest.cmn" + control_file_parameters = VelestControlFile( + ref_lat=search.octree.location.lat, + ref_lon=-search.octree.location.lon, + n_earthquakes=n_earthquakes, + max_distance_station=200, + min_depth=-0.2, + allow_low_velocity=False, + velocity_damping=1.0, + station_correction_damping=0.1, + use_elevation=False, + use_station_correction=False, + model_file="model.mod", + phase_file="phase_velest.pha", + mainout_file="main.out", + outcheck_file="log.out", + finalcnv_file="final.cnv", + stacorrection_file="stacor.dat", + isingle=False, + vp_vs_ratio=1.65, + iteration_number=99, + invertratio=0, + ) + control_file_parameters.write_config_file(control_file) + # export velocity model file + dep = search.ray_tracers.root[0].earthmodel.layered_model.profile("z") + vp = search.ray_tracers.root[0].earthmodel.layered_model.profile("vp") + vs = search.ray_tracers.root[0].earthmodel.layered_model.profile("vs") + dep_velest = [] + vp_velest = [] + vs_velest = [] + for i, d in enumerate(dep): + if float(d) / 1000 not in dep_velest: + dep_velest.append(float(d) / 1000) + vp_velest.append(float(vp[i]) / 1000) + vs_velest.append(float(vs[i]) / 1000) + velmod_file = outdir / "model.mod" + self.make_velmod_file(velmod_file, vp_velest, vs_velest, dep_velest) + + export_info = outdir / "export_info.json" + export_info.write_text(self.model_dump_json(indent=2)) + return outdir + + def export_phases_slim( + self, + outfile: Path, + event: EventDetection, + observed_arrivals: list[tuple[Receiver, PhaseDetection]], + ): + mag = event.magnitude.average if event.magnitude is not None else 0.0 + lat = event.effective_lat + lon = event.effective_lon + if lat < 0: + vsn = "S" + lat = abs(lat) + else: + vsn = "N" + if lon < 0: + vew = "W" + lon = abs(lon) + else: + vew = "E" + with outfile.open("a") as file: + file.write( + f"{event.time.strftime('%y%m%d %H%M')} {event.time.second:2d}.{str(event.time.microsecond)[:2]} {lat:7.4f}{vsn:1s} {lon:8.4f}{vew:1s} {event.depth/1000:7.2f} {mag:5.2f}\n" + ) + count_p = 0 + count_s = 0 + for rec, dectection in observed_arrivals: + if dectection.observed.detection_value < 0.4: + quality_weight = 3 + elif dectection.observed.detection_value < 0.6: + quality_weight = 2 + elif dectection.observed.detection_value < 0.8: + quality_weight = 1 + else: + quality_weight = 0 + if dectection.phase.endswith("P"): + phase = "P" + count_p += 1 + else: + phase = "S" + count_s += 1 + traveltime = (dectection.observed.time - event.time).total_seconds() + file.write( + f" {rec.station:6s} {phase:1s} {quality_weight:1d} {traveltime:7.2f}\n" + ) + file.write("\n") + if count_p == 0 and count_s == 0: + logging.warning("Warning:No phases obesered for event{event.time}, removed") + with outfile.open("r") as file: + lines = file.readlines() + with outfile.open("w") as file: + file.writelines(lines[:-2]) + + return count_p, count_s + + @staticmethod + def export_station(stations: list[Station], filename: Path) -> None: + with filename.open("w") as fpout: + fpout.write("(a6,f7.4,a1,1x,f8.4,a1,1x,i4,1x,i1,1x,i3,1x,f5.2,2x,f5.2)\n") + station_index = 1 + for station in stations: + lat = station.lat + lon = station.lon + sta = station.station + elev = station.elevation + if lat < 0: + vsn = "S" + lat = abs(lat) + else: + vsn = "N" + if lon < 0: + vew = "W" + lon = abs(lon) + else: + vew = "E" + fpout.write( + f"{sta:6s}{lat:7.4f}{vsn} {lon:8.4f}{vew} {int(elev):4d} 1 {station_index:3d} 0.00 0.00\n" + ) + station_index += 1 + fpout.write("\n") + + @staticmethod + def make_velmod_file(modname: Path, vp: list, vs: list, dep: list): + nlayer = len(dep) + vdamp = 1.0 + with modname.open("w") as fp: + fp.write("initial 1D-model for velest\n") + # the second line - indicate the number of layers for Vp + fp.write( + f"{nlayer} vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n" + ) + # vp model + for i, v in enumerate(vp): + fp.write(f"{v:5.2f} {dep[i]:7.2f} {vdamp:7.3f}\n") + # vs model + fp.write("%3d\n" % nlayer) + for i, v in enumerate(vs): + fp.write(f"{v:5.2f} {dep[i]:7.2f} {vdamp:7.3f}\n") From 769cd6f1cba4c65d22a654381c73b2b0ab9d7b82 Mon Sep 17 00:00:00 2001 From: Hao Date: Fri, 12 Jul 2024 12:41:20 +0000 Subject: [PATCH 07/10] modified according to reviews --- src/qseek/exporters/velest.py | 112 +++++++++++++--------------------- 1 file changed, 41 insertions(+), 71 deletions(-) diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index 7492de43..28fb7e55 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -5,11 +5,11 @@ from typing import NamedTuple import rich -from rich.prompt import FloatPrompt +from rich.prompt import FloatPrompt, IntPrompt from qseek.exporters.base import Exporter from qseek.models.detection import EventDetection, PhaseDetection, Receiver -from qseek.models.station import Station +from qseek.models.station import Location, Station from qseek.search import Search logger = logging.getLogger(__name__) @@ -45,23 +45,23 @@ class VelestControlFile(NamedTuple): ref_lat: float ref_lon: float # should be negative for East n_earthquakes: int - isingle: bool - max_distance_station: float - min_depth: float - allow_low_velocity: bool - vp_vs_ratio: float - velocity_damping: float # Damping parameter for the velocity - station_correction_damping: float # Damping parameter for the station - use_elevation: bool - use_station_correction: bool - iteration_number: int - invertratio: int - model_file: str - phase_file: str - mainout_file: str - outcheck_file: str - finalcnv_file: str - stacorrection_file: str + isingle: bool = True + min_depth: float = -0.2 + vp_vs_ratio: float = 1.65 + iteration_number: int = 99 + invertratio: int = 0 + model_file: str = "model.mod" + phase_file: str = "phase_velest.pha" + mainout_file: str = "main.out" + outcheck_file: str = "log.out" + finalcnv_file: str = "final.cnv" + stacorrection_file: str = "stacor.dat" + velocity_damping: float = 1.0 # Damping parameter for the velocity + station_correction_damping: float = 0.1 # Damping parameter for the station + max_distance_station: float = 200.0 + use_elevation: bool = False + use_station_correction: bool = False + allow_low_velocity: bool = False def write_config_file(self, file: Path): with file.open("w") as fp: @@ -83,7 +83,7 @@ class Velest(Exporter): async def export(self, rundir: Path, outdir: Path) -> Path: rich.print("Exporting qseek search to VELEST project folder") min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.2) - min_receivers_number = FloatPrompt.ask( + min_receivers_number = IntPrompt.ask( "Minimum number of receivers (P phase)", default=10 ) min_p_phase_confidence = FloatPrompt.ask( @@ -164,23 +164,6 @@ async def export(self, rundir: Path, outdir: Path) -> Path: ref_lat=search.octree.location.lat, ref_lon=-search.octree.location.lon, n_earthquakes=n_earthquakes, - max_distance_station=200, - min_depth=-0.2, - allow_low_velocity=False, - velocity_damping=1.0, - station_correction_damping=0.1, - use_elevation=False, - use_station_correction=False, - model_file="model.mod", - phase_file="phase_velest.pha", - mainout_file="main.out", - outcheck_file="log.out", - finalcnv_file="final.cnv", - stacorrection_file="stacor.dat", - isingle=False, - vp_vs_ratio=1.65, - iteration_number=99, - invertratio=0, ) control_file_parameters.write_config_file(control_file) # export velocity model file @@ -202,6 +185,17 @@ async def export(self, rundir: Path, outdir: Path) -> Path: export_info.write_text(self.model_dump_json(indent=2)) return outdir + def velest_location(self, location: Location) -> str: + if location.effective_lat < 0: + velest_lat = f"{location.effective_lat:7.4f}S" + else: + velest_lat = f"{location.effective_lat:7.4f}N" + if location.effective_lon < 0: + velest_lon = f"{location.effective_lon:8.4f}W" + else: + velest_lon = f"{location.effective_lon:8.4f}E" + return velest_lat, velest_lon + def export_phases_slim( self, outfile: Path, @@ -209,21 +203,11 @@ def export_phases_slim( observed_arrivals: list[tuple[Receiver, PhaseDetection]], ): mag = event.magnitude.average if event.magnitude is not None else 0.0 - lat = event.effective_lat - lon = event.effective_lon - if lat < 0: - vsn = "S" - lat = abs(lat) - else: - vsn = "N" - if lon < 0: - vew = "W" - lon = abs(lon) - else: - vew = "E" + lat, lon = self.velest_location(event) with outfile.open("a") as file: file.write( - f"{event.time.strftime('%y%m%d %H%M')} {event.time.second:2d}.{str(event.time.microsecond)[:2]} {lat:7.4f}{vsn:1s} {lon:8.4f}{vew:1s} {event.depth/1000:7.2f} {mag:5.2f}\n" + f"{event.time:%y%m%d %H%M %S.%f}"[:-4] + + f" {lat} {lon} {event.depth/1000:7.2f} {mag:5.2f}\n" ) count_p = 0 count_s = 0 @@ -256,28 +240,14 @@ def export_phases_slim( return count_p, count_s - @staticmethod - def export_station(stations: list[Station], filename: Path) -> None: + def export_station(self, stations: list[Station], filename: Path) -> None: with filename.open("w") as fpout: fpout.write("(a6,f7.4,a1,1x,f8.4,a1,1x,i4,1x,i1,1x,i3,1x,f5.2,2x,f5.2)\n") station_index = 1 for station in stations: - lat = station.lat - lon = station.lon - sta = station.station - elev = station.elevation - if lat < 0: - vsn = "S" - lat = abs(lat) - else: - vsn = "N" - if lon < 0: - vew = "W" - lon = abs(lon) - else: - vew = "E" + lat, lon = self.velest_location(station) fpout.write( - f"{sta:6s}{lat:7.4f}{vsn} {lon:8.4f}{vew} {int(elev):4d} 1 {station_index:3d} 0.00 0.00\n" + f"{station.station:6s}{lat} {lon} {int(station.elevation):4d} 1 {station_index:3d} 0.00 0.00\n" ) station_index += 1 fpout.write("\n") @@ -293,9 +263,9 @@ def make_velmod_file(modname: Path, vp: list, vs: list, dep: list): f"{nlayer} vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n" ) # vp model - for i, v in enumerate(vp): - fp.write(f"{v:5.2f} {dep[i]:7.2f} {vdamp:7.3f}\n") + for vel, depth in zip(vp, dep, strict=False): + fp.write(f"{vel:5.2f} {depth:7.2f} {vdamp:7.3f}\n") # vs model fp.write("%3d\n" % nlayer) - for i, v in enumerate(vs): - fp.write(f"{v:5.2f} {dep[i]:7.2f} {vdamp:7.3f}\n") + for vel, depth in zip(vs, dep, strict=False): + fp.write(f"{vel:5.2f} {depth:7.2f} {vdamp:7.3f}\n") From f7b435e58d3745d09a329fc48de5add29ddf6924 Mon Sep 17 00:00:00 2001 From: Hao Date: Mon, 15 Jul 2024 07:39:51 +0000 Subject: [PATCH 08/10] modified according to reviews --- src/qseek/exporters/velest.py | 89 ++++++++++++++++++----------------- 1 file changed, 47 insertions(+), 42 deletions(-) diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index 28fb7e55..7f597bf0 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -65,27 +65,31 @@ class VelestControlFile(NamedTuple): def write_config_file(self, file: Path): with file.open("w") as fp: - fp.write(CONTROL_FILE_TPL.format(**self._asdict())) + fp.write_text(CONTROL_FILE_TPL.format(**self._asdict())) class Velest(Exporter): """Crate a VELEST project folder for 1D velocity model estimation.""" - min_pick_semblance: float = 0.2 + min_event_semblance: float = 0.2 min_receivers_number: int = 10 min_p_phase_confidence: float = 0.3 min_s_phase_confidence: float = 0.3 max_traveltime_delay: float = 2.5 + distance_border: float = 500.0 n_picks_p: int = 0 n_picks_s: int = 0 n_events: int = 0 async def export(self, rundir: Path, outdir: Path) -> Path: rich.print("Exporting qseek search to VELEST project folder") - min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.2) + min_event_semblance = self.min_event_semblance min_receivers_number = IntPrompt.ask( "Minimum number of receivers (P phase)", default=10 ) + min_distance_to_border = FloatPrompt.ask( + "Minimum distance to border(m)", default=500.0 + ) min_p_phase_confidence = FloatPrompt.ask( "Minimum pick probability for P phase", default=0.3 ) @@ -95,7 +99,6 @@ async def export(self, rundir: Path, outdir: Path) -> Path: max_traveltime_delay = FloatPrompt.ask( "Maximum difference between theoretical and observed arrival", default=2.5 ) - self.min_pick_semblance = min_pick_semblance self.min_receivers_number = min_receivers_number self.min_p_phase_confidence = min_p_phase_confidence self.min_s_phase_confidence = min_s_phase_confidence @@ -121,10 +124,12 @@ async def export(self, rundir: Path, outdir: Path) -> Path: phase_file = outdir / "phase_velest.pha" n_earthquakes = 0 for event in catalog: - if event.semblance < min_pick_semblance: + if event.semblance < min_event_semblance: continue if event.receivers.n_observations(phase_p) < min_receivers_number: continue + if event.distance_border < min_distance_to_border: + continue observed_arrivals: list[tuple[Receiver, PhaseDetection]] = [] @@ -173,29 +178,19 @@ async def export(self, rundir: Path, outdir: Path) -> Path: dep_velest = [] vp_velest = [] vs_velest = [] - for i, d in enumerate(dep): - if float(d) / 1000 not in dep_velest: - dep_velest.append(float(d) / 1000) - vp_velest.append(float(vp[i]) / 1000) - vs_velest.append(float(vs[i]) / 1000) + km = 1000.0 + for d, vpi, vsi in zip(dep, vp, vs, strict=False): + if float(d) / km not in dep_velest: + dep_velest.append(float(d) / km) + vp_velest.append(float(vpi) / km) + vs_velest.append(float(vsi) / km) velmod_file = outdir / "model.mod" - self.make_velmod_file(velmod_file, vp_velest, vs_velest, dep_velest) + make_velmod_file(velmod_file, vp_velest, vs_velest, dep_velest) export_info = outdir / "export_info.json" export_info.write_text(self.model_dump_json(indent=2)) return outdir - def velest_location(self, location: Location) -> str: - if location.effective_lat < 0: - velest_lat = f"{location.effective_lat:7.4f}S" - else: - velest_lat = f"{location.effective_lat:7.4f}N" - if location.effective_lon < 0: - velest_lon = f"{location.effective_lon:8.4f}W" - else: - velest_lon = f"{location.effective_lon:8.4f}E" - return velest_lat, velest_lon - def export_phases_slim( self, outfile: Path, @@ -203,7 +198,7 @@ def export_phases_slim( observed_arrivals: list[tuple[Receiver, PhaseDetection]], ): mag = event.magnitude.average if event.magnitude is not None else 0.0 - lat, lon = self.velest_location(event) + lat, lon = velest_location(event) with outfile.open("a") as file: file.write( f"{event.time:%y%m%d %H%M %S.%f}"[:-4] @@ -232,7 +227,7 @@ def export_phases_slim( ) file.write("\n") if count_p == 0 and count_s == 0: - logging.warning("Warning:No phases obesered for event{event.time}, removed") + logger.warning("Removing event {event.time}: No phases observed") with outfile.open("r") as file: lines = file.readlines() with outfile.open("w") as file: @@ -245,27 +240,37 @@ def export_station(self, stations: list[Station], filename: Path) -> None: fpout.write("(a6,f7.4,a1,1x,f8.4,a1,1x,i4,1x,i1,1x,i3,1x,f5.2,2x,f5.2)\n") station_index = 1 for station in stations: - lat, lon = self.velest_location(station) + lat, lon = velest_location(station) fpout.write( f"{station.station:6s}{lat} {lon} {int(station.elevation):4d} 1 {station_index:3d} 0.00 0.00\n" ) station_index += 1 fpout.write("\n") - @staticmethod - def make_velmod_file(modname: Path, vp: list, vs: list, dep: list): - nlayer = len(dep) - vdamp = 1.0 - with modname.open("w") as fp: - fp.write("initial 1D-model for velest\n") - # the second line - indicate the number of layers for Vp - fp.write( - f"{nlayer} vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n" - ) - # vp model - for vel, depth in zip(vp, dep, strict=False): - fp.write(f"{vel:5.2f} {depth:7.2f} {vdamp:7.3f}\n") - # vs model - fp.write("%3d\n" % nlayer) - for vel, depth in zip(vs, dep, strict=False): - fp.write(f"{vel:5.2f} {depth:7.2f} {vdamp:7.3f}\n") + +def velest_location(location: Location) -> str: + if location.effective_lat < 0: + velest_lat = f"{location.effective_lat:7.4f}S" + else: + velest_lat = f"{location.effective_lat:7.4f}N" + if location.effective_lon < 0: + velest_lon = f"{location.effective_lon:8.4f}W" + else: + velest_lon = f"{location.effective_lon:8.4f}E" + return velest_lat, velest_lon + + +def make_velmod_file(modname: Path, vp: list[float], vs: list[float], dep: list[float]): + nlayer = len(dep) + vdamp = 1.0 + with modname.open("w") as fp: + fp.write("initial 1D-model for velest\n") + # the second line - indicate the number of layers for Vp + fp.write(f"{nlayer} vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n") + # vp model + for vel, depth in zip(vp, dep, strict=False): + fp.write(f"{vel:5.2f} {depth:7.2f} {vdamp:7.3f}\n") + # vs model + fp.write("%3d\n" % nlayer) + for vel, depth in zip(vs, dep, strict=False): + fp.write(f"{vel:5.2f} {depth:7.2f} {vdamp:7.3f}\n") From 6318e00b241ffb53afb61c307f25d4d6cc7b10e2 Mon Sep 17 00:00:00 2001 From: Hao Date: Tue, 16 Jul 2024 14:21:00 +0000 Subject: [PATCH 09/10] modified according to reviews --- src/qseek/exporters/velest.py | 63 +++++++++++++++-------------------- 1 file changed, 27 insertions(+), 36 deletions(-) diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index 7f597bf0..442429e7 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -4,6 +4,7 @@ from pathlib import Path from typing import NamedTuple +import numpy as np import rich from rich.prompt import FloatPrompt, IntPrompt @@ -64,8 +65,7 @@ class VelestControlFile(NamedTuple): allow_low_velocity: bool = False def write_config_file(self, file: Path): - with file.open("w") as fp: - fp.write_text(CONTROL_FILE_TPL.format(**self._asdict())) + file.write_text(CONTROL_FILE_TPL.format(**self._asdict())) class Velest(Exporter): @@ -96,13 +96,10 @@ async def export(self, rundir: Path, outdir: Path) -> Path: min_s_phase_confidence = FloatPrompt.ask( "Minimum pick probability for S phase", default=0.3 ) - max_traveltime_delay = FloatPrompt.ask( - "Maximum difference between theoretical and observed arrival", default=2.5 - ) + max_traveltime_delay = self.max_traveltime_delay self.min_receivers_number = min_receivers_number self.min_p_phase_confidence = min_p_phase_confidence self.min_s_phase_confidence = min_s_phase_confidence - self.max_traveltime_delay = max_traveltime_delay outdir.mkdir() search = Search.load_rundir(rundir) @@ -199,39 +196,33 @@ def export_phases_slim( ): mag = event.magnitude.average if event.magnitude is not None else 0.0 lat, lon = velest_location(event) - with outfile.open("a") as file: - file.write( - f"{event.time:%y%m%d %H%M %S.%f}"[:-4] - + f" {lat} {lon} {event.depth/1000:7.2f} {mag:5.2f}\n" - ) - count_p = 0 - count_s = 0 - for rec, dectection in observed_arrivals: - if dectection.observed.detection_value < 0.4: - quality_weight = 3 - elif dectection.observed.detection_value < 0.6: - quality_weight = 2 - elif dectection.observed.detection_value < 0.8: - quality_weight = 1 - else: - quality_weight = 0 - if dectection.phase.endswith("P"): - phase = "P" - count_p += 1 - else: - phase = "S" - count_s += 1 - traveltime = (dectection.observed.time - event.time).total_seconds() - file.write( - f" {rec.station:6s} {phase:1s} {quality_weight:1d} {traveltime:7.2f}\n" + write_out = ( + f"{event.time:%y%m%d %H%M %S.%f}"[:-4] + + f" {lat} {lon} {event.depth/1000:7.2f} {mag:5.2f}\n" + ) + count_p = 0 + count_s = 0 + for rec, dectection in observed_arrivals: + quality_weight = ( + np.digitize( + dectection.observed.detection_value, [1.0, 0.8, 0.6, 0.4, 0.0] ) - file.write("\n") + - 1 + ) + if dectection.phase.endswith("P"): + phase = "P" + count_p += 1 + else: + phase = "S" + count_s += 1 + traveltime = (dectection.observed.time - event.time).total_seconds() + write_out += f" {rec.station:6s} {phase:1s} {quality_weight:1d} {traveltime:7.2f}\n" + write_out += "\n" if count_p == 0 and count_s == 0: logger.warning("Removing event {event.time}: No phases observed") - with outfile.open("r") as file: - lines = file.readlines() - with outfile.open("w") as file: - file.writelines(lines[:-2]) + if count_p or count_s: + with outfile.open("a") as file: + file.write(write_out) return count_p, count_s From 9a64a5b8e0369d51dffe9df8a258c325417ce2e2 Mon Sep 17 00:00:00 2001 From: Hao Date: Tue, 16 Jul 2024 14:59:30 +0000 Subject: [PATCH 10/10] modified exporter velest --- src/qseek/exporters/velest.py | 94 ++++++++++++++++++++--------------- 1 file changed, 55 insertions(+), 39 deletions(-) diff --git a/src/qseek/exporters/velest.py b/src/qseek/exporters/velest.py index 442429e7..4db07856 100644 --- a/src/qseek/exporters/velest.py +++ b/src/qseek/exporters/velest.py @@ -14,6 +14,9 @@ from qseek.search import Search logger = logging.getLogger(__name__) +KM = 1000.0 + +CONFIDENCE_QUALITY_BINS = [1.0, 0.8, 0.6, 0.4, 0.0] CONTROL_FILE_TPL = """velest parameters are below, please modify according to their documents {ref_lat} {ref_lon} 0 0.0 0 0.00 1 @@ -83,31 +86,38 @@ class Velest(Exporter): async def export(self, rundir: Path, outdir: Path) -> Path: rich.print("Exporting qseek search to VELEST project folder") - min_event_semblance = self.min_event_semblance - min_receivers_number = IntPrompt.ask( - "Minimum number of receivers (P phase)", default=10 + self.min_event_semblance = FloatPrompt.ask( + "Minimum event semblance", + default=self.min_event_semblance, + ) + self.min_receivers_number = IntPrompt.ask( + "Minimum number of receivers (P phase)", + default=self.min_receivers_number, + ) + self.min_distance_to_border = FloatPrompt.ask( + "Minimum distance to border(m)", + default=self.distance_border, ) - min_distance_to_border = FloatPrompt.ask( - "Minimum distance to border(m)", default=500.0 + self.min_p_phase_confidence = FloatPrompt.ask( + "Minimum pick probability for P phase", + default=self.min_p_phase_confidence, ) - min_p_phase_confidence = FloatPrompt.ask( - "Minimum pick probability for P phase", default=0.3 + self.min_s_phase_confidence = FloatPrompt.ask( + "Minimum pick probability for S phase", + default=self.min_s_phase_confidence, ) - min_s_phase_confidence = FloatPrompt.ask( - "Minimum pick probability for S phase", default=0.3 + self.max_traveltime_delay = FloatPrompt.ask( + "Maximum travel time delay", + default=self.max_traveltime_delay, ) - max_traveltime_delay = self.max_traveltime_delay - self.min_receivers_number = min_receivers_number - self.min_p_phase_confidence = min_p_phase_confidence - self.min_s_phase_confidence = min_s_phase_confidence outdir.mkdir() search = Search.load_rundir(rundir) phases = search.image_functions.get_phases() for phase in phases: - if "P" in phase: + if phase.endswith("P"): phase_p = phase - if "S" in phase: + if phase.endswith("S"): phase_s = phase catalog = search.catalog @@ -121,11 +131,11 @@ async def export(self, rundir: Path, outdir: Path) -> Path: phase_file = outdir / "phase_velest.pha" n_earthquakes = 0 for event in catalog: - if event.semblance < min_event_semblance: + if event.semblance < self.min_event_semblance: continue - if event.receivers.n_observations(phase_p) < min_receivers_number: + if event.receivers.n_observations(phase_p) < self.min_receivers_number: continue - if event.distance_border < min_distance_to_border: + if event.distance_border < self.min_distance_to_border: continue observed_arrivals: list[tuple[Receiver, PhaseDetection]] = [] @@ -137,26 +147,26 @@ async def export(self, rundir: Path, outdir: Path) -> Path: observed = detection.observed if ( detection.phase == phase_p - and observed.detection_value <= min_p_phase_confidence + and observed.detection_value <= self.min_p_phase_confidence ): continue if ( detection.phase == phase_s - and observed.detection_value <= min_s_phase_confidence + and observed.detection_value <= self.min_s_phase_confidence ): continue if ( detection.traveltime_delay.total_seconds() - > max_traveltime_delay + > self.max_traveltime_delay ): continue observed_arrivals.append((receiver, detection)) - countp, counts = self.export_phases_slim( + count_p, count_s = self.export_phases_slim( phase_file, event, observed_arrivals ) - self.n_picks_p += countp - self.n_picks_s += counts + self.n_picks_p += count_p + self.n_picks_s += count_s n_earthquakes += 1 self.n_events = n_earthquakes @@ -175,12 +185,11 @@ async def export(self, rundir: Path, outdir: Path) -> Path: dep_velest = [] vp_velest = [] vs_velest = [] - km = 1000.0 - for d, vpi, vsi in zip(dep, vp, vs, strict=False): - if float(d) / km not in dep_velest: - dep_velest.append(float(d) / km) - vp_velest.append(float(vpi) / km) - vs_velest.append(float(vsi) / km) + for d, vpi, vsi in zip(dep, vp, vs, strict=True): + if float(d) / KM not in dep_velest: + dep_velest.append(float(d) / KM) + vp_velest.append(float(vpi) / KM) + vs_velest.append(float(vsi) / KM) velmod_file = outdir / "model.mod" make_velmod_file(velmod_file, vp_velest, vs_velest, dep_velest) @@ -205,7 +214,8 @@ def export_phases_slim( for rec, dectection in observed_arrivals: quality_weight = ( np.digitize( - dectection.observed.detection_value, [1.0, 0.8, 0.6, 0.4, 0.0] + dectection.observed.detection_value, + CONFIDENCE_QUALITY_BINS, ) - 1 ) @@ -218,11 +228,12 @@ def export_phases_slim( traveltime = (dectection.observed.time - event.time).total_seconds() write_out += f" {rec.station:6s} {phase:1s} {quality_weight:1d} {traveltime:7.2f}\n" write_out += "\n" - if count_p == 0 and count_s == 0: - logger.warning("Removing event {event.time}: No phases observed") + if count_p or count_s: with outfile.open("a") as file: file.write(write_out) + else: + logger.warning("Event {event.time}: No phases observed") return count_p, count_s @@ -251,17 +262,22 @@ def velest_location(location: Location) -> str: return velest_lat, velest_lon -def make_velmod_file(modname: Path, vp: list[float], vs: list[float], dep: list[float]): - nlayer = len(dep) +def make_velmod_file( + modname: Path, + velocity_p: list[float], + velocity_s: list[float], + depths: list[float], +): + nlayer = len(depths) vdamp = 1.0 with modname.open("w") as fp: fp.write("initial 1D-model for velest\n") # the second line - indicate the number of layers for Vp fp.write(f"{nlayer} vel,depth,vdamp,phase (f5.2,5x,f7.2,2x,f7.3,3x,a1)\n") - # vp model - for vel, depth in zip(vp, dep, strict=False): + + for vel, depth in zip(velocity_p, depths, strict=True): fp.write(f"{vel:5.2f} {depth:7.2f} {vdamp:7.3f}\n") - # vs model + fp.write("%3d\n" % nlayer) - for vel, depth in zip(vs, dep, strict=False): + for vel, depth in zip(velocity_s, depths, strict=True): fp.write(f"{vel:5.2f} {depth:7.2f} {vdamp:7.3f}\n")