I’m trying to interpolate a 2D surface. surface is working well, but now I’ve got some additional gradient constraints and gave greenspline a try. My call is
where depth.dat has the form latitude, longitude, depth, e.g.
-125.523 49.469 26.8897
-123.988 49.227 40.2961
-124.697 49.000 27.7391
-125.498 48.961 19.0781
-124.262 48.901 33.8943
and grad.dat has the form latitude, longitude, dip-direction, gradient (i.e. tan(dip)), e.g.:
-125.523 49.469 52.1 0.160174
-123.988 49.227 28.9 0.249328
-124.697 49.000 43.2 0.151236
-125.498 48.961 47.6 0.277325
-124.262 48.901 55.8 0.367928
I get the error message:
greenspline [ERROR]: Input data have 3 column(s) but at least 4 are needed
N.B.: I’m on gmt 6.4.0 and without the gradient file, it works just fine:
gmt greenspline depth.dat -Gout.grd -St0.25 -R-126/-122/43/50 -I0.1 -Z2
Does anyone have a clue about what’s going wrong here?
Thanks for looking into this swiftly, looking forward to getting it via conda. Would it be possible that the fix be included in a release before the end of the year?
Unfortunately the above command fails on my machine with:
UnsatisfiableError: The following specifications were found to be incompatible with each other:
Output in format: Requested package -> Available versions
Package libstdcxx-ng conflicts for:
gmt=6.5.0.dev2+b0c2f6d -> geos[version='>=3.11.0,<3.11.1.0a0'] -> libstdcxx-ng[version='>=10.3.0|>=12|>=7.5.0|>=9.3.0']
conda-forge/linux-64::python==3.10.5=h582c2e5_0_cpython -> libffi[version='>=3.4.2,<3.5.0a0'] -> libstdcxx-ng[version='>=11.2.0|>=7.5.0']
Package libgcc-ng conflicts for:
gmt=6.5.0.dev2+b0c2f6d -> fftw[version='>=3.3.10,<4.0a0'] -> libgcc-ng[version='>=10.3.0|>=9.4.0|>=7.5.0|>=9.3.0|>=11.2.0']
gmt=6.5.0.dev2+b0c2f6d -> libgcc-ng[version='>=12']
Package libzlib conflicts for:
gmt=6.5.0.dev2+b0c2f6d -> hdf5[version='>=1.12.2,<1.12.3.0a0'] -> libzlib[version='>=1.2.12,<1.3.0a0']
gmt=6.5.0.dev2+b0c2f6d -> libzlib[version='>=1.2.13,<1.3.0a0']
Package zlib conflicts for:
gmt=6.5.0.dev2+b0c2f6d -> libnetcdf[version='>=4.8.1,<4.8.2.0a0'] -> zlib[version='>=1.2.12,<1.3.0a0']
gmt=6.5.0.dev2+b0c2f6d -> zlib[version='>=1.2.13,<1.3.0a0'
I see that you might not be able to fix this before the next point release, but still wanted to let you know. Perhaps I can come across some conda magic to fix this.
Hmm, those conflicts might be due to some other package in your existing conda environment. Usually I would end up deleting the entire conda environment and reinstalling the packages from scratch, but not sure if you want to go with that
Note that we are currently dealing with some issues with GEOS 3.11.1 and GDAL 3.6.0 at https://github.com/conda-forge/gmt-feedstock/pull/224, and that might help resolve some of the dependency conflicts, but it’ll take a few days to get ironed out. Edit: the migrations are done (thanks @seisman!) so you could give the conda install command a try again.
Thank you for your help again. I now got the gmt=6.5.0+dev2 running. For the record, for a environment with nothing installed at all, one needs to do the double step:
greenspline [WARNING]: Can only remove/restore mean z in mixed {z, grad(z)} data sets
and gmt grdinfo grd.out indicates all 0 z-values:
...
out.grd: v_min: 0 v_max: 0 name: z
...
The same issue existed when I tried a comparable command with gmt=5.x. Without gradient information, there is some output, though I didn’t check if it makes sense physically:
If you can post your depth.dat and grad.dat files we might be able to see what is going on. That the z and grads share lon/lat should not technically be a problem but may turn out to be a clue to a bug.
OK, I think it is safe to say at this point that the -A option does not work. We have no tests for it either. I am looking at the simpler 1-d spline and trying to impose a slope constraint and it does affect the curve but not as I expected (which I suspect has to do with normalisations), but I also found that if the slope constraint is given at a data point then I get no solution at all. So at the least you would have to perturb your slope locations away from the given points by some noise, but no guarantee that the solution is good. I will work on the 1-D problem to understand it first.
Some more notes on this: The -A option in green spline always assumed that these would be additional locations where slopes are available. We did not consider the case (yours, basically) where you have both value and slope given at specific locations. In that case, there should only be one parameter for that point (strength of the Green function there) and the same parameter is used with the directional derivative. In the -A case we have, these points are assumed separate and each gets their own unknown parameter to solve for. So we have N data plus M slope constraints = N + M parameters. Your case is actually N parameters with two constraints per point. This would require either green spline to determine which slope constraints occur at the same locations as data constrains, or a flag to provide both constraints on the same data record. Not straightforward in either case.
I can confirm that moving the support points of the gradient measurements away from the support points of the values produces erratic results. I tried to move them parallel and perpendicular to the azimuth, with small and larger increments.
What would be a possible timeline for a fix and can I help to advance it? I understand @pwessel is currently at AGU, and so am I. Perhaps we have a chance to discuss things briefly in person? If so, I’m more than happy to schedule a brief meeting. Please drop me an E-Mail in that case: wbloch@eoas.ubc.ca.
