Problem using Verde to project ETOPO file to cartesian coordinates

[hwil055]

Kia ora everyone,
I'm a geophysicist based at AUT (Auckland University of Technology) and this is my first time posting a question here :-)
I'm still in the early stages of using Python (and any of the Fatiando tools) so please excuse my ignorance....
I'm using the Transform 21 Harmonica tutorial to guide me through processing gravity data over Papua New Guinea and have gotten stuck trying to (use Verde to) project an ETOPO (2022) tile (.nc format, from NCEI/NOAA, loaded with xarray) to cartesian coordinates.
When I follow the steps in the tutorial (i.e. I run topo_plain = vd.project_grid(topography, projection), my RAM gets maxed out and my system hangs there and never resolves (I have to crash out).
I get the same result on another PC so I wonder if I'm not using the .nc file correctly. It's ~23 MB and I have no problem plotting it up before trying to project but it does appear to have a different format (e.g. variable names) to that used in the tutorial.
Is anyone who is familiar with the tutorial, able to offer some suggestions please?
Many thanks in advance,
Helen

[santisoler]

Hi @hwil055! Thanks for opening this question.

The vd.project_grid() function projects the grid points and then performs an interpolation into a new regular grid on the projected coordinates. When working with large grids these interpolations can demand a huge amount of RAM. This is specially true if you use the "linear" or "cubic" method, where linear and cubic splines are using for the interpolation, respectively. By default vd.project_grid() uses method="linear".

The good news is that for dense grids these interpolation methods are usually an overkill. You could use nearest neighbours instead and the results will still be accurate enough. So, for your case, pass method="nearest" in your vd.project_grid() function call:

vd.project_grid(topography, projection, method="nearest")

You can check the documentation of the vd.project_grid function here: https://www.fatiando.org/verde/latest/api/generated/verde.project_grid.html

Hope it works!

[hwil055]

Hello and thank you for your advice @santisoler ! Changing the method to nearest neighbour worked after about 4 hours of run time - at least on the machine with 16 Gb RAM - and I was able to finish the tutorial. I'm keen to explore the rest of the Fatiando ecosystem but I think I might need to upgrade my own laptop (12 Gb RAM) to do this properly, using the datasets/areas I'm interested in :-)

Kind regards,
Helen

[santisoler]

I'm glad it worked out! Although 4 hours seem a lot of time, but I imagine your grid is quite big since it's stored in a ~23MB netCDF file (which could include compression).

For the nearest method you don't need too much RAM, you really do with cubic and linear interpolators. I suspect that what is taking too long is the block reduction operation carried to avoid antialiasing. You could experiment with disabling the antialiasing passing antialias=False.

Interpolating large datasets it's a tricky problem. Hopefully we'll have a fast solution in Verde in the near future 🙂