How to interpolate geoid grid for reference heights

[geojoeyphillips]

Hello Fatiando group!

I really enjoy the python packages you all have setup for processing gravity data. I have been going through the documentation and using your tools to work on my own data. The trouble I'm having is referencing my measurement points to the ellipsoid. From your old tutorial I'm talking about this step:

topography_geometric = topography_proj + geoid_proj

My region is rather small so when I do:

region_pad = vd.pad_region(rvdc_region,pad=5)
data = data[vd.inside((data.Longitude,data.Latitude),rvdc_region)]
geoid = geoid.sel(longitude=slice(*region_pad[:2]), latitude=slice(*region_pad[2:]))
topography = topography.sel(lon=slice(*region_pad[:2]), lat=slice(*region_pad[2:]))

My topography has dimensions (601,541) while the geoid is (60,60) leading to (0,0) for topography_geometric. Is there a way to interpolate the geoid grid to (601,541) or interpolate my topography grid to (600,540) so I can interpolate the geoid (6010,609)?

Thank you for your help!

[leouieda]

Hi @geojoeyphillips thanks for posting this!

topography_geometric = topography_proj + geoid_proj only works if the two grids are co-registered/co-located. If the coordinates don't line up, xarray can't do the operation and wind up with an empty grid.

The way to do this is to interpolate the geoid array onto the same coordinates as the topography array. This could be done with Verde but for grids the easiest way is to use the xarray interpolation methods:

# fill_value used to extrapolate and avoid having NaNs in the grid
geoid_interp = geoid_proj.interp_like(topography_proj, kwargs={"fill_value": None}) 
topography_geometric = topography_proj + geoid_interp

Be sure to plot the geoid_interp to verify that it looks sensible.

Hope this helps! Let me know the outcome 😀

[geojoeyphillips]

Thank you so much! Much easier than my failed attempt:

f = interp2d(geoid_proj.easting,geoid_proj.northing,geoid_proj,kind = 'cubic')
geoid_new = f(topography_proj.easting,topography_proj.northing)

geoid_xr = xr.DataArray(geoid_new,coords = {'northing':topography_proj.northing,'easting':topography_proj.easting},dims=['northing','easting'])

Thank you again! Seems like with your geoid_interp method I can project back to lat,long then plot with pygmt.

If I can ask another question, can I grid the gravity_disturbance_mgal = data.gravity - data.normal_gravity_mgal the same way? Or do I need to go through https://www.fatiando.org/verde/latest/gallery/scipygridder.html after making gravity_disturbance_mgal into a grid?

Thank you!

[leouieda]

Thank you so much! Much easier than my failed attempt:

I always think xarray shines when you're given an xarray.DataArray (either coming from Verde or loading from a file). Creating them is always a bit of a pain and not very intuitive.

can I grid the gravity_disturbance_mgal = data.gravity - data.normal_gravity_mgal the same way?

You can but you probably don't want to do that (and don't have to). The normal gravity is a calculated thing so instead of calculating on certain points and then interpolating, you can calculate it directly on the points you need. This usually means on the point where you have observed gravity. If you're using a gravity grid from a spherical harmonic model, then you can pass in latitude and height from the grid to boule.WGS84.normal_gravity to get an array of normal gravity. See https://www.fatiando.org/boule/v0.4.1/user_guide/gravity_disturbance.html

---

PS: Here are some general tips for using this forum and GitHub 🚀

1. You can format things with Markdown syntax. That's how I make my text look like code blocks and things like this = code formatted. It really helps everyone to read a code block if it's formatted as such. See this tutorial from GitHub. And don't worry about missing it or forgetting to do it. We're here to help and not to judge 🙂 (I edited you original post to include that)
2. If an answer solves your problem, please click on "Mark as answer" above ☝🏾 to let people know that your question as answered and what the answer was.