Added a script to export FWI models to csv

seismic/ASDFdatabase/FederatedASDFViewer.py

@@ -32,7 +32,7 @@
 from scipy.stats import circmean as cmean
 from collections import defaultdict
 from shapely import geometry
-import traceback
+from seismic.misc import print_exception
 
 class CustomPPSD(PPSD):
     def __init__(self, stats, skip_on_gaps=False,
@@ -134,13 +134,6 @@ def _plot_histogram(self, fig, draw=False, filename=None):
         return fig    
 # end class
 
-def printException(e: Exception):
-    exc_type, exc_obj, exc_tb = sys.exc_info()
-    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
-    print(exc_type, e, fname, exc_tb.tb_lineno)
-    print(traceback.format_exc())
-# end func
-
 class FigureImage(gui.Image):
     def __init__(self, **kwargs):
         super(FigureImage, self).__init__("/%s/get_image_data?update_index=0" % id(self), **kwargs)
@@ -340,7 +333,7 @@ def setMapImage(nc, minLon=None, minLat=None, maxLon=None, maxLat=None):
                     children['latBoundsBox'].children['max'].set_value(maxLat)
 
             except Exception as e:
-                printException(e)
+                print_exception(e)
             # end try
         # end func
 
@@ -552,7 +545,7 @@ def setTraceImage(nc, sc, lc, cc, st=None, et=None):
                 ti = FigureImage(fig=fig)
                 self.rowContainer.children[key].children['rightContainer'].children['plot'] = ti
             except Exception as e:
-                printException(e)
+                print_exception(e)
             # end try
         # end func
 
@@ -612,7 +605,7 @@ def startStepChanged(emitter, value=None):
                                      args=(nc, sc, lc, cc, st, et))
                 t.start()
             except Exception as e:
-                printException(e)
+                print_exception(e)
             # end try
         # end func
 
@@ -627,7 +620,7 @@ def removeWidget(emitter):
                                                    self.rowWidgetCount * ROW_WIDGET_HEIGHT * PADDING_FACTOR +
                                                    MAP_WIDGET_PADDING)
                 except Exception as e:
-                    printException(e)
+                    print_exception(e)
                 # end try
             # end if
         # end func

seismic/ASDFdatabase/xdmf2csv.py

@@ -0,0 +1,100 @@
+"""
+Description:
+    Exports contents of XDMF to CSV.
+
+References:
+
+CreationDate:   02/05/24
+Developer:      [email protected]
+
+Revision History:
+    LastUpdate:     02/05/24   RH
+"""
+
+import click
+import vtk
+from vtk.numpy_interface import dataset_adapter as dsa
+import numpy as np
+from seismic.misc import print_exception
+import pyproj
+import pandas as pd
+import os
+
+def tranforms_coords(x:np.ndarray, y:np.ndarray, z:np.ndarray) -> \
+    (np.ndarray, np.ndarray, np.ndarray):
+
+    transformer = pyproj.Transformer.from_crs(
+        {"proj": 'geocent', "ellps": 'WGS84', "datum": 'WGS84'},
+        {"proj": 'latlong', "ellps": 'WGS84', "datum": 'WGS84'})
+
+    xyz2lonlatalt = lambda x, y, z: np.vstack(transformer.transform(x, y, z,
+                                                                    radians=False)).T
+
+    result = xyz2lonlatalt(x, y, z)
+    result[:, -1] *= -1 # convert to depth
+
+    return result[:, 0], result[:, 1], result[:, 2] # lons, lats, depths
+# end func
+
+CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help'])
+@click.command(context_settings=CONTEXT_SETTINGS)
+@click.argument('input-file', required=True,
+                type=click.Path(exists=True))
+@click.argument('output-folder', required=True,
+                type=click.Path(exists=True))
+def process(input_file, output_folder):
+    reader = vtk.vtkXdmfReader()
+
+    try:
+        print('Reading input file: {}'.format(input_file))
+        reader.SetFileName(input_file)
+        reader.UpdateInformation()
+    except Exception as e:
+        print_exception(e)
+    # end func
+
+
+    pd2cd = vtk.vtkPointDataToCellData()
+    pd2cd.SetInputConnection(reader.GetOutputPort())
+    pd2cd.PassPointDataOn()
+
+    pd2cd.Update()
+
+    data = dsa.WrapDataObject(pd2cd.GetOutput())
+
+    data.PointData.keys()
+
+    xyz = data.GetPoints()
+
+    rkeys = ['VP', 'VSH', 'VSV']
+    valsDict = {}
+
+    for key in rkeys:
+        valsDict[key] = data.PointData[key]
+    # end for
+
+    # tranform coordinates from geocentric xyz to geographic (both in wgs84)
+    print('Transforming coordinates..')
+    lons, lats, depths = tranforms_coords(xyz[:, 0], xyz[:, 1], xyz[:, 2])
+
+    # write csv
+    ofn = os.path.splitext(os.path.basename(input_file))[0] + '.csv'
+    ofn = os.path.join(output_folder, ofn)
+
+    df = pd.DataFrame()
+
+    df['lons'] = lons
+    df['lats'] = lats
+    df['depths_km'] = depths
+    for key in rkeys:
+        df[key] = valsDict[key]
+    # end for
+
+    print('Writing output file: {}'.format(ofn))
+    df.to_csv(ofn, header=True, index=False)
+# end func
+
+
+if __name__=="__main__":
+    process()
+# end if

seismic/misc.py

@@ -16,6 +16,7 @@
 import os, glob, fnmatch, sys
 import numpy as np
 import logging
+import traceback
 logging.basicConfig()
 
 def setup_logger(name, log_file=None, level=logging.INFO, propagate=False):
@@ -83,4 +84,11 @@ def rtp2xyz(r, theta, phi):
     xout[:, 1] = rst * np.sin(phi)
     xout[:, 2] = r * np.cos(theta)
     return xout
+# end func
+
+def print_exception(e: Exception):
+    exc_type, exc_obj, exc_tb = sys.exc_info()
+    fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
+    print(exc_type, e, fname, exc_tb.tb_lineno)
+    print(traceback.format_exc())
 # end func