Problem with extracting specific contour from a ContourPlot

Question

I am trying to exctract contours of specific value (2.6*10^-6) out of a ContourPlot and measure their area.

Below you can see the code. It is noticable that not all of the relevant contours showed in the left plot are shown on the right one as well. Any suggestion of what am I doing wrong?

data = Import["C:\\Users\\doa\\Documents\\private\\school\\PhD\\thesistex\\sigma_limit_difference.xlsx", {"Sheets", "synthetic clusters"}];
data = data[[2 ;; 3002, 2 ;; 3]];
\[ScriptCapitalD] = SmoothKernelDistribution[data, 9.34];
contour = ContourPlot[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0, 
1024}, PlotRange -> {Automatic}, PlotPoints -> 50, Contours -> 10, ColorFunction -> "DarkRainbow", AxesStyle -> {30, 20}, ImageSize -> Scaled[0.7]];
contour2 =ContourPlot[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0, 
1024}, PlotRange -> {Automatic}, PlotPoints -> 60, Contours -> {2.6*10^-6}, ColorFunction -> "DarkRainbow", AxesStyle -> {30, 20}, ImageSize -> Scaled[0.7],ContourLabels -> True];
{contour, contour2}

I am not sure if there is a need to insert the excel file (in case it is needed please explain how and I will)

thanks in advance

Doron

All thank you very much for the help:

I tried again, this time using some random data so everyone can folow up.

The solution of changing the PlotRange parameter to "All" did'nt solve the problem. It is not just the "White areas" are seen on the left plot with plotrange full, but most impoprtant is whay dont we see them all in the right Plot? Maybe it is related to very small contours?

As for the contour area: Now I get an answer (Thanks...) but still get this warning message:

"NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>"

here is the code I am using:

SeedRandom[5]; data = RandomReal[1024, {3000, 2}];
\[ScriptCapitalD] = SmoothKernelDistribution[data, 9.34];
contour = ContourPlot[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0, 
1024}, PlotRange -> {All}, PlotPoints -> 50, Contours -> 10, ColorFunction -> "DarkRainbow", AxesStyle -> {30, 20}, ImageSize -> Scaled[0.4]];
contour2 =ContourPlot[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y,0,1024}, PlotRange -> {All}, PlotPoints -> 60, Contours -> {2.52*10^-6}, ColorFunction -> "DarkRainbow", AxesStyle -> {30, 20}, ImageSize -> Scaled[0.4], ContourLabels -> True];
{contour, contour2}
GG = NIntegrate[Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024}]

Answer 1

You can play with the method of integration. The default, "GlobalAdaptive", is not very good in this case. You should also consider how much precision is desired. For instance, the example distribution is based on a distribution of 3000 points and the resulting PDF is interpolated on 64 by 64 grid. I'll seek a precision of three digits.

"MultiPeriodic" works without any (helpful) warnings, but perhaps is not as accurate as requested (see further below):

NIntegrate[
  Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024}, 
  PrecisionGoal -> 3, 
  Method -> {"MultiPeriodic", "SymbolicProcessing" -> 0}] // Timing

{0.497470, 291727.}

Increase the requested precision to 4 and we get a warning:

NIntegrate[
  Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024}, 
  PrecisionGoal -> 4, 
  Method -> {"MultiPeriodic", "SymbolicProcessing" -> 0}] // Timing

NIntegrate::maxp: The integral failed to converge after 155092 integrand evaluations.
 NIntegrate obtained 292514.50755723054` and 748.1766581075457` for the integral and
 error estimates. >>

{2.205138, 292515.}

Note that the first answer is not within interval indicated by the second answer and error estimate.

With "LocalAdaptive", one can specify a partitioning of the domain. (No warning messages.)

NIntegrate[
  Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024}, 
  PrecisionGoal -> 3, 
  Method -> {"LocalAdaptive", "Partitioning" -> 128, "SymbolicProcessing" -> 0}] // Timing

{4.386049, 292851.}

Monte Carlo methods can work well in such cases, if high precision is not requried.

NIntegrate[
  Boole[PDF[\[ScriptCapitalD], {x, y}] > 11.52*10^-7], {x, 0, 1024}, {y, 0, 1024}, 
  Method -> {"QuasiMonteCarlo", "Partitioning" -> 64, 
    "MaxPoints" -> 10, "SymbolicProcessing" -> 0}] // Timing

{6.179461, 292700.}

Compare with SEngstrom's answer, which is basically a lower-endpoint Riemann sum:

(data2 = Table[Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1023}, {y, 0, 1023}];
  Select[Flatten[data2], # > 11.52*10^-7 &] // Length) // Timing

{8.264056, 292877}

All these give roughly the same answer (to three digits), the first one seeming a sort of outlier.

---

Another approach to the contour plot, building on Simon Wood's answer, is to use the interpolation grid itself for the plot-points grid. Unfortunately, I have to mask the plot range by hand, for the frame to look nice. :/

With[{ifunc = 
   First @ Cases[PDF[\[ScriptCapitalD], {x, y}], _InterpolatingFunction, -1, Heads -> True]},
 Show[
  ContourPlot[ifunc[x, y], 
   Evaluate[Sequence @@ MapThread[Prepend, {ifunc["Domain"], {x, y}}]], 
   PlotRange -> {{0, 1024}, {0, 1024}}, 
   PlotRangePadding -> Scaled[0.02], 
   PlotPoints -> Most@Dimensions@ifunc["Grid"], MaxRecursion -> 3, 
   Contours -> {2.52*10^-6}],
  Graphics[{White, Rectangle[{1024, -21}, {1045, 1045}], 
    Rectangle[{-21, 1024}, {1045, 1045}], 
    Rectangle[{-21, -21}, {0, 1045}], 
    Rectangle[{-21, -21}, {1045, 0}]}]
  ]
 ]

Answer 2

To force a numerical integration over the discrete Boole[] seems unsuitable to this problem. Why not sample the distribution back to a non-continuous representation and approximate the integral there? You can get fancier with the integration approximation if you need to.

data2 = Table[
   Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0, 1024}];
Plus @@ Select[Flatten[data2], # > 11.52*10^-7 &]

(* 0.441366 *)

As pointed out by Michael E2 the requested property, the area, has a zeroth order approximation of

Select[Flatten[data2], # > 11.52*10^-7 &] // Length

Answer 3

Here's one way to get a contour plot showing all the regions without going to very large numbers of plot points.

The PDF of the SmoothKernelDistribution is just an InterpolatingFunction, so we can examine the underlying data and extract the grid positions at which the function value exceeds the contour value. Those grid positions can then be fed into the contour plot as additional sample points.

pts = With[{ifunc = First @ 
 Cases[PDF[\[ScriptCapitalD], {x, y}], _InterpolatingFunction, -1, Heads -> True]},
   Extract[ifunc["Grid"], 
    Position[ifunc["ValuesOnGrid"], z_ /; z >= 2.52*10^-6]]];

ContourPlot[
 Evaluate@PDF[\[ScriptCapitalD], {x, y}], {x, 0, 1024}, {y, 0, 1024}, 
 PlotRange -> All, PlotPoints -> {50, pts}, Contours -> {2.52*10^-6}]