Updated Volume Resampler

* Added a csv exporter for inversion outputs
* Fixed minor bugs in the volume resampler

seismic/pick_harvester/export_xyz.py

@@ -0,0 +1,93 @@
+#!/bin/env python
+"""
+Description:
+    Generate an xyz file from tomographic inversion output
+
+References:
+
+CreationDate:   01/06/23
+Developer:      rakib.hassan@ga.gov.au
+
+Revision History:
+    LastUpdate:     01/06/23   RH
+    LastUpdate:     dd/mm/yyyy  Who     Optional description
+"""
+
+import numpy as np
+import click
+import traceback
+
+
+def generate_xyz(input_file_name, output_file_name):
+    fh = open(input_file_name, 'r')
+
+    ncx = ncy = ncz = dx = dy = 0
+    xs = ys = zs = vals = None
+
+    #===============================================================
+    # Read data
+    #===============================================================
+    try:
+        ncx, ncy, ncz, dx, dy, romi, rami, r1, r2, r3 = map(float, fh.readline().split())
+        ncx, ncy, ncz = map(int, [ncx, ncy, ncz])
+
+        cz = np.array(list(map(float, fh.readline().split())))
+
+        cx = np.linspace(romi + dx/2, romi + ncx * dx - dx/2, ncx)
+        cy = np.linspace(rami + dy/2, rami + ncy * dy - dy/2, ncy)[::-1]
+        cx[cx > 180.] -= 360
+
+        vals = []
+        for line in fh.readlines():
+            vals.append(list(map(float, line.split())))
+        # end for
+        vals = np.array(vals).flatten()
+
+        assert len(vals) == ncx * ncy * ncz
+    except Exception as e:
+        print(traceback.format_exc())
+        raise RuntimeError('Failed to read file: {}..'.format(input_file_name))
+    # end try
+    fh.close()
+
+    #===============================================================
+    # Generate grids
+    #===============================================================
+
+    gz, gy, gx = np.meshgrid(cz, cy, cx)
+    gy = gy.transpose(1, 0, 2)
+    gz = gz.transpose(1, 0, 2)
+
+    assert np.all(gx[0, 0, :] == cx)
+    assert np.all(gy[0, :, 0] == cy)
+    assert np.all(gz[:, 0, 0] == cz)
+
+    gvals = np.zeros(gz.shape)
+
+    gvals[:, :, :] = vals.reshape((ncz, ncy, ncx))
+
+    np.savetxt(output_file_name,
+               np.array([gx.flatten(), gy.flatten(), gz.flatten(), gvals.flatten()]).T,
+               delimiter=',', header='lon, lat, depth(km), perturbation')
+# end func
+
+CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help'])
+@click.command(context_settings=CONTEXT_SETTINGS)
+@click.argument('input_file', type=click.Path(exists=True, dir_okay=False),
+                required=True)
+@click.argument('output_file', type=click.Path(exists=False, dir_okay=False),
+                required=True)
+def process(input_file, output_file):
+    """
+    INPUT_FILE: output of tomographic inversion
+    OUTPUT_FILE_NAME: name of output file, with columns: lon, lat, depth(km), perturbation
+    """
+
+    ofn = output_file if (output_file[-4:] == '.csv') else output_file + '.csv'
+
+    generate_xyz(input_file, ofn)
+# end func
+
+if __name__=="__main__":
+    process()
+# end if

utils/resample_volume.py

@@ -24,6 +24,9 @@
 import click
 import os.path as path
 from seismic.xcorqc.utils import rtp2xyz
+from pyproj import Proj
+
+DEFAULT_PROJ = 3577  # AU Albers
 
 # define utility functions
 def rtp2xyz(r, theta, phi):
@@ -54,7 +57,7 @@ def __init__(self, fname, lonf, latf, zf, nn=1, p=4):
         self._zf = zf
 
         # file format [lon, lat, z(depth km), perturbation]
-        data = np.loadtxt(self._fname)
+        data = np.loadtxt(self._fname, delimiter=',')
         self._gridLLZ = data[:, :3]
         self._values = data[:, -1]
 
@@ -84,15 +87,17 @@ def __init__(self, fname, lonf, latf, zf, nn=1, p=4):
         lato = self._gridLLZ[:self._nx*self._ny:self._nx, 1] # original lats
         zo = self._gridLLZ[::self._nx*self._ny, 2] # original zs
 
