Skip to content

Commit

Permalink
Change default interpolation to "nearest" to avoid logical error with…
Browse files Browse the repository at this point in the history 
… nan. Thumbs up of the merge - DF
  • Loading branch information
@fourndo
fourndo committed Jun 27, 2017
1 parent bdb4f2a commit 4890374
Showing 1 changed file with 1 addition and 1 deletion.
2 changes: 1 addition & 1 deletion SimPEG/Utils/modelutils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
import numpy as np
from scipy.interpolate import griddata, interp1d

def surface2ind_topo(mesh, topo, gridLoc='CC', method='cubic', fill_value=np.nan):
def surface2ind_topo(mesh, topo, gridLoc='CC', method='nearest', fill_value=np.nan):
"""
Get active indices from topography
"""
Expand Down

3 comments on commit 4890374

@thast
Copy link
Member

@thast thast commented on 4890374 Jun 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@fourndo : I am unsure of what is the best default for topogrpahy.
I chose cubic for the following reasons:

  • it is usually the go-to methods for interpolating topography in GIS...
  • the logical nan error is probably due to the fact that the topography does not cover the whole mesh. I think it can be fair to display an error when the provided topography does not cover fully the mesh (maybe arguable over the padding cells but for now...). Then the user knows it and can choose either to use the nearest interpolation method, or change its topography to fully cover the mesh.
  • nearest can give very blocky results if the topodata are very spaced out compared to the mesh cell sizes

@fourndo
Copy link
Member Author

@fourndo fourndo commented on 4890374 Jun 27, 2017 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I totally get your point @thast . And yes the issue was that the topo didn't cover all the corners of the mesh.

Would be great if griddata would also extrapolate on cubic and linear, but it doesn't. Issue with not having extrapolation is that nan are bad for the logical, and IF the topo is to short we end up with walls. Maybe we could add a print statement at some point to tell the user that the topo is too short.
Nearest is fast too, predictable and it has the advantage that you can technically specify a flat topo with a single point.

It's just for default anyway, user can be fancier if needed.

My two cents

@thast
Copy link
Member

@thast thast commented on 4890374 Jun 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a good point, and a different way to look at it too than mine.
@lheagy brought a good point: if we add a proper documentation, we do not have to worry that much about what the users are going to use as input in the function, and then fixing the default based on this prediction. Users can be fancier, but we need to provide the documentation for that. I will add it soon.

Please sign in to comment.