Add shift to discrete topography for NSEM

[dccowan]

Summary

The PR adds a utility for shifting locations relative to discrete surface topography.

Squash and Merge Summary

Add new shift_to_discrete_topography function that shifts locations relative to discrete surface topography. Add also a new get_discrete_topography function that generates discrete topography locations out of a mesh and its active cells. We measure electric fields at the Earth's surface when performing MT surveys. Like DC/IP, the original measurement locations of the electric fields can end up in air cells when we define discrete surface topography. These functions allow the user to shift locations relative to discrete topography on Tensor and Tree meshes. For Airborne NSEM, they also allow the user to preserve the original flight heights. Deprecate the following functions: gettopoCC, drapeTopotoLoc, and closestPointsGrid.

PR Checklist

• If this is a work in progress PR, set as a Draft PR
• Linted my code according to the style guides.
• Added tests to verify changes to the code.
• Added necessary documentation to any new functions/classes following the
  expect style.
• Marked as ready for review (if this is was a draft PR), and converted
  to a Pull Request
• Tagged @simpeg/simpeg-developers when ready for review.

What does this implement/fix?

We measure electric fields at the Earth's surface when performing MT surveys. Like DC/IP, the original measurement locations of the electric fields can end up in air cells when we define discrete surface topography. This utility allows the user to shift locations relative to discrete topography on Tensor and Tree meshes. For Airborne NSEM, this utility also allows the user to preserve the original flight heights.

Additional information

Unlike the DC/IP utility that alters a survey object, this utility is inputs the original locations and outputs the shifted ones. I could see us needing similar functionality for NSEM survey objects if they are created directly by loading some type of data file (like DC/IP).

[codecov]

Codecov Report

❌ Patch coverage is 88.77551% with 22 lines in your changes missing coverage. Please review.
✅ Project coverage is 81.45%. Comparing base (1845349) to head (317c3de).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
simpeg/utils/mesh_utils.py	80.30%	9 Missing and 4 partials ⚠️
...mpeg/electromagnetics/static/utils/static_utils.py	70.00%	6 Missing ⚠️
tests/utils/test_mesh_utils.py	96.92%	1 Missing and 1 partial ⚠️
...mpeg/electromagnetics/static/resistivity/survey.py	85.71%	0 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1683      +/-   ##
==========================================
- Coverage   81.47%   81.45%   -0.02%     
==========================================
  Files         418      419       +1     
  Lines       54763    54932     +169     
  Branches     5200     5224      +24     
==========================================
+ Hits        44617    44746     +129     
- Misses       8753     8786      +33     
- Partials     1393     1400       +7     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[jcapriot]

I feel like this should just be moved to a utility function outside of any of the static or natural_source modules. Having the same function defined in many locations goes against the "zen" of python.

[dccowan]

I feel like this should just be moved to a utility function outside of any of the static or natural_source modules. Having the same function defined in many locations goes against the "zen" of python.

agreed

[dccowan]

I noticed something important about these projections functions. When originally developed for DC/IP, the projection ensured the final electrode locations didn't end up on an 'exposed' edge or corner of a mesh cell. This was accomplished by shifting the horizontal position of the original locations to lie over the nearest cell centers. This is somewhat necessary for electric field measurements but not necessarily for magnetic field measurements.

Does anyone have any opinions on this?

[jcapriot]

I noticed something important about these projections functions. When originally developed for DC/IP, the projection ensured the final electrode locations didn't end up on an 'exposed' edge or corner of a mesh cell. This was accomplished by shifting the horizontal position of the original locations to lie over the nearest cell centers. This is somewhat necessary for electric field measurements but not necessarily for magnetic field measurements.

Does anyone have any opinions on this?

It’s definitely useful for the CC dc simulation, as shifting it to the centers is consistent with the 0 Neumann BC at the top. For the nodal formulation, this is less of an issue as the nodes (which are where the interpolated values live) lay directly on the boundary faces.

[dccowan]

I noticed something important about these projections functions. When originally developed for DC/IP, the projection ensured the final electrode locations didn't end up on an 'exposed' edge or corner of a mesh cell. This was accomplished by shifting the horizontal position of the original locations to lie over the nearest cell centers. This is somewhat necessary for electric field measurements but not necessarily for magnetic field measurements.
Does anyone have any opinions on this?

