Modem updates

[kujaku11]

Update to functionality of Modem tools and updated documentation.

Description

• Added functions to Data
  • add/remove a station
  • remove component
  • add error to components
  • swap phase
  • remove static shift
  • write VTK file in right-hand or left-hand coordinate system
  • Added ocean_bottom argument in project_stations_on_topography to place stations at bottom of ocean
• Added functions to Model
  • write files for Geosoft and Leapfrog
  • Added arguments to add_topography_to_model2 to shift east and north incase the projection is off (often when UTM grids are crossed)
• All station locations taken care of with Station
• added function Station.write_shp_file to write a shapefile of station locations
• Fixed bug in Station when stations cross UTM grids, will now project all stations to the center point UTM zone, can specify that zone or its calculated.
• Updated import of evtk to pyevtk
• Updated documentation
• Updated tests
• Added functionality to gui/modem_plot_response_qt5.py to include the added functionality in Data
• Updated tests to use Path

Motivation and Context

Adding functionality

How Has This Been Tested?

Added test for Data.

Screenshots (if appropriate):

Types of changes

• [ ] Bug fix (non-breaking change which fixes an issue)
• [ x] New feature (non-breaking change which adds functionality)
• [ ] Breaking change (fix or feature that would cause existing functionality to change)

Checklist:

• [x ] My code follows the code style of this project.
• [ x] My change requires a change to the documentation.
• [x ] I have updated the documentation accordingly.
• [x ] I have added tests to cover my changes.
• [x ] All new and existing tests passed
• Well, Almost the rotation of mesh test fails, its only off by a few meters. So either recalculate or update default files.

[alkirkby]

ModEM RMS maps seems to be broken on this branch now? I tried running the example /examples/scripts/ModEM_PlotRMS.py and it doesn't work on modem_updates but works on develop:

Traceback (most recent call last):

File "C:\mtpywin\mtpy\examples\scripts\ModEM_PlotRMS.py", line 54, in plot_elements='both',

File "C:\mtpywin\mtpy\mtpy\modeling\modem\plot_rms_maps.py", line 240, in init self.plot()

File "C:\mtpywin\mtpy\mtpy\modeling\modem\plot_rms_maps.py", line 427, in plot rms, filt = self._calculate_rms(p_dict)

File "C:\mtpywin\mtpy\mtpy\modeling\modem\plot_rms_maps.py", line 259, in _calculate_rms :, self.period_index, ii, jj

IndexError: only integers, slices (:), ellipsis (...), numpy.newaxis (None) and integer or boolean arrays are valid indices

[alkirkby]

I also get an error when trying to create a datafile: running example mtpy\examples\scripts\ModEM_build_inputfiles.py:

ERROR:mtpy.modeling.modem.data.Data:'utf-32-le' codec can't decode bytes in position 0-3: code point not in range(0x110000) Traceback (most recent call last):

File "C:\mtpywin\mtpy\mtpy\modeling\modem\data.py", line 895, in fill_data_array data_array = self.get_relative_station_locations(mt_dict, data_array)

File "C:\mtpywin\mtpy\mtpy\modeling\modem\data.py", line 490, in get_relative_station_locations stations_obj.get_station_locations(mt_list)

File "C:\mtpywin\mtpy\mtpy\modeling\modem\station.py", line 267, in get_station_locations self.calculate_rel_locations()

File "C:\mtpywin\mtpy\mtpy\modeling\modem\station.py", line 284, in calculate_rel_locations self.station_locations["rel_east"] = self.east - self.center_point.east[0]

File "C:\mtpywin\mtpy\mtpy\modeling\modem\station.py", line 359, in center_point self.logger.info(f"Projecting lat, lon to UTM zone {center_location['zone'][0]}")

UnicodeDecodeError: 'utf-32-le' codec can't decode bytes in position 0-3: code point not in range(0x110000)

During handling of the above exception, another exception occurred:

Traceback (most recent call last):

File "C:\mtpywin\mtpy\examples\scripts\ModEM_build_inputfiles.py", line 78, in do.write_data_file()

File "C:\mtpywin\mtpy\mtpy\modeling\modem\data.py", line 1255, in write_data_file longitude_format=longitude_format,

File "C:\mtpywin\mtpy\mtpy\modeling\modem\data.py", line 906, in fill_data_array raise ValueError(error)

ValueError: 'utf-32-le' codec can't decode bytes in position 0-3: code point not in range(0x110000)

[kujaku11]

@alkirkby RMS Plots fails because there is an option to set the attribute period_index to "all".

From the logic it looks _calculate_rms only wants an integer to return the rms array. So the question is what is meant by "all"? Is it to plot each period separately, or is it supposed to be a single plot where all periods are averaged or something else?

If its plot a summarization of all periods in one plot, I added a fix using np.nanmean() and it works.

[kujaku11]

Could not replicate your errors in ModEM_build_inputfiles, but I did update that example to import the directories for EDI files and where the output files go. Also updated to plot the topography.

[alkirkby]

@alkirkby RMS Plots fails because there is an option to set the attribute period_index to "all".

