Overall issue: I am finding the intersection of two planes. I have the values of both planes at ungridded points, so I subtract their values (call the difference z). Then I have x,y,z values and I find the contour for z=0. I then plan to take this contour and map it back to my original surfaces. There are several ways to accomplish that last step.
Specific issue: This is the code I have so far:
# random data to simulate z
f(x, y) = sin(x + 10sin(x)) + cos(y)
xx = rand(10000);
yy = rand(10000);
d = hcat(hcat(xx,yy),f.(xx,yy));
# create Delauny triangularization
using GMT
tri = GMT.triangulate(d) #the output are index numbers of the original points based on line
# create contour at z=0
cont = GMT.contour(d,index=tri, dump=true, contour=[0]) #this gives x,y,z values
Given a specific point in cont, say cont[1][1,1:2], is there an easy way to know which triangle/edge it is in/on?
I believe that I can code a search through all the triangles, but I feel like there should be an easier way. I can also take the x,y points in the contour and plug them back into my functions (not necessarily what i want to do though).
Most of the documentation I’m reading has to do with plotting or output of tri and cont, but I don’t know how to effectively use the structure of tri in anything but plotting.
Do you want to plot a trisurf, i.e. a triangulated surface of equation z=f(x,y), (x,y)\in [a,b]\times[c,d]?
If the answer is positive, then I don’t think that GMT provides a function to plot it. Instead you can use DelaunayTriangulations.jl to triangulate the rectangle of definition.
Namely, you should define a mesgrid over that rectangle
and triangulate the set of meshgrid points, (x,y), stored in a vector pts.
Then lifting the vertices of the Delaunay triangulation to the surface, i.e. to
3d points, (x,y, f(x,y)), and keeping the triangles resulted from the planar points triangulation, we get a surface triangulation that can be plotted as Plots.mesh3d:
Example:
using DelaunayTriangulation
using Plots, ColorSchemes
plotlyjs();
f(x, y) = sin(x) + cos(y)
m = 75
n = 75 #n can be different from m
xl = range(0, 2pi, n) # or xl= 2pi*rand(n)
yl=range(0, 2pi, m) # yl= 2pi*rand(m)
x = vec(ones(m)*xl')
y = vec(yl *ones(n)')
z = f.(x, y)
pts = [[xe, ye] for (xe, ye) in zip(x, y)]
tri = triangulate(pts)
cm = stack(tri.triangles, dims=1) #cm is a Matrix of size (N, 3)
plt = plot(seriestype=:mesh3d, x, y, z,
connections=tuple([(c .-1) for c in eachcol(cm)]...), #mesh3d works with tri inds starting from 0
legend=false, intensity=z,
seriescolor = reverse(cgrad(ColorSchemes.matter.colors)), size=(400, 400))
inwhichpolygon works, but only if my input is Vector{GMTdataset}. I can get this when I triangulate using triplot. (Yay solution)
On the other hand, I also need to compute the contours. Currently, I am using GMT.contour for this which, if I want to tell it which triangularization to use, requires the input be a GMTdataset of indices. This can be accomplished with triangulate instead of triplot. Thus, I end up calculating the Delauney triangularization of the same data set twice. Though it is not a problem for me right now, it is annoying that (1) I can’t be certain that the triangualizations are the same and (2) I have to use extra computation power to do it twice.
The only solution I can forsee is if there is a method of computing contours on a nonuniform grid that accepts triplot as compatible triangularization. As it stands, I’m just going to compute the triangularization twice.
