# Plotting surface without interpolating at discontinuity

**URL:** https://discourse.julialang.org/t/plotting-surface-without-interpolating-at-discontinuity/95279
**Category:** General Usage
**Tags:** question, plotting, plots
**Created:** [February 27, 2023, 3:58pm UTC](https://discourse.julialang.org/t/plotting-surface-without-interpolating-at-discontinuity/95279 "2023-02-27T15:58:04Z")
**Posts on this page:** 6
**Page:** 1

<div class="post-metadata">

### Author: ![NoFishLikeIan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nofishlikeian/32/20917_2.png) [@NoFishLikeIan](https://discourse.julialang.org/u/NoFishLikeIan)
#### Post date: [February 27, 2023, 3:58pm UTC](https://discourse.julialang.org/t/plotting-surface-without-interpolating-at-discontinuity/95279/1 "2023-02-27T15:58:04Z")

</div>

I am attempting to plot a function that has a discontinuity, consider for example,

```julia
function f(x, y)
    if x^2 + y^2 ≤ 1/2
        return 1.
    else
        return 0.
    end
end

```

If I use the `surface` function from `Plots.jl` I obtain the following plot

```julia
surface(-1:0.01:1, -1:0.01:1, f)

```

![test](https://global.discourse-cdn.com/julialang/original/3X/c/3/c300dfa83eb6e9c31b8840e254b65551e6db772d.png)

Is there a way to plot the function without `surface` automatically interpolating the discontinuity points?

---

<div class="post-metadata">

### Author: ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [February 27, 2023, 5:56pm UTC](https://discourse.julialang.org/t/plotting-surface-without-interpolating-at-discontinuity/95279/2 "2023-02-27T17:56:20Z")

</div>

I don’t know of an automatic method, but in general, to obtain the same result, `NaN`s are inserted around the discontinuities to prevent interpolation:

 ![Plots_plotlyjs_surface_discontinuity](https://global.discourse-cdn.com/julialang/original/3X/5/b/5b346ed2466aee1b1d76344f646190da5087d0b2.png)

> **Plots.jl code**
>
> ```julia
> using Plots; plotlyjs()
> 
> function f(x, y, δ)
> r = x^2 + y^2
> if r ≤ 1/2 - δ
> return 1.0
> elseif r ≥ 1/2 + δ
> return 0.0
> end
> return NaN
> end
> 
> δ = 0.01
> x = y = -1.0:δ:1.0
> Plots.surface(x, y, (x,y)->f(x,y,δ))
> 
> # QC mask of NaN:
> heatmap(x, y, (x,y)->f(x,y,δ), ratio=1, lims=(-1,1))
> 
> ```

---

<div class="post-metadata">

### Author: ![NoFishLikeIan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nofishlikeian/32/20917_2.png) [@NoFishLikeIan](https://discourse.julialang.org/u/NoFishLikeIan)
#### Post date: [February 27, 2023, 6:12pm UTC](https://discourse.julialang.org/t/plotting-surface-without-interpolating-at-discontinuity/95279/3 "2023-02-27T18:12:16Z")

</div>

Thank you for your answer!

Yeah, I thought about this, but in my case the discontinuity points are hard to compute. They are determined by an implicit equation, which I would have to solve with a root finder. I can do that, but it is not that elegant.

---

<div class="post-metadata">

### Author: ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [February 27, 2023, 7:30pm UTC](https://discourse.julialang.org/t/plotting-surface-without-interpolating-at-discontinuity/95279/4 "2023-02-27T19:30:42Z")

</div>

Ok, check if this automatic method could work in your case:

```julia
using Plots; plotlyjs(size=(800,600), dpi=600)

f(x, y) = x^2 + y^2/4 ≤ 1/2 ? 1. : 0.

x = -1.0:0.01:1.0
y = -2.0:0.01:2.0
z = f.(x,y')

ϵ = 0.5
∇x = vcat(diff(z, dims=1), zeros(1, length(y)))
∇y = hcat(diff(z, dims=2), zeros(length(x)))
z[hypot.(∇x, ∇y) .> ϵ] .= NaN

surface(x, y, z', xlabel="x", ylabel="y", zlabel="z")

```

_ **NB:** changed the circular disc to an ellipse for easier debugging_

 ![Plots_plotlyjs_discontinuity_automatic](https://global.discourse-cdn.com/julialang/original/3X/6/4/64ed07e84f4cd559195eb9c68876ad51c5db4186.png)

---

<div class="post-metadata">

### Author: ![jheinen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jheinen/32/2787_2.png) [@jheinen](https://discourse.julialang.org/u/jheinen)
#### Post date: [March 1, 2023, 9:49am UTC](https://discourse.julialang.org/t/plotting-surface-without-interpolating-at-discontinuity/95279/5 "2023-03-01T09:49:22Z")

</div>

… should also work with GR (or Plots with the default backend):

 ![Screenshot 2023-03-01 at 10.47.35](https://global.discourse-cdn.com/julialang/original/3X/5/d/5d88c619b0b060e09a01ed24db34eb86cfb730ee.png)

---

<div class="post-metadata">

### Author: ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [March 1, 2023, 10:15am UTC](https://discourse.julialang.org/t/plotting-surface-without-interpolating-at-discontinuity/95279/6 "2023-03-01T10:15:26Z")

</div>

@jheinen, Plots.jl with gr() backend seems to exhibit several bugs on this example.

Bug on `aspect_ratio=:equal`, which is not applied to the z-axes:

```julia
surface(x, y, z', ratio=:equal)

```

Bug on the color scale limits if we set `zlims` to fix the `aspect_ratio` bug above:

```julia
surface(x, y, z', ratio=:equal, zlims=(-2,2))

```
