Hi everyone. My data has 4 columns: lat, Lon, depth, velocity, and I am trying to interpolate velocity for different depth layers. I used gmt surface but it was unsuccessful. The reason is after using the surface command, the output has only 3 columns.
I tried this:
gmt blockmean MIT_myarea.csv -: -R-105/-80/5/25 -I0.5 -W > mean.txt #the out put has 4 col
gmt surface mean.txt -Gseis.grd -I0.5 -: -R-105/-80/5/25
gmt grd2xyz seis.grd > interpolated_velocity.txt #the output has 3 col
Why is the 4th column gone? Any help would be appreciated. Thanks for your time.
The only 3-D interpolator in GMT is greenspline. Are your data already equidistant (so you just want to reformat into a grid) or do you truly need interpolation because your points are more or less randomly organized? Anyway, see greenspline for 3-D gridding.
Depends. If you intend to just reformat the data into a 3-D netcdf data cube with variable depth spacing then there is no interpolation needed. If you instead want to have layers every dz km then you need interpolation. GMT is not yet well suited to do this. We just recently added some support for 3-D cubes and have grdinterpolate which can work on cubes and grdinfo understand cubes. If you do need interpolation then greenspline should do that but it probably will need the number of data points to be managable (10000 or less) since it solves a nxn matrix.
Thank you so much.
I need to find velocities in depths that the model doesn’t have. So, I think I have to use greenspline (maybe several times to cover all data (~ 20500)).
I’d like to ask another question, Sir.
The greenspline’s output is not an ASCII file (I thought I should use grd2xyz to convert it to a text file but I am not sure if this is the right way). Could you please let me know how I can see the results? I need the interpolated data in a text file to enter as an output in a tomography package. Thanks.
Id like to ask a question. To my understanding, if I want to extract a slice from the output of greenspline at a desired depth (let’s call it D), I should assign D to -T. This way, I have a slice of earth at D. I was wondering if you could let me know if I am correct. Thanks.
and the numbers of the column including the layer (velocity value) do not go beyond 7.8 after implementing “grdinterpolate” which is incorrect. For deeper layers it should be more than 9 km/s. Could you please share your thoughts on this with me? thanks in advance for your help.
I don’t have experience with greenspline to catch troubles just by eye. But in your first post you said that you already had your data interpolated in lat, Lon, depth, is that right?
Thanks for your reply. I have my data with equidistant x and y in different layers (not equidistant depths) with velocity values at these different layers. And I want to acquire velocity values in arbitrary depths. That is why I need interpolation on my data.
OK, so that’s what I said here (or in another similar thread). Not knowing why the greensplineis failing (it with many points it’s a very demanding method) the alternatives are.
With plain GMT: create individual grids, one by each depth layer, and stack them together with grdinterpolate
Use the GMT.jl Julia wrapper. With it it would be (hopefully) just do something like Cube = xyzw2cube("/path/to/yourfile"); Next, either save it with gmtwrite(Cube, "cube.nc") and go back to plain GMT. Or use G = cubeslice(Cube, newdepth) and get that desired (linear) interpolated layer.
This is an important suggestion. I have gridded files at several pressure layers from ERA5 and I want to interpolate at the desired pressure level. Joaquim’s suggestion could work also for me.