@@ -79,7 +79,7 @@ def print_model_info(fout, cmd, clargs):
7979 print ('\t {:20s} -> {:}' .format (key , clargs [key ]), file = fout )
8080 return
8181
82- def compute__ (grd , sta_list_utm , utmzone , fout , fstats , vprint_fun , ** dargs ):
82+ def compute__ (igrd , sta_list_utm , utm_lcm , fout , fstats , vprint_fun , ** dargs ):
8383 """ Function to perform the bulk of a Strain Tensor estimation.
8484 For each of the grid cells, a ShenStrain object will be created, using
8585 the list of stations and the **dargs options.
@@ -111,23 +111,27 @@ def compute__(grd, sta_list_utm, utmzone, fout, fstats, vprint_fun, **dargs):
111111 not flushed before returning or something). Anyway, always close the
112112 streams before exiting.
113113 """
114+ #print('--> Thread given grid : X:{:}/{:}/{:} Y:{:}/{:}/{:}'.format(igrd.x_min, igrd.x_max, igrd.x_step, igrd.y_min, igrd.y_max, igrd.y_step))
114115 node_nr , nodes_estim = 0 , 0
115- for x , y in grd :
116+ for x , y in igrd :
116117 clat , clon = radians (y ), radians (x )
117- N , E , ZN , _ = ell2utm (clat , clon , Ellipsoid ("wgs84" ), utmzone )
118- assert ZN == utmzone
118+ #print('--> computing tensor at lon {:}, lat {:}'.format(x, y))
119+ N , E , ZN , lcm = ell2utm (clat , clon , Ellipsoid ("wgs84" ), utm_lcm )
120+ #assert ZN == utmzone
121+ assert utm_lcm == lcm
119122 vprint_fun ('[DEBUG] Grid point at {:+8.4f}, {:8.4f} or E={:}, N={:}' .format (
120123 x , y , E , N ))
121124 if not dargs ['multiproc_mode' ]:
122125 print ('[DEBUG] {:5d}/{:7d}' .format (node_nr + 1 , grd .xpts * grd .ypts ), end = "\r " )
123126 ## Construct the Strain instance, with all args (from input)
124- sstr = ShenStrain (E , N , sta_list_utm , ** dargs )
127+ # sstr = ShenStrain(E, N, sta_list_utm, **dargs)
128+ sstr = ShenStrain (E , N , clat < 0e0 , sta_list_utm , ** dargs )
125129 ## check azimouth coverage (aka max β angle)
126130 if degrees (max (sstr .beta_angles ())) <= dargs ['max_beta_angle' ]:
127131 try :
128132 sstr .estimate ()
129133 vprint_fun ('[DEBUG] Computed tensor at {:+8.4f} {:+8.4f} for node {:3d}/{:3d}' .format (x , y , node_nr + 1 , grd .xpts * grd .ypts ))
130- sstr .print_details_v2 (fout , utmzone )
134+ sstr .print_details_v2 (fout , utm_lcm )
131135 if fstats : print ('{:+9.4f} {:+10.4f} {:6d} {:14.2f} {:10.2f} {:12.3f}' .format (x ,y ,len (sstr .__stalst__ ), sstr .__options__ ['d_coef' ],sstr .__options__ ['cutoff_dis' ], sstr .__sigma0__ ), file = fstats )
132136 nodes_estim += 1
133137 except RuntimeError :
@@ -156,8 +160,8 @@ class myFormatter(
156160 Dionysos Satellite Observatory\n
157161Send bug reports to:
158162 Xanthos Papanikolaou, xanthos@mail.ntua.gr
159- Dimitris Anastasiou,danast@mail.ntua.gr
160- November, 2017 '''
163+ Dimitris Anastasiou,dganastasiou@gmail.com
164+ September, 2021 '''
161165
162166parser .add_argument ('-i' , '--input-file' ,
163167 default = argparse .SUPPRESS ,
@@ -405,15 +409,16 @@ class myFormatter(
405409 ## TODO is this mean_lon the optimal?? or should it be the region's mean longtitude
406410 ##
407411 mean_lon = degrees (sum ([ x .lon for x in sta_list_ell ]) / len (sta_list_ell ))
408- utm_zone = floor (mean_lon / 6 )+ 31
409- utm_zone = utm_zone + int (utm_zone <= 0 )* 60 - int (utm_zone > 60 )* 60
410- vprint ('[DEBUG] Mean longtitude is {} deg.; using Zone = {} for UTM' .format (mean_lon , utm_zone ))
412+ #utm_zone = floor(mean_lon/6)+31
413+ #utm_zone = utm_zone + int(utm_zone<=0)*60 - int(utm_zone>60)*60
414+ lcm = radians (floor (mean_lon ))
415+ #print('[DEBUG] Mean longtitude is {} deg.; using Zone = {} for UTM'.format(mean_lon, utm_zone))
411416 sta_list_utm = deepcopy (sta_list_ell )
412417 for idx , sta in enumerate (sta_list_utm ):
413- N , E , Zone , lcm = ell2utm (sta .lat , sta .lon , Ellipsoid ("wgs84" ), utm_zone )
418+ N , E , Zone , lcm = ell2utm (sta .lat , sta .lon , Ellipsoid ("wgs84" ), lcm )
414419 sta_list_utm [idx ].lon = E
415420 sta_list_utm [idx ].lat = N
416- assert Zone == utm_zone , "[ERROR] Invalid UTM Zone."
421+ # assert Zone == utm_zone, "[ERROR] Invalid UTM Zone."
417422 vprint ('[DEBUG] Station list transformed to UTM.' )
418423
419424 ## Open file to write Strain Tensor estimates; write the header
@@ -426,12 +431,12 @@ class myFormatter(
426431 if args .one_tensor :
427432 print ('[DEBUG] Estimating Strain Tensor at region\' s barycentre.' )
428433 if args .method == 'shen' :
429- sstr = ShenStrain (0e0 , 0e0 , sta_list_utm , ** dargs )
434+ sstr = ShenStrain (0e0 , 0e0 , False , sta_list_utm , ** dargs )
430435 else :
431- sstr = ShenStrain (0e0 , 0e0 , sta_list_utm , weighting_function = 'equal_weights' )
436+ sstr = ShenStrain (0e0 , 0e0 , False , sta_list_utm , weighting_function = 'equal_weights' )
432437 sstr .set_to_barycenter ()
433438 sstr .estimate ()
434- sstr .print_details (fout , utm_zone )
439+ sstr .print_details (fout , utm_lcm )
435440 fout .close ()
436441 write_station_info (sta_list_ell )
437442 print ('[DEBUG] Total running time: {:10.2f} sec.' .format ((time .time () - start_time )))
@@ -458,6 +463,7 @@ class myFormatter(
458463 ##+ coordinates in lon/lat pairs, in degrees!
459464 if args .multiproc_mode :
460465 grd1 , grd2 , grd3 , grd4 = grd .split2four ()
466+ print ('--> grid split to four!' )
461467 fout1 = open (".out.thread1" , "w" )
462468 fout2 = open (".out.thread2" , "w" )
463469 fout3 = open (".out.thread3" , "w" )
@@ -470,10 +476,14 @@ class myFormatter(
470476 else :
471477 fstats1 = fstats2 = fstats3 = fstats4 = None
472478 print ('[DEBUG] Estimating strain tensors in multi-threading mode' )
473- p1 = multiprocessing .Process (target = compute__ , args = (grd1 , sta_list_utm , utm_zone , fout1 , fstats1 , vprint ), kwargs = dargs )
474- p2 = multiprocessing .Process (target = compute__ , args = (grd2 , sta_list_utm , utm_zone , fout2 , fstats2 , vprint ), kwargs = dargs )
475- p3 = multiprocessing .Process (target = compute__ , args = (grd3 , sta_list_utm , utm_zone , fout3 , fstats3 , vprint ), kwargs = dargs )
476- p4 = multiprocessing .Process (target = compute__ , args = (grd4 , sta_list_utm , utm_zone , fout4 , fstats4 , vprint ), kwargs = dargs )
479+ #print('--> Thread will be given grid : X:{:}/{:}/{:} Y:{:}/{:}/{:}'.format(grd1.x_min, grd1.x_max, grd1.x_step, grd1.y_min, grd1.y_max, grd1.y_step))
480+ #print('--> Thread will be given grid : X:{:}/{:}/{:} Y:{:}/{:}/{:}'.format(grd2.x_min, grd2.x_max, grd2.x_step, grd2.y_min, grd2.y_max, grd2.y_step))
481+ #print('--> Thread will be given grid : X:{:}/{:}/{:} Y:{:}/{:}/{:}'.format(grd3.x_min, grd3.x_max, grd3.x_step, grd3.y_min, grd3.y_max, grd3.y_step))
482+ #print('--> Thread will be given grid : X:{:}/{:}/{:} Y:{:}/{:}/{:}'.format(grd4.x_min, grd4.x_max, grd4.x_step, grd4.y_min, grd4.y_max, grd4.y_step))
483+ p1 = multiprocessing .Process (target = compute__ , args = (grd1 , sta_list_utm , lcm , fout1 , fstats1 , vprint ), kwargs = dargs )
484+ p2 = multiprocessing .Process (target = compute__ , args = (grd2 , sta_list_utm , lcm , fout2 , fstats2 , vprint ), kwargs = dargs )
485+ p3 = multiprocessing .Process (target = compute__ , args = (grd3 , sta_list_utm , lcm , fout3 , fstats3 , vprint ), kwargs = dargs )
486+ p4 = multiprocessing .Process (target = compute__ , args = (grd4 , sta_list_utm , lcm , fout4 , fstats4 , vprint ), kwargs = dargs )
477487 [ p .start () for p in [p1 , p2 , p3 , p4 ]]
478488 [ p .join () for p in [p1 , p2 , p3 , p4 ]]
479489 for fl in [fout1 , fout2 , fout3 , fout4 ]:
@@ -498,7 +508,7 @@ class myFormatter(
498508 os .remove (".sta.thread" + str (fnr ))
499509
500510 else :
501- compute__ (grd , sta_list_utm , utm_zone , fout , fstats , vprint , ** dargs )
511+ compute__ (grd , sta_list_utm , lcm , fout , fstats , vprint , ** dargs )
502512 else :
503513 ## Using veis method. Compute delaunay triangles and estimate one tensor
504514 ##+ per triangle centre
@@ -515,9 +525,9 @@ class myFormatter(
515525 cy = (sta_list_utm [trng [0 ]].lat + sta_list_utm [trng [1 ]].lat + sta_list_utm [trng [2 ]].lat )/ 3e0
516526 ## Construct a strain instance, at the triangle's barycentre, with only
517527 ##+ 3 points (in UTM) and equal_weights weighting scheme.
518- sstr = ShenStrain (cx , cy , [sta_list_utm [trng [0 ]], sta_list_utm [trng [1 ]], sta_list_utm [trng [2 ]]], weighting_function = 'equal_weights' )
528+ sstr = ShenStrain (cx , cy , cy < 0e0 , [sta_list_utm [trng [0 ]], sta_list_utm [trng [1 ]], sta_list_utm [trng [2 ]]], weighting_function = 'equal_weights' )
519529 sstr .estimate ()
520- sstr .print_details (fout , utm_zone )
530+ sstr .print_details (fout , lcm )
521531 ## Print the triangle in the corresponding file (ellipsoidal crd, degrees)
522532 print ('> {:}, {:}, {:}' .format (sta_list_utm [trng [0 ]].name , sta_list_utm [trng [1 ]].name , sta_list_utm [trng [2 ]].name ), file = dlnout )
523533 print ('{:+8.5f} {:8.5f}\n {:+8.5f} {:8.5f}\n {:+8.5f} {:8.5f}\n {:+8.5f} {:8.5f}' .format (* [ degrees (x ) for x in [sta_list_ell [trng [0 ]].lon , sta_list_ell [trng [0 ]].lat , sta_list_ell [trng [1 ]].lon , sta_list_ell [trng [1 ]].lat , sta_list_ell [trng [2 ]].lon , sta_list_ell [trng [2 ]].lat , sta_list_ell [trng [0 ]].lon , sta_list_ell [trng [0 ]].lat ]]), file = dlnout )
0 commit comments