# Locate Points within GMT Delauney triangularization

**URL:** <https://discourse.julialang.org/t/locate-points-within-gmt-delauney-triangularization/110967>\
**Category:** General Usage\
**Tags:** gmt, voronoidelaunay\
**Created:** [February 29, 2024, 4:52pm UTC](https://discourse.julialang.org/t/locate-points-within-gmt-delauney-triangularization/110967 "2024-02-29T16:52:03Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![morningglory05](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/morningglory05/32/43282_2.png) [@morningglory05](https://discourse.julialang.org/u/morningglory05)\
**Post date:** [February 29, 2024, 4:52pm UTC](https://discourse.julialang.org/t/locate-points-within-gmt-delauney-triangularization/110967/1 "2024-02-29T16:52:03Z")

</div>

**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:

```julia
# 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.

---

<div class="post-metadata">

**Author:** ![empet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/empet/32/221303_2.png) [@empet](https://discourse.julialang.org/u/empet)\
**Post date:** [February 29, 2024, 8:17pm UTC](https://discourse.julialang.org/t/locate-points-within-gmt-delauney-triangularization/110967/2 "2024-02-29T20:17:20Z")

</div>

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:

```julia
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)) 

```

![triangulatedsurf](https://global.discourse-cdn.com/julialang/original/3X/9/e/9ee6a769866bd5424271566eb94aa90aa4d7e863.png)

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [February 29, 2024, 10:28pm UTC](https://discourse.julialang.org/t/locate-points-within-gmt-delauney-triangularization/110967/3 "2024-02-29T22:28:38Z")

</div>

> [@empet](#):
>
> If the answer is positive, then I don’t think that GMT provides a function to plot it.

Actually it does

```julia-auto
x,y,z = GMT.peaks(N=45, grid=false);
trisurf([x[:] y[:] z[:]], pen=0.5, show=true)

```

 ![GMTjl_j-fs8](https://global.discourse-cdn.com/julialang/original/3X/d/8/d8401db07e0f6a1131f7992884d204c3ebbf3cc0.png)

> [@morningglory05](#):
>
> is there an easy way to know which triangle/edge it is in/on?

Maybe [inwhichpolygon](https://www.generic-mapping-tools.org/GMTjl_doc/documentation/utilities/inwhichpolygon/index.html#inwhichpolygon)?

---

<div class="post-metadata">

**Author:** ![morningglory05](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/morningglory05/32/43282_2.png) [@morningglory05](https://discourse.julialang.org/u/morningglory05)\
**Post date:** [April 30, 2024, 3:49pm UTC](https://discourse.julialang.org/t/locate-points-within-gmt-delauney-triangularization/110967/4 "2024-04-30T15:49:56Z")

</div>

Sorry for the delayed response.

`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.

Thanks for the solution!

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [April 30, 2024, 5:53pm UTC](https://discourse.julialang.org/t/locate-points-within-gmt-delauney-triangularization/110967/5 "2024-04-30T17:53:09Z")

</div>

`triplot` calls `triangulate` so they are both using the same triangulation (from Shewchuk [Triangle](http://www.cs.berkeley.edu/~jrs))

You can compute contours directly from the triangulation. See the module [contour](https://www.generic-mapping-tools.org/GMTjl_doc/documentation/modules/contour/index.html#contour)