+        lono[lono < 0] += 360
         # preserve distribution of nodes, especially for z nodes,
         # which become sparser as we go deeper
         lonio = interp1d(np.linspace(0, len(lono), len(lono)), lono)
         latio = interp1d(np.linspace(0, len(lato), len(lato)), lato)
         zio = interp1d(np.linspace(0, len(zo), len(zo)), zo)
 
-        lonr = lonio(np.linspace(0, len(lono), len(lono) * self._lonf))
-        latr = latio(np.linspace(0, len(lato), len(lato) * self._latf))
-        zr = zio(np.linspace(0, len(zo), len(zo) * self._zf))
+        lonr = lonio(np.linspace(0, len(lono), int(len(lono) * self._lonf)))
+        latr = latio(np.linspace(0, len(lato), int(len(lato) * self._latf)))
+        zr = zio(np.linspace(0, len(zo), int(len(zo) * self._zf)))
+        lonr[lonr > 180] -= 360
 
         self._rgridLon, self._rgridLat, self._rgridZ = np.meshgrid(lonr, latr, zr, indexing='ij')
 
@@ -216,41 +221,54 @@ def WriteNetcdf(self, fnamestem):
     # end func
 
     def WriteSGrid(self, fnamestem):
-        def write_header(hfn, dfn, mname):
+        def write_header(hfn, dfn, pname):
             hlines = [  'GOCAD SGrid 1',
                         'HEADER {',
-                        'name:%s'%(mname),
+                        'name:%s'%(pname),
                         '*volume:true',
                         'ascii:on',
                         '*painted:on',
-                        '*painted*variable:p',
+                        '*painted*variable:{}'.format(pname),
                         'last_selected_folder:Volumes',
                         '*volume*grid:true',
                         '*volume*solid:true',
                         '*cube*solid:false',
                         '*cube*grid:false',
                         '*volume*points:false',
                         '}',
+                        'GOCAD_ORIGINAL_COORDINATE_SYSTEM',
+                        'NAME Default',
+                        'PROJECTION Unknown',
+                        'DATUM Unknown',
+                        'AXIS_NAME "X" "Y" "Z"',
+                        'AXIS_UNIT "m" "m" "m"',
+                        'ZPOSITIVE Depth',
+                        'END_ORIGINAL_COORDINATE_SYSTEM',
                         'AXIS_N %d %d %d' % (self._rnx, self._rny, self._rnz),
                         'PROP_ALIGNMENT POINTS',
                         'ASCII_DATA_FILE %s' % path.basename(dfn),
-                        'PROPERTY 1 "p"',
-                        'PROPERTY_CLASS 1 "p"',
-                        'END]' ]
+                        'PROPERTY 1 "{}"'.format(pname),
+                        'PROPERTY_CLASS 1 "{}"'.format(pname),
+                        'END' ]
 
             fo = open(hfn, 'w+')
             for l in hlines: fo.write(l+'\n')
             fo.close()
         # end func
 
         def write_data(dfn):
+            p = Proj(DEFAULT_PROJ)  # AU Albers
+
+
             lons, lats, ds = self._rgridLon.flatten(), self._rgridLat.flatten(), \
                              self.RADIUS - self._rgridZ.flatten()
 
             kxyz, jxyz, ixyz = np.meshgrid(np.arange(self._rnz), np.arange(self._rny),
                                            np.arange(self._rnx), indexing='ij')
 
-            np.savetxt(dfn, np.array([lons, lats, -ds*1e3, self._img,
+            px, py = p(lons, lats)
+
+            np.savetxt(dfn, np.array([px, py, ds*1e3, self._img,
                                       ixyz.flatten(), jxyz.flatten(),
                                       kxyz.flatten()]).T, fmt='%4.4f %4.4f %4.4f %5.7f %d %d %d')
         # end func
@@ -259,7 +277,7 @@ def write_data(dfn):
         dfn = fnamestem + '.ascii.txt'
         modelname = path.basename(fnamestem)
 
-        write_header(hfn, dfn, modelname)
+        write_header(hfn, dfn, 'p')
         write_data(dfn)
     # end func
 # end class
@@ -286,7 +304,7 @@ def process(input, output_file_stem, output_type, resampling_factors,
     """
     Script for super- or sub-sampling a 3D volume using inverse-distance-weighted interpolation.
 
-    INPUT: Path to xyz file, with nodal values on a regular 3D grid\n
+    INPUT: Path to four-columned (lon, lat, depth(km), perturbation) csv file\n
 
     Example Usage: python resample_volume.py P-wave_1x1.it3.xyz --resampling-factors 2x2x2 --output-file-stem /tmp/tomo
     """