GeoscienceAustralia / mtpy2

This Repository has been moved to https://github.com/MTgeophysics/mtpy
GNU General Public License v3.0
5 stars 1 forks source link

Convert model for visualisation software (ParaView, Gocad) #2

Open basaks opened 8 years ago

basaks commented 8 years ago

Jingming to provide clarification/code samples.

basaks commented 8 years ago

I have made some progress on this on git # 4be65f51e161035f273df3404ebbd0a324d3208c. Now I can create a .dat file as demonstrated in the main section of mtpy.modeling.modem_new.py. Edit: Since then modem_new.py was renamed to modem.py.

I use the .edi files from MTPy_development/test_data as the input edi_file_list.

    md = Data(edi_list=edi_file_list,
                period_min=.1,
                period_max=300,
                max_num_periods=12)

    md.write_data_file(
        save_path=MTPY_DEVELOPMENT_TEST_DATA,
        fn_basename='test.dat'
    )

The generated test.data is attached here (renamed with a .txt extension as github does not allow .dat upload): test.dat.txt

However, I need clarification on how to define abs_err for tip data.

Here is the relevant code block: lines 856-880 in mtpy.modeling.modem_new.

if comp.find('t') == 0:
    rel_err = self.error_tipper
    abs_err = 0.0  # TODO: fix abs_err for tip
elif comp.find('z') == 0:
    if self.error_type == 'floor':
        rel_err = self.data_array[ss][c_key+'_err'][ff, z_ii, z_jj]/\
                  abs(zz)
        if rel_err < self.error_floor/100.:
            rel_err = self.error_floor/100.
        abs_err = rel_err*abs(zz)
    elif self.error_type == 'value':
        abs_err = abs(zz)*self.error_value/100.

    elif self.error_type == 'egbert':
        d_zxy = self.data_array[ss]['z'][ff, 0, 1]
        d_zyx = self.data_array[ss]['z'][ff, 1, 0]
        abs_err = np.sqrt(abs(d_zxy*d_zyx))*\
                  self.error_egbert/100.
    elif self.error_type == 'floor_egbert':
        abs_err = self.data_array[ss][c_key+'_err'][ff, z_ii, z_jj]
        d_zxy = self.data_array[ss]['z'][ff, 0, 1]
        d_zyx = self.data_array[ss]['z'][ff, 1, 0]
        if abs_err < np.sqrt(abs(d_zxy*d_zyx))*self.error_egbert/100.:
            abs_err = np.sqrt(abs(d_zxy*d_zyx))*self.error_egbert/100.
    else:
        raise
if abs_err == 0.0:
    abs_err = 1e3  # TODO: clarify why
    print ('error at {0} is 0 for period {1}'.format(
            sta, per)+'set to 1e3')  # SB: why?

This line assigns an arbitrary value of 1e3 for abs_err=0, which seems very peculiar.

basaks commented 8 years ago

Now can write_vtk_station_file, a .vtu file extension. Also tested in revision # 1d7a848bae2dc8123897814ad1e2aa0547d0127b.

basaks commented 8 years ago

Need a ws3dinv data file to test convert_ws3dinv_data_file in Data class.

basaks commented 8 years ago

Please review modem_stations.vtu.zip

This was created using the modeling.modem.Data.write_vtk_station_file and the .edi files from MTPy_development/test_data as the input edi_file_list.

basaks commented 8 years ago

Please review ModEM_Model.ws.zip

This was created using the modeling.modem.Model.write_model_file and the .edi files from MTPy_development/test_data as the input edi_file_list.

Tests are available in mtpy/tests/modeling/test_modem.py.