# \[question\] Finding number of closed potential contours and area enclosed by them

**URL:** <https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034>\
**Category:** Visualization\
**Tags:** question, plotting, pluto, makie, glmakie\
**Created:** [December 20, 2024, 6:07pm UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034 "2024-12-20T18:07:51Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![hisacro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hisacro/32/202631_2.png) [@hisacro](https://discourse.julialang.org/u/hisacro)\
**Post date:** [December 20, 2024, 6:07pm UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/1 "2024-12-20T18:07:51Z")

</div>

Hello, I’m trying to calculate the number of closed contours in a potential as well the area enclosed by closed,

i.e. In this image below at certain energy level there are 4 closed contours. I’m trying to figure a way to calculate the number of closed contours (4 here) and the area enclosed by them, given the energy levels.

 ![con_1](https://global.discourse-cdn.com/julialang/original/3X/a/7/a7958f02718c45c8fc1aa63e77c38436e84f4d65.png)

sample code for above (adopted from Pluto.jl notebook example),

```julia
using GLMakie

len = 50
x = range(-1, 1, length=len)
y = range(-1, 1, length=len)
four_wells(x,y) = x^4 - x^2 + y^4 - y^2
z = [four_wells(x[i], y[j]) for i in 1:len, j in 1:len]

# energy
height = -0.4

f = Figure(;resolution=(680,400))
ax1 = Axis3(f[1, 1])
ax2 = Axis(f[1, 2])

# Draw the contour in the 3D graph
plane = [height for i in 1:length(x), j in 1:length(y)];

# Drawing the function and the contour on the 3G graph
surface!(ax1, x, y, z, colormap=:viridis)
surface!(ax1, x, y, plane; transparency=true, colormap=[:white,:white])

# Drawing the function and the countour on the heatmap
hm = contourf!(ax2, x, y, z)
Colorbar(f[:, 3], hm; vertical=true)
contour!(ax2, x, y, z; linewidth=3, levels=[height], color=:white)
f

```

> **Pluto notebook file for the same**
>
> ```julia
> ### A Pluto.jl notebook ###
> # v0.19.35
> 
> using Markdown
> using InteractiveUtils
> 
> # This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
> macro bind(def, element)
> quote
> local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
> local el = $(esc(element))
> global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
> el
> end
> end
> 
> # ╔═╡ ca18be6d-8a9b-4c53-bd2e-553ea01ebb13
> using PlutoUI, GLMakie
> 
> # ╔═╡ 503aa0a0-ece8-4acc-b437-3119055d22bc
> k = @bind height PlutoUI.Slider(-1:0.1:1; default=0, show_value=true)
> 
> # ╔═╡ 2dab2fcd-77bf-4651-8223-6854863f425a
> md"""
> # Appendix
> """
> 
> # ╔═╡ 58d43864-9a74-42da-96d5-afb59f22cd0f
> begin
> len = 50
> x = range(-1, 1, length=len)
> y = range(-1, 1, length=len)
> four_wells(x,y) = x^4 - x^2 + y^4 - y^2
> z = [four_wells(x[i], y[j]) for i in 1:len, j in 1:len]
> end
> 
> # ╔═╡ e39373bd-8e12-469f-bd66-2d005ff6cff2
> begin
> # Setup of the figure
> f = Figure(;resolution=(680,400))
> ax1 = Axis3(f[1, 1])
> ax2 = Axis(f[1, 2])
> 
> # Draw the contour in the 3D graph
> plane = [height for i in 1:length(x), j in 1:length(y)];
> 	
> # Drawing the function and the contour on the 3G graph
> surface!(ax1, x, y, z, colormap=:viridis)
> surface!(ax1, x, y, plane; transparency=true, colormap=[:white,:white])
> 
> # Drawing the function and the countour on the heatmap
> hm = contourf!(ax2, x, y, z)
> Colorbar(f[:, 3], hm; vertical=true)
> contour!(ax2, x, y, z; linewidth=3, levels=[height], color=:white)
> end;
> 
> # ╔═╡ 9c17c0d8-e4cb-41c0-b180-b18a3c87e265
> # Displaying the function
> f
> 
> # ╔═╡ Cell order:
> # ╠═503aa0a0-ece8-4acc-b437-3119055d22bc
> # ╟─9c17c0d8-e4cb-41c0-b180-b18a3c87e265
> # ╠═e39373bd-8e12-469f-bd66-2d005ff6cff2
> # ╟─2dab2fcd-77bf-4651-8223-6854863f425a
> # ╠═ca18be6d-8a9b-4c53-bd2e-553ea01ebb13
> # ╠═58d43864-9a74-42da-96d5-afb59f22cd0f
> 
> ```

There’s a python implementation on [stackoverflow](https://stackoverflow.com/questions/22678990/how-can-i-calculate-the-area-within-a-contour-in-python-using-the-matplotlib), any pointers on which packages to look on Julia will be more helpful.

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [December 20, 2024, 8:58pm UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/2 "2024-12-20T20:58:20Z")

</div>

I’m not familiar with the internals of Makie’s contour code, so I can’t offer a solution as elegant as the python equivalent. Nevertheless, here is an idea for how to do a the calculation using other Julia tools:

- To obtain a representation of the contour, MDBM.jl will be able to sample the iso-surface of the potential
- The list of points given by MDBM.jl will not distinguish the number of closed contours, so you will have to cluster the points into disconnected contours. It should be fairly easy to write a function for this if there aren’t any packages available.
- With the list of points in each contour, you can sort these points in a counterclockwise fashion and calculate their area, e.g. using the [shoelace formula](https://en.wikipedia.org/wiki/Shoelace_formula) as in the python example

I hope these ideas help

---

<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:** [December 20, 2024, 11:16pm UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/3 "2024-12-20T23:16:41Z")

</div>

Funny, this seems the fill pits steps in DEMs. Found this package [WhereTheWaterFlows](https://github.com/mauro3/WhereTheWaterFlows.jl) that has to solve it as an intermediate step.

---

<div class="post-metadata">

**Author:** ![hisacro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hisacro/32/202631_2.png) [@hisacro](https://discourse.julialang.org/u/hisacro)\
**Post date:** [December 22, 2024, 3:20am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/4 "2024-12-22T03:20:47Z")

</div>

Another implementation in python, [shapely · PyPI](https://pypi.org/project/shapely/).

```julia
from shapely.geometry import Polygon
ct = [(0,0), (0,1), (1,1), (1,0)]
polygon = Polygon(ct)
>>> polygon.area
1.0

```

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [December 22, 2024, 10:11am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/5 "2024-12-22T10:11:28Z")

</div>

The DEM-filling in WhereTheWaterflows will fill the local minima completely. In above example it would fill those dips with “water” until their shores touch. So, probably not the right thing. But it should be possible to extract the contours from the plotting datastructure somehow… Ah, Makie calls into [GitHub - JuliaGeometry/Contour.jl: Calculating contour curves for 2D scalar fields in Julia](https://github.com/JuliaGeometry/Contour.jl) .

---

<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:** [December 22, 2024, 11:48am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/6 "2024-12-22T11:48:10Z")

</div>

Yes, but if we then subtract the filled and non filled grids the zero contour of that difference grid will be what is seemed here.

---

<div class="post-metadata">

**Author:** ![hisacro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hisacro/32/202631_2.png) [@hisacro](https://discourse.julialang.org/u/hisacro)\
**Post date:** [December 22, 2024, 6:06pm UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/7 "2024-12-22T18:06:41Z")

</div>

> [@mauro3](#):
>
> . In above example it would fill those dips with “water” until their shores touch

In fact I’m trying to get that point. So, what I understood from DEM-filling is, if I starting filling water in a local minima and it fills gradually till it reaching the brim and then a sudden increase in volume.

Will I be able to plot volume filled with respect to height capturing the above behaviour?

---

<div class="post-metadata">

**Author:** ![maltezfaria](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maltezfaria/32/7728_2.png) [@maltezfaria](https://discourse.julialang.org/u/maltezfaria)\
**Post date:** [December 23, 2024, 8:18am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/8 "2024-12-23T08:18:38Z")

</div>

For the area you could use ImplicitIntegration.jl:

```julia
using GLMakie
using ImplicitIntegration

len = 50
x = range(-1, 1; length = len)
y = range(-1, 1; length = len)
four_wells(x, y) = x^4 - x^2 + y^4 - y^2
z = [four_wells(x[i], y[j]) for i in 1:len, j in 1:len]

# energy
height = -0.4

f = Figure(; resolution = (680, 400))
ax1 = Axis3(f[1, 1])
ax2 = Axis(f[1, 2])

# Draw the contour in the 3D graph
plane = [height for i in 1:length(x), j in 1:length(y)];

# Drawing the function and the contour on the 3G graph
surface!(ax1, x, y, z; colormap = :viridis)
surface!(ax1, x, y, plane; transparency = true, colormap = [:white, :white])

# Drawing the function and the countour on the heatmap
hm = contourf!(ax2, x, y, z)
Colorbar(f[:, 3], hm; vertical = true)
contour!(ax2, x, y, z; linewidth = 3, levels = [height], color = :white)

# compute a quadrature for the implicit region ϕ<0
ϕ = (x) -> four_wells(x[1],x[2]) - height # level-set
quad, logger = quadgen(ϕ, (-1, -1), (1, 1); order = 3)
println("Estimated area: ", integrate(x -> 1, quad))
# visualize the quadrature nodes
scatter!(ax2, quad.coords; markersize = 6)

f

```

But I don’t think you can easily obtain the connected components using ImplicitIntegration.jl

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [December 24, 2024, 3:26am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/9 "2024-12-24T03:26:18Z")

</div>

Here’s how you can get the line coordinates with Contours.jl:

```julia
using Contours

len = 50
x = range(-1, 1; length = len)
y = range(-1, 1; length = len)
four_wells(x, y) = x^4 - x^2 + y^4 - y^2
z = [four_wells(x[i], y[j]) for i in 1:len, j in 1:len]

# energy
height = -0.4

cls = contours(x, y, z, [height]) # compute the contours at manually specified levels
level = only(levels(cls)) # extract the only result (since we only asked for one level)

points = [l.vertices for l in lines(level)] # ::Vector{Vector{Tuple{Float64, Float64}}}

```

Now that we have the lines, we need to check that they’re closed. The nature of Contour.jl’s algorithm makes this simple - we only need to know whether the first and last points are the same.

```julia
closed_ring_points = filter(points) do pointvec
    pointvec[begin] == pointvec[end] # you can also make this approximate, if you want
end

```

Now that you have polygon-like constructs, it’s easy enough to solve for its area. There are multiple solutions, I’m using GeoInterface.jl and GeometryOps.jl here.

```julia
import GeoInterface as GI, GeometryOps as GO
# construct geometry from vectors of points
rings = GI.LinearRing.(closed_ring_points)
polygons = [GI.Polygon([ring]) for ring in rings]

areas = GO.area.(polygons) # 4x ~0.1694

# if you want to be efficient / nonallocating
total_area = sum(GO.area, polygons) # 0.6775172177537117

```

Hope that helps! If you are worried about intersecting polygons, you can also use GeometryOps.jl to fix that:

```julia
mp = GI.MultiPolygon(polygons)

fixed_mp = GO.fix(mp; corrections = [GO.UnionIntersectingPolygons()])

```

---

<div class="post-metadata">

**Author:** ![hisacro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hisacro/32/202631_2.png) [@hisacro](https://discourse.julialang.org/u/hisacro)\
**Post date:** [December 24, 2024, 4:40am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/10 "2024-12-24T04:40:05Z")

</div>

Mind linking Contours.jl repo, the link you provided isn’t pointing to the package.

Edit: It seems discourse auto searches the package, again in my local environment too I can find only `Contour` not `Contours` and it doesn’t have the function `contours`.

---

<div class="post-metadata">

**Author:** ![raman\_kumar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raman_kumar/32/26782_2.png) [@raman\_kumar](https://discourse.julialang.org/u/raman_kumar)\
**Post date:** [December 24, 2024, 5:01am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/11 "2024-12-24T05:01:40Z")

</div>

This is correct github link of [Contour.jl](https://github.com/JuliaGeometry/Contour.jl) and this is link of [Contour.jl documentation](https://juliageometry.github.io/Contour.jl/stable/). See contours() function in this image inside _levels()_.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/d/9d95c6535ad8a26aabfec624e4025e11c6ae14de.png)

---

<div class="post-metadata">

**Author:** ![hisacro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hisacro/32/202631_2.png) [@hisacro](https://discourse.julialang.org/u/hisacro)\
**Post date:** [December 24, 2024, 5:03am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/12 "2024-12-24T05:03:27Z")

</div>

but that doesnt provide the function `contours`

```julia
julia> Contour.
ContourCollection ContourLevel Curve2 E
EN ES EW N
NE NESW NS NW
NWSE S SE SN
SW W WE WN
WS _get_case advance_edge chase!
contour contourlevels contours coordinates
edge_LUT eval get_first_crossing get_level_cells
get_next_edge! include interpolate level
levels lines next_edge next_map
trace_contour vertices

```

edit 1,

```julia
julia> @which contours(x, y, z, [height])
contours(x, y, z, levels; VT)
     @ Contour ~/.julia/packages/Contour/vefVR/src/Contour.jl:94

```

seems to be [Contour.jl/src/Contour.jl at master · JuliaGeometry/Contour.jl · GitHub](https://github.com/JuliaGeometry/Contour.jl/blob/master/src/Contour.jl)

I have manually included the above Contour.jl and interpolate.jl, seems to work now but I wonder why `using Contour` is not providing the function.

---

<div class="post-metadata">

**Author:** ![asinghvi17](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/asinghvi17/32/8272_2.png) [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)\
**Post date:** [December 25, 2024, 7:28am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/13 "2024-12-25T07:28:03Z")

</div>

Ah yeah my bad, it should have been Contour.jl.

Can you confirm the version of the package? It should be the latest.

---

<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:** [January 5, 2025, 8:37pm UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/14 "2025-01-05T20:37:12Z")

</div>

Still interested on this? With GMT.jl master we can now do.

```julia
G = peaks(N=128);
G2 = fillsinks(G);
viz(G2, shade=true)

```

 ![filada](https://global.discourse-cdn.com/julialang/original/3X/b/b/bb9f929b423ca204e28a3b3aff874335b5ebded6.jpeg)

---

<div class="post-metadata">

**Author:** ![hisacro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hisacro/32/202631_2.png) [@hisacro](https://discourse.julialang.org/u/hisacro)\
**Post date:** [January 6, 2025, 5:22am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/15 "2025-01-06T05:22:27Z")

</div>

Can you elaborate a bit on how G is defined and how to get area by using fillsinks.

I couldn’t find that function on the [documentation](https://www.generic-mapping-tools.org/GMTjl_doc/search/index.html?q=fillsinks) page.

---

<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:** [January 6, 2025, 11:09am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/16 "2025-01-06T11:09:28Z")

</div>

The HTML documentation has not yet been generated but the online dons do exist. Just do the normal `? fillsinks`  
This function is meant to fill the _sinks_, not to measure the area. But we can that as well. There is an option to also save the contours (which are not exactly contours, see the docs about the 256 levels) in a global variable.

```julia
G = peaks(N=128);
G2 = fillsinks(G, saco=true);
D = get(GMT.SACO[1], "saco", nothing);

gmtspatial(D, area=true)
BoundingBox: [-1.831659459383318, 0.31258763887390784, -1.923463287251611, 0.2138927839499514, 0.3258726517453034, 4.3948167896335795]

3×3 GMTdataset{Float64, 2}
 Row │ centroid_x centroid_y area
─────┼──────────────────────────────────
   1 │ 0.312588 0.213893 0.325873
   2 │ -1.83166 0.201386 3.74307
   3 │ 0.293074 -1.92346 4.39482

```

> [@hisacro](#):
>
> Can you elaborate a bit on how G is defined

`G` is just a simple demo grid made with the _peaks_ function. Note that for you to use this you must convert your array into a `GMTgrid` type. The function `mat2grid` will help on that.

---

<div class="post-metadata">

**Author:** ![hisacro](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hisacro/32/202631_2.png) [@hisacro](https://discourse.julialang.org/u/hisacro)\
**Post date:** [January 10, 2025, 5:38am UTC](https://discourse.julialang.org/t/question-finding-number-of-closed-potential-contours-and-area-enclosed-by-them/124034/17 "2025-01-10T05:38:20Z")

</div>

Found a cool guide about visualization of contour area using GLMakie.

> **[Julia and Projective Geometric Algebra](https://olarth.medium.com/julia-and-projective-geometric-algebra-e2d5abb11519)**
>
> experience the geometry superpowers of PGA