It’s definitely useful for the CC dc simulation, as shifting it to the centers is consistent with the 0 Neumann BC at the top. For the nodal formulation, this is less of an issue as the nodes (which are where the interpolated values live) lay directly on the boundary faces.

In the original implementation, forcing shifted locations over cell centers only happened for tree meshes. I added a "shift_horizontal" option to the new utilities. By default it is set to True since we mostly use this for electric dipole receivers.

[dccowan]

I'm not sure I like the kwarg 'option' for setting discrete surface topography base on the top or center of cells. In the original utilities it was fine because the number of input arguments was small. We may want to choose something more descriptive? Although it might mean I need to change the deprecation of the original functionality. Thoughts?

[santisoler]

Thanks for putting this together, @dccowan. I left a few comments below. Nothing too big. Let me know what do you think.

[santisoler]

By reading the code I noticed that for a TensorMesh, the function returns a mesh and an array. While for the TreeMesh case, it returns two arrays.

Would it be possible to stay with only one type of output? Also, we should update the docstrings so the returns match the output of the function.

[dccowan]

For the TensorMesh, it looks like we are leveraging the closest_points_index method. But the TreeMesh class also appears to have the method closest_points_index. I'll need to think a little bit more about what we want the outputs to be for functionality that users will interact with.

[santisoler]

I'm not sure I like the kwarg 'option' for setting discrete surface topography base on the top or center of cells. In the original utilities it was fine because the number of input arguments was small. We may want to choose something more descriptive? Although it might mean I need to change the deprecation of the original functionality. Thoughts?

What about naming it to something like shift_vertically_on? So a function call would read as:

shift_to_discrete_topography(
        mesh,
        pts,
        active_cells,
        shift_vertically_on="top",
        shift_horizontal=shift_horizontal,
        heights=heights,
    )

If we are going for a change in the signature of the function (and not just renaming it), I think we could avoid using the deprecate_function. Instead we can keep the old function (not delete it) and decorate it with @deprecated. See this as an example of how it's being used in SimPEG:

simpeg/simpeg/potential_fields/gravity/survey.py

Lines 109 to 133 in cba34e9

	@property
	@deprecated(
	"The `components` property is deprecated, "
	"and will be removed in SimPEG v0.25.0. "
	"Within a gravity survey, receivers can contain different components. "
	"Iterate over the sources and receivers in the survey to get "
	"information about their components.",
	category=FutureWarning,
	)
	def components(self):
	"""Number of components measured at each receiver.
	.. deprecated:: 0.24.0
	The `components` property is deprecated, and will be removed in
	SimPEG v0.25.0. Within a gravity survey, receivers can contain
	different components. Iterate over the sources and receivers in the
	survey to get information about their components.
	Returns
	-------
	int
	Number of components measured at each receiver.
	"""
	return self.source_field.receiver_list[0].components

Just note that we import it through warnings or typing_extensions depending on the Python version:

simpeg/simpeg/potential_fields/gravity/survey.py

Lines 5 to 11 in cba34e9

	try:
	from warnings import deprecated
	except ImportError:
	# Use the deprecated decorator provided by typing_extensions (which
	# supports older versions of Python) if it cannot be imported from
	# warnings.
	from typing_extensions import deprecated

[santisoler]

@dccowan, I just pushed a few changes. Mostly related to the deprecation warnings, admonitions and tests. We remove deprecated bits of code after two minor releases since the deprecation is released. So, since these functions will be deprecated in v0.25.0 (upcoming release), we should remove them in v0.27.0.

I also restored the old ind_active argument in the deprecated gettopoCC function. Renaming arguments introduces a breaking change. For example, if any user has a code calling that function like: gettopoCC(mesh=mesh, ind_active=active_cells), it will fail if we rename ind_active to active_cells. Since we are going to remove that function in two minor releases anyways I don't think we should care about that.

I also bumped the minimum required version of discretize to v0.12, since we need to use the point2index method here.

I left a few comments below. Check them out and let me know what do you think.