From the logic it looks _calculate_rms only wants an integer to return the rms array. So the question is what is meant by "all"? Is it to plot each period separately, or is it supposed to be a single plot where all periods are averaged or something else?

If its plot a summarization of all periods in one plot, I added a fix using np.nanmean() and it works.

It's not a straight mean it's an RMS - I did some checking to make sure the overall RMS matches what is in the logfile in ModEM and it needs to be calculated as an RMS so i do the same for the by station RMS. There should be a bunch of properties in the rms_array in teh residual object. THe one you need for component is rms_array['rms_z_component']. Check out _calculate_rms, line 240 of rms_maps in develop

[kujaku11]

It's not a straight mean it's an RMS - I did some checking to make sure the overall RMS matches what is in the logfile in ModEM and it needs to be calculated as an RMS so i do the same for the by station RMS. There should be a bunch of properties in the rms_array in teh residual object. THe one you need for component is rms_array['rms_z_component']. Check out _calculate_rms, line 240 of rms_maps in develop

Sounds right, I updated accordingly.

[alkirkby]

I am still having trouble getting modem_build_inputfiles.py to work. Same error as before. It seems to have a problem with center_location['zone']. Will keep looking into it

[alkirkby]

@kujaku11 when you run the example code, what do you get for center_location['zone']? for me it is empty (None) and I think that's where the error is coming in but I'm not completely sure. It's frustrating as sometimes it gives an error when I try and print center_location['zone'] so really hard to debug

[kujaku11]

@kujaku11 when you run the example code, what do you get for center_location['zone']? for me it is empty (None) and I think that's where the error is coming in but I'm not completely sure. It's frustrating as sometimes it gives an error when I try and print center_location['zone'] so really hard to debug

@alkirkby That is odd. When I run mtpy\examples\scripts\ModEM_build_inputfiles.py I get the following:

ModEM Data Object:
    Number of stations: 28
    Number of periods:  26
    Period range:  
        Min: 0.0017782794100389228 s
        Max: 3162.2776601683795 s
    Rotation angle:     0.0
    Data center:        
         latitude:  -20.5100 deg
         longitude: 138.0100 deg
         Elevation: -516.0 m
         Easting:   188142.4792 m
         Northing:  7729227.6114 m
         UTM zone:  54K
    Model EPSG:         28354
    Model UTM zone:     None
    Impedance data:     True
    Tipper data:        True

And get the following plots:

Have you pulled the current version?

[alkirkby]

Yep I pull the latest version each day, if you print center_location['zone'] what do you get? is it None?

[kujaku11]

Yep I pull the latest version each day, if you print center_location['zone'] what do you get? is it None?

do.center_point
Out[9]: 
rec.array([ (-20.51,  138.01,  188142.47916241,  7729227.61144812, -516., '54K')], 
          dtype=[('lat', '<f8'), ('lon', '<f8'), ('east', '<f8'), ('north', '<f8'), ('elev', '<f8'), ('zone', '<U4')])

do.station_locations.center_point
Out[10]: 
rec.array([ (-20.51,  138.01,  188142.47916241,  7729227.61144812, -516., '54K')], 
          dtype=[('lat', '<f8'), ('lon', '<f8'), ('east', '<f8'), ('north', '<f8'), ('elev', '<f8'), ('zone', '<U4')])

[kujaku11]

And the error you are getting is:

ValueError: 'utf-32-le' codec can't decode bytes in position 0-3: code point not in range(0x110000)

That sounds like an encoding issue. What type of machine are you on?

[kujaku11]

What do you get if you print do.station_locations.utm_zone?

[alkirkby]

the do.station_locations object is None. Perhaps we should catch up on skype to try and sort this out?

[kujaku11]

I've finished most of the updates to the ModEM module. Merged with the most recent version of develop.

[kujaku11]

There was a merge conflict that I messed up in utils/calculator.py with rounding numbers. Updated to what is on the develop branch.

[alkirkby]

just running the example ModEM_plot_slices_on_basemap.py and it produces nothing in the plot, e.g.:

Compare with output from develop branch:

I will continue testing and look into it at the end.

[alkirkby]

examples/scripts/ModEM_PlotPTmap.py needed a small tweak (on this branch and develop) so that files save correctly. I have made the change on the develop branch.

[alkirkby]

Everything else works fine for me so once basemap plot is resolved should be good to go. Sorry I introduced conflicts - I have updated occam2d to fix issue #135 so I think that just needs to be brought across but check you are happy with it.

[kujaku11]

Everything else works fine for me so once basemap plot is resolved should be good to go. Sorry I introduced conflicts - I have updated occam2d to fix issue #135 so I think that just needs to be brought across but check you are happy with it.

Any resolution on this, I'm having issues installing mpl_toolkits.basemap and haven't been able to test it.

[alkirkby]

Ah sorry no I haven't looked into it. You need to install basemap e.g. conda install basemap (I think) then you should be able to get the modules from mpl_toolkits. I will look at it now

[alkirkby]

I found the bug but need to go now unfortunately so haven't managed to fix it. It's plotting in the wrong spot because the center point is projected from lat/long to east/north using the utm zone of that center point, not the model_epsg that is optionally passed into the data object. It should just be a minor tweak (hopefully)