How to make 2D density plot?

Share icon

Tue, 01/22/2013 - 09:14 pm

I have a 3-column matrix attached, and I am trying to make a density plot of it. The first column is the x coordinates, the second column is the y coordinates, and the third column is the corresponding intensity at (x,y). Since my intensity range is pretty large, it is better to make it in log10 scale.

I have also attached an desired graph I made from this data using Mathematica. Could you help me with making this kind of graph using Igor？

Thanks a lot. phyllo-bs.txt bz.png
Hi!

It sounds like you need an imageplot, but that takes a 2D wave (matrix) of Z values, where the X and Y values are either given by the wave scaling, or by separate waves.

If your data is evenly spaced along the X and Y axis, you could reorder the data into a matrix. If it is unevenly spaced, you could interpolate with ImageInterpolate. Which works, but is a tad slow for large data sets.

Alternatively you could do a contour plot, which also takes data in you format of an XYZ-wave triplet.

January 23, 2013 at 05:55 am - Permalink

some general comments:
(1) using a Voronoi interpolation works OK, but to get a fine grid takes a lot of time. Also you will get a spurious level in the "doughnut hole" that you may want to remove (e.g. set to NaN)
(2) A Gizmo scatter plot will work, giving a z-axis view resembling a 2D image. You will have to play with the marker and size.
(3) A Gizmo parametric triangular surface will look best but is complicated to make.
(4) My recommendation for a fast and simple method:

Load the delimited text data as waves wx, wy, wz
Convert to log: wz=log(wz)
Make a graph and display wy vs wx as dots; use the Plan mode for equal x and y scaling.
Using the Modify Trace Appearance dialog, Set as f(z).. color from wz
Modify dot size and colorscale to suit your needs.

Note that your Mathematica display is also doing some form of smoothing or interpolation to achieve its blended appearance.

January 23, 2013 at 06:21 am - Permalink

It appears that your data are sampled with sufficient density at equal radial increments (0.003806). This suggests that with a little bit of effort you can write a function that would produce an approximation using a bilinear-like interpolation value for an arbitrary x and y pair. The following is a quick illustration; it is not complete (need to handle the wrap-around at the end of the cycle) but it does most of what you want. To try it, create a destination wave, e.g.,
•make/n=(100,100) ddd  // use desired resolution here.
•SetScale/I x -1,1,"", ddd
•SetScale/I y -1,1,"", ddd


Paste the following function in your procedure window:
Function myGetZ(inx,iny,wave0,wave1,wave2)
    Variable inx,iny
    Wave wave0,wave1,wave2
   
    Variable rad=sqrt(inx*inx+iny*iny)
    if(rad< 0.284939 || rad>0.95344)
        return NaN
    endif
   
    // find the angle of the interpolated point relative to the origin
    Variable angle=atan2(iny,inx)
   
    // this is the angle of the data relative to the origin
    matrixop/o aa=atan2(wave1,wave0)
    // radial distance
    matrixop/o rr=sqrt(wave0*wave0+wave1*wave1)
    // find the lower angle in the data
    FindLevel/Q aa,(angle)
    Variable startAnglePoint=trunc(V_LevelX)
   
    // finding the next angle can be tricky.  In general it appears that intervals are on the order of
    // 0.015 rad or so.
    FindLevel/Q aa,(angle+0.015)
    Variable nextAnglePoint=trunc(V_LevelX)
    if(nextAnglePoint==startAnglePoint)
        nextAnglePoint+=1
    endif
   
    // add check here if we need to wrap around
   
    // find the point below the required radius
    Variable firstAlpha=aa[startAnglePoint]
    FindLevel/Q/R=[ startAnglePoint,inf] rr,(rad)
    Variable startRadPoint=trunc(V_LevelX)
    Variable firstDataValue=wave2[startRadPoint]
    Variable secondDataValue=wave2[startRadPoint+1]
    Variable dr=rr[startRadPoint+1]-rr[startRadPoint]
    Variable alpha=(rad-rr[startRadPoint])/dr
    Variable firstSideInterpolation=firstDataValue*(1-alpha)+alpha*secondDataValue

    // find the second interpolated value
    FindLevel/Q/R=[ nextAnglePoint,inf] rr,(rad)
    startRadPoint=trunc(V_LevelX)
    firstDataValue=wave2[startRadPoint]
    secondDataValue=wave2[startRadPoint+1]
    dr=rr[startRadPoint+1]-rr[startRadPoint]
    alpha=(rad-rr[startRadPoint])/dr
    Variable secondSideInterpolation=firstDataValue*(1-alpha)+alpha*secondDataValue
    Variable secondAlpha=aa[nextAnglePoint]

    // now interpolate in angle:
    Variable da=secondAlpha-firstAlpha
    if(da==0)
        return firstSideInterpolation
    endif
   
    Variable da2=(angle-firstalpha)/da
    return firstSideInterpolation*(1-da2)+da2*secondSideInterpolation
End


and finally execute:
•ddd=myGetZ(x,y,wave0,wave1,wave2)
•NewImage/K=0 ddd


I hope this helps,

A.G.
WaveMetrics, Inc.

January 23, 2013 at 12:25 pm - Permalink