Problem writing rfile - Githubissues

Hi Shahar - I'm having trouble with a very simple example of reading the berkeley.rfile and writing it out verbatim. I need to do some editing of rfiles and tried this case as a test exercise to verify reading and writing.

The script is attached. I have no problem with the topography block and if I set nblock=1 it can write the NEW rfile (fnew) correctly. However, when I try to write the volume blocks the block headers have the wrong parameters (hh, hv, z0, nc, ni, nj, nk). I have counted bytes and the resulting file has the correct size, but the data is clearly wrong.

It's probably a user error and not related to pySW4, but I'm wondering if you see anything wrong with my code. You may have come across such a case ...

Thanks! Artie

# usgs_rfile_edit.py
"""
This is an exercise for reading rfiles with pySW4 and writing out the same
This test case reads berkeley.rfile 
Later will need to edit material porperties
Notes: ordering of material properties is: rho, vp, vs, qp, qs
"""
import numpy as np
import pySW4

f = './berkeley.rfile'

#(magic, precision, attenuation, az, lon0, lat0, mlen, proj_str, nb) = pySW4.prep.rfileIO.read_hdr(f)
rfile = pySW4.prep.rfileIO.read(f, block_number='all', verbose=True)

# For testing, number of blocks for the output NEW rfile
# Set nblock = 1 to only write topography and lines 32-33
#nblocks = 1

# Open file handle the the NEW rfile
fnew = open('./test.rfile', 'wb')
# write header
pySW4.prep.rfileIO.write_hdr(fnew, 
                             rfile.magic, 
                             rfile.precision, 
                             rfile.attenuation, 
                             rfile.az, 
                             rfile.lon0, 
                             rfile.lat0, 
                             rfile.proj_str,
                             #nblocks)
                             rfile.nb)
rfile_hdr = [rfile.magic, rfile.precision, rfile.attenuation, rfile.az, rfile.lon0, rfile.lat0, rfile.mlen, rfile.proj_str, rfile.nb]

bytes_rfile_hdr = 0 
for item in rfile_hdr:
    bytes_rfile_hdr += item.nbytes
bytes_total = bytes_rfile_hdr
print('bytes_rfile_hdr: ', bytes_rfile_hdr)

# Loop over blocks
for ib in range(rfile.nb):
#for ib in range(nblocks):
    print()
    print('ib: ', ib)
    print('hh, hv, z0: ', rfile.blocks[ib].hh, rfile.blocks[ib].hv, rfile.blocks[ib].z0)
    print(' ni, nj, nk, nc:', rfile.blocks[ib].ni, rfile.blocks[ib].nj, rfile.blocks[ib].nk, rfile.blocks[ib].nc)

    # write_block_hdr(f, hh, hv, z0, nc, ni, nj, nk)
    pySW4.prep.rfileIO.write_block_hdr(fnew, 
                                       rfile.blocks[ib].hh, 
                                       rfile.blocks[ib].hv, 
                                       rfile.blocks[ib].z0, 
                                       rfile.blocks[ib].nc, 
                                       rfile.blocks[ib].ni, 
                                       rfile.blocks[ib].nj, 
                                       rfile.blocks[ib].nk)

    # Write block header without pySW4
    # fnew.write(np.float64(rfile.blocks[ib].hh))
    # fnew.write(np.float64(rfile.blocks[ib].hv))
    # fnew.write(np.float64(rfile.blocks[ib].z0))
    # fnew.write(np.int32(rfile.blocks[ib].nc))
    # fnew.write(np.int32(rfile.blocks[ib].ni))
    # fnew.write(np.int32(rfile.blocks[ib].nj))
    # fnew.write(np.int32(rfile.blocks[ib].nk))

    bytes_block_hdr = int(3 * 8 + 4 * 4)
    print('bytes_block_hdr: ', bytes_block_hdr)
    bytes_total += bytes_block_hdr

    # For topography ib == 0 and rfile.blocks[ib].nk == 1
    if ib == 0:
        print('write topography')
        ibytes = 0
        for i in range(rfile.blocks[ib].ni):
            for j in range(rfile.blocks[ib].nj):
                fnew.write(np.float32(rfile.blocks[ib].data[i,j]))
                ibytes += 4
        bytes_block = int (rfile.blocks[ib].ni * rfile.blocks[ib].nj * 4)
        bytes_total += bytes_block
        print('ibytes: ', ibytes)

    # For volumetric properties the block has nk depth levels and nc=5 material properties
    if ib > 0:
        print('write volume block')
        ibytes = 0
        for i in range(rfile.blocks[ib].ni):
            for j in range(rfile.blocks[ib].nj):
                for k in range(rfile.blocks[ib].nk):
                    for c in range(rfile.blocks[ib].nc):
                        fnew.write(np.float32(rfile.blocks[ib].data[i,j,k,c]))
                        ibytes += 4
        bytes_block = int (rfile.blocks[ib].ni * rfile.blocks[ib].nj * rfile.blocks[ib].nk * rfile.blocks[ib].nc * 4)
        bytes_total += bytes_block
        print('ibytes: ', ibytes)
    print(ib, bytes_block_hdr, bytes_block)
print(bytes_total)
fnew.close()

# Finally, read the NEW rfile, header and block headers should be identical
rfile1 = pySW4.prep.rfileIO.read('./test.rfile', block_number='all', verbose=True)

shaharkadmiel / pySW4

Problem writing rfile #19