Contour plot with polar axis
J-E Petit
Thu, 02/05/2015 - 02:44 pm
I am looking for a way to plot a 2-D wave on a polar graph, which doesn't seem to be straight forward on Igor.
I have already gone through the 2 threads dealing with the subject (http://www.igorexchange.com/node/2606 & http://www.igorexchange.com/node/5469), but I must say I am fairly new to programming with Igor (usually use Python) and I can barely find/understand what works best in my case...
Ok, so, I have a 2-D wave, with angles as rows and wind speed as columns (or the other way around, just need to transpose the matrix if needed). Each cells contains concentration values.
How would you (in concrete terms) convert the "rectangle" image to a polar projection?
Many thanks for your help!!
J-E
http://www.igorexchange.com/node/6288
With the resulting polar image, simply add a contour plot of the image.
--Jim Prouty
Software Engineer, WaveMetrics, Inc.
February 5, 2015 at 03:11 pm - Permalink
Your description is somewhat lacking in detail. If I understand you correctly, you have a 2D wave where the rows correspond to some angles and the cols correspond to some speed and the actual cells represent some scalar values. I assume from this that you want to display the scalar values on a coordinate system where the "angles" correspond to say polar angles and the speed would correspond to the magnitude of the polar radius. If that is the case:
Wave inWave
Variable rows=dimsize(inWave,0)
Variable cols=dimSize(inWave,1)
Variable points=rows*cols
Make/O/D/N=(points,3) polarTriplet
Variable i,j,angle,count=0,rad
// assuming angle is in radians measured CC from the positive x-axis:
for(i=0;i<cols;i+=1)
rad=dimOffset(inWave,1)+i*dimDelta(inWave,1)
for(j=0;j<rows;j+=1)
angle=dimOffset(inWave,0)+j*dimDelta(inWave,0)
polarTriplet[count][0]=rad*cos(angle)
polarTriplet[count][1]=rad*sin(angle)
polarTriplet[count][2]=inWave[j][i]
count+=1
endfor
endfor
End
When you have converted to polarTriplet simply display a contour of this triplet wave.
HTH,
A.G.
WaveMetrics, Inc.
February 5, 2015 at 03:21 pm - Permalink
Thanks veru much for your answers!!! I tried both suggestions, but before, I'd like to be more thorough on the description of the data and the objective of my procedure :)
First, I enclose the polar plot I get from my Python script, which I would like to reproduce on Igor (file Polar plot Python). the polar angles are obviously associated with wind direction, and radius axis with wind speed. f(z) colors are linked to concentration values at each wind direction & speed.
Then, I enclose what I have in Igor (already reproduce the calculation to compute mx values). Image "Image_mx". My 2D wave, called "mx", is composed of 720 rows, describing angles from 0° to 359.5° with an angular resolution of 0.5°, and 40 columns representing wind speed from 0 to 20 km/h and a speed resolution of 0.5 km/h.
The idea is then to "tranform" the rectangular image plot to polar, just like the graph from Python.
I tried the convertToPolar function, but the result doesn't look like satisfactory. (image "Contour_plt_convertToPolar"). I used 100 levels, without labels.
I tried the DisplayPolarImage macro, and got a polar-like graph, which is a good news! But the data used doesn't look like correct. I executed "DisplayPolarImage("mx",2,angle_rad)"; where angle-rad is a 1D wave containing 720 rows describing the angles in radian.
Maybe I didn't use your functions correctly?
Thanks again for your help
J-E
February 6, 2015 at 02:47 am - Permalink
I am not sure exactly what you tried there. Remember, posting images is not nearly as useful as posting a small IGOR experiment with the relevant data.
Before you use convertToPolar() you need to make sure that your wave scaling is indeed in radians. I suspect it might be in degrees which could explain the strange result. After you change either the wave scaling or add the appropriate factor to the "angle" variable you should get a polarTriplet wave which you can now interpolate at any resolution you want. The pixelated image you show in one of your results indicates insufficient sampling. If you use ImageInterpolate (with the keyword Voronoi) directly, you have the freedom to adjust the parameters to the /S flag and determine the boundaries and sampling in both directions. Your command would look something like:
The results are stored in the wave M_InterpolatedImage in the current data folder.
HTH,
A.G.
WaveMetrics, Inc.
February 6, 2015 at 09:45 am - Permalink
wave scaling is a concept which I'm not fully confident with. But before bothering you, I'll go through the manual; I'll probably find all I need :)
February 6, 2015 at 11:44 am - Permalink
I think I got my wave scaling right; but I realized that your function assumes anles (in rad) calculated counterclockwise from the positive x-axis. And the polar graph would need clockwise angles from the positive y axis (0° corresponds to the North). This would probably explain the weird countour plot that I get.
I guess the "most simple" way to overcome this is to modify the converttopolar() procedure; but I would be unable to do it myself...
Any hint?
If necessary, I can provide the 2D wave I'm using(?!)
Thanks!
February 9, 2015 at 07:04 am - Permalink
I'm not a big fan of wave scaling myself so it is fine by me if you don't use it. In that case you need to assume something about your sampling, and state it. Here is how I'd change the code:
Wave inWave
Variable rows=dimsize(inWave,0)
Variable cols=dimSize(inWave,1)
Variable points=rows*cols
Make/O/D/N=(points,3) polarTriplet
Variable i,j,angle,count=0,rad
// assuming angle is in degrees measured clockwise from the y-axis
Variable da=pi/360 // 2*pi/720
for(i=0;i<cols;i+=1)
rad=0.5*i // wind speed in 0.5 km/hr increments
for(j=0;j<rows;j+=1)
angle=j*da //
polarTriplet[count][0]=rad*sin(angle)
polarTriplet[count][1]=rad*cos(angle)
polarTriplet[count][2]=inWave[j][i]
count+=1
endfor
endfor
End
HTH,
AG
February 9, 2015 at 09:49 am - Permalink
In a similar way, I created 3 waves so that for each angle, there is 200 corresponding wind speed values; and for each (angle,speed) couple there is one concentration value. It is just a matter of wave redimensionning (code is below)
wave inWave
Variable rows=dimsize(inWave,0)
Variable cols=dimSize(inWave,1)
Variable points=rows*cols
Make/O/D/N=(points) xpolar
Make/O/D/N=(points) ypolar
wave mx
duplicate mx mxT
MatrixTranspose mxT
variable i,count=0,j=0
wave angle, speed
for (i=cols-1;i<points;i+=cols)
xpolar[i-cols+1,i]=angle[count]
count+=1
endfor
for (i=0;i<rows;i+=1)
Concatenate/O/NP {speed, ypolar}, temp
ypolar=temp
endfor
KillWaves temp
duplicate mxT zpolar
Redimension/E=1/N=(points) zpolar
End Function
Function plot()
wave ypolar, xpolar, zpolar
WMNewPolarGraph("", "")
WMPolarAppendTrace("PolarGraph0",ypolar, xpolar,360)
ModifyGraph mode=2,lsize(polarY0)=2,zColor(polarY0)={zpolar,*,*,Rainbow,1}
SetDataFolder root:Packages:WMPolarGraphs:PolarGraph0:
variable/G angleDirection
variable/G zeroAngleWhere
angleDirection=1
zeroAngleWhere=90
End Function
Then, I have a nice surface plot on polar axes, but the axes are all below the trace (see enclosed image). I checked the "draw on Top of Graph" button for the Vert and HorizCrossing axes, but it didn't do anything. Is there an other way to overcome this?
Thanks!
February 16, 2015 at 01:03 am - Permalink
Please do; I think that using the Polar Graphing package for your data is not a good idea (the polar graphing macros don't provide for drawing the axes in front of the traces).
Creating an image plot and your own polar axes is more workable.
--Jim Prouty
Software Engineer, WaveMetrics, Inc.
February 16, 2015 at 11:54 am - Permalink
"mx" consists of 720 rows (corresponding to angles from 0 to 360 with a resolution of 0.5) and 200 columns (corresponding to wind speed from 0 to 20 km/h with a resolution of 0.1 km/h)
I also enclosed the "angle" and "speed" waves, just in case.
February 17, 2015 at 01:18 am - Permalink
SetScale/P y 0,0.1,"kmph", mx // 100 meters per hour per column
Then I opened Polar Image.ipf (saved from http://www.igorexchange.com/node/6288) as a Procedure file, and then closed it (without killing it), which automatically compiles all procedures.
Then I chose DisplayPolarImage from the Macros menu, which displays first one dialog where you choose the matrix, choose degrees for angle units, and choose X scaling for Angle from.
Click Continue to get to the second dialog, use conservative numbers like 200x200 rows x columns for the rectangular polar image because this part is dog-slow. It took about an hour to complete the Voronoi interpolation.
If you chose Yes for Display matrix as image? then you see a grayscale image plot with reversed Y axis. Ignore this, we can do better.
There is another way to display the image by using the attached experiment which modifies the New Polar Graphs code by having these Override definitions in the main Procedure window (the Overrides cannot be moved elsewhere):
// Overrides work from only the main Procedure window.
Override StrConstant sLayerGrid = "ProgFront"
Override StrConstant sLayerAxes = "ProgFront"
Override StrConstant sLayerAxisLabels = "UserFront"
Macro MakePolarGraphForImage(matrixName)
String matrixName
Prompt matrixName, "polar Image", popup, WaveList("*Image",";","DIMS:2")
fMakePolarGraphForImage($matrixName)
End
Function fMakePolarGraphForImage(matrix)
Wave matrix
// create a polar pair that spans the matrix extent
Variable dx= DimDelta(matrix,0)
Variable xmin= DimOffset(matrix,0)
Variable xmax= xmin + dx*(DimSize(matrix,0)-1)
Variable dy= DimDelta(matrix,1)
Variable ymin= DimOffset(matrix,1)
Variable ymax= ymin + dy*(DimSize(matrix,1)-1)
String matName= NameOfWave(matrix)
String radiusname= CleanupName(matName[0,28]+"_r",1)
String angleName= CleanupName(matName[0,28]+"_a",1) // radians
Make/O/D/N=0 $radiusname/WAVE=radius
Make/O/D/N=0 $angleName/WAVE=radians
Variable/C topLeft= r2polar(cmplx(xmin,ymin))
Variable/C bottomRight= r2polar(cmplx(xmax,ymax))
radius= {real(topLeft), real(bottomRight)}
radians= {imag(topLeft), imag(bottomRight)}
String graphName= UniqueName("PolarImage",6,0)
graphName= WMNewPolarGraph("_default_", graphName)
String traceName= WMPolarAppendTrace(graphName,radius, radians, 2*pi)
ModifyGraph/W=$graphName lsize($traceName)=0
AppendImage/W=$graphName/L=VertCrossing/B=HorizCrossing matrix
ModifyImage/W=$graphName $NameOfWave(matrix) ctab= {*,*,Rainbow256,0}
WMPolarAxesRedrawGraphNow(graphName)
WMPolarGraphs(1) // show the panel, switch to Range tab
WMPolarSwitchTab(1)
End
Then call the MakePolarGraphForImage macro with mxImage to get polar axes on top of the polar image.
You can set the range in the Polar Graphs Range tab and the Angle Axes Thickness to mask the ragged edge a bit. See the resulting image and attached experiment.
--Jim Prouty
Software Engineer, WaveMetrics, Inc.
February 17, 2015 at 01:12 pm - Permalink
February 23, 2015 at 02:47 am - Permalink
That was motivated by changes I made to the MakePolarGraphForImage procedure which uses Polar Image.ipf.
This revision still requires you to set the resulting polar graph's angle rotation direction and angle=0 location to match the values given to DisplayPolarImage().
The revised procedures prompt the user for those parameters. Really the only change to this code is:
For convenience, here is the entire revised Procedure code (this code MUST be in the main Procedure window):
#include <New Polar Graphs>
// Overrides work from only the main Procedure window.
Override StrConstant sLayerGrid = "ProgFront"
Override StrConstant sLayerAxes = "ProgFront"
Override StrConstant sLayerAxisLabels = "UserFront"
Macro MakePolarGraphForImage(matrixName)
String matrixName
Prompt matrixName, "polar Image", popup, WaveList("*Image*",";","DIMS:2")
fMakePolarGraphForImage($matrixName)
End
Function fMakePolarGraphForImage(matrix)
Wave matrix
// create a polar pair that spans the matrix extent
Variable dx= DimDelta(matrix,0)
Variable xmin= DimOffset(matrix,0)
Variable xmax= xmin + dx*(DimSize(matrix,0)-1)
Variable dy= DimDelta(matrix,1)
Variable ymin= DimOffset(matrix,1)
Variable ymax= ymin + dy*(DimSize(matrix,1)-1)
String matName= NameOfWave(matrix)
String radiusname= CleanupName(matName[0,28]+"_r",1)
String angleName= CleanupName(matName[0,28]+"_a",1) // radians
Make/O/D/N=0 $radiusname/WAVE=radius
Make/O/D/N=0 $angleName/WAVE=radians
Variable/C topLeft= r2polar(cmplx(xmin,ymin))
Variable/C bottomRight= r2polar(cmplx(xmax,ymax))
radius= {real(topLeft), real(bottomRight)}
radians= {imag(topLeft), imag(bottomRight)}
String graphName= UniqueName("PolarImage",6,0)
graphName= WMNewPolarGraph("_default_", graphName)
String traceName= WMPolarAppendTrace(graphName,radius, radians, 2*pi)
ModifyGraph/W=$graphName lsize($traceName)=0
AppendImage/W=$graphName/L=VertCrossing/B=HorizCrossing matrix
ModifyImage/W=$graphName $NameOfWave(matrix) ctab= {*,*,Rainbow256,0}
WMPolarAxesRedrawGraphNow(graphName)
WMPolarGraphs(1) // show the panel, switch to Range tab
WMPolarSwitchTab(1)
End
--Jim Prouty
Software Engineer, WaveMetrics, Inc.
March 11, 2015 at 04:23 pm - Permalink
I took a little bit of time to get back to your procedure. It works nicely, and what is useful is that you can change the settings of the graph through the polargraph panel. However, the calculation of mxImage take way too much time...
So as I tried earlier in this post, I plotted as a surface plot my triplet wave on a polar graph. And as the axes were beneath the surface plot, I tried to manually draw them on the graph. I don't really know if it is the most suitable solution, because one needs to interact with the "axes" through the procedure.... but it's works pretty well and it's fast....
wave inWave
Variable rows=dimsize(inWave,0)
Variable cols=dimSize(inWave,1)
Variable points=rows*cols
Make/O/D/N=(points) xpolar
Make/O/D/N=(points) ypolar
variable i,count=0,j=0
wave angle, speed
for (i=cols-1;i<points;i+=cols)
xpolar[i-cols+1,i]=angle[count]
count+=1
endfor
for (i=0;i<rows;i+=1)
Concatenate/O/NP {speed, ypolar}, temp
ypolar=temp
endfor
KillWaves temp
duplicate inWave zpolar
MatrixTranspose zpolar
Redimension/E=1/N=(points) zpolar
End Function
Function plot(MxToPlot)
wave MxToPlot
polar(MxToPlot)
wave ypolar, xpolar, zpolar
Make/O/N=3600 edge_angle, edge_line0, edge_line1, edge_line2,edge_line3
edge_angle=p*0.1
edge_line0=Wavemax(ypolar)
edge_line1=edge_line0/4
edge_line2=edge_line1*2
edge_line3=edge_line1*3
Make/O/N=2 x_axis1, x_axis2, x_axis3, x_axis4, y_axis
y_axis=Wavemax(ypolar)
x_axis1=180*p
x_axis2=x_axis1+45
x_axis3=x_axis2+45
x_axis4=x_axis3+45
WMNewPolarGraph("", "")
SetDataFolder root:Packages:WMPolarGraphs:PolarGraph0:
variable/G angleDirection
variable/G zeroAngleWhere
variable/G tickLabelBlue
variable/G tickLabelGreen
variable/G tickLabelRed
variable/G radiusAxisColorBlue
variable/G radiusAxisColorGreen
variable/G radiusAxisColorRed
variable/G angleAxisColorBlue
variable/G angleAxisColorGreen
variable/G angleAxisColorRed
angleDirection=1
zeroAngleWhere=90
tickLabelBlue=65535
tickLabelGreen=65535
tickLabelRed=65535
radiusAxisColorBlue=65535
radiusAxisColorGreen=65535
radiusAxisColorRed=65535
angleAxisColorBlue=65535
angleAxisColorGreen=65535
angleAxisColorRed=65535
SetDataFolder root:
WMPolarAppendTrace("PolarGraph0",ypolar, xpolar,360)
WMPolarAppendTrace("PolarGraph0",edge_line0, edge_angle,360)
WMPolarAppendTrace("PolarGraph0",edge_line1, edge_angle,360)
WMPolarAppendTrace("PolarGraph0",edge_line2, edge_angle,360)
WMPolarAppendTrace("PolarGraph0",edge_line3, edge_angle,360)
WMPolarAppendTrace("PolarGraph0",y_axis, x_axis1,360)
WMPolarAppendTrace("PolarGraph0",y_axis, x_axis2,360)
WMPolarAppendTrace("PolarGraph0",y_axis, x_axis3,360)
WMPolarAppendTrace("PolarGraph0",y_axis, x_axis4,360)
ModifyGraph width=566.929,height=566.929
ModifyGraph mode=2,lsize(polarY0)=5,zColor(polarY0)={zpolar,*,*,SeaLandAndFire,0}
ModifyGraph mode(polarY1)=0,lsize(polarY1)=3,rgb(polarY1)=(0,0,0)
ModifyGraph lstyle(polarY2)=2,rgb(polarY2)=(65535,65535,65535),mode(polarY2)=0
ModifyGraph lstyle(polarY3)=2,rgb(polarY3)=(65535,65535,65535),mode(polarY3)=0
ModifyGraph lstyle(polarY4)=2,rgb(polarY4)=(65535,65535,65535),mode(polarY4)=0
ModifyGraph lstyle(polarY5)=2,rgb(polarY5)=(65535,65535,65535),mode(polarY5)=0
ModifyGraph lstyle(polarY6)=2,rgb(polarY6)=(65535,65535,65535),mode(polarY6)=0
ModifyGraph lstyle(polarY7)=2,rgb(polarY7)=(65535,65535,65535),mode(polarY7)=0
ModifyGraph lstyle(polarY8)=2,rgb(polarY8)=(65535,65535,65535),mode(polarY8)=0
TextBox/C/B=1/F=0/A=MC/N=text0/X=0.18/Y=49.28 "\\Z18\\f01N"
TextBox/C/B=1/F=0/A=MC/N=text1/X=36.27/Y=34.91 "\\Z18\\f01NE"
TextBox/C/B=1/F=0/A=MC/N=text2/X=48.24/Y=0.46 "\\Z18\\f01E"
TextBox/C/B=1/F=0/A=MC/N=text3/X=37.44/Y=-34.87 "\\Z18\\f01SE"
TextBox/C/B=1/F=0/A=MC/N=text4/X=0.95/Y=-49.39 "\\Z18\\f01S"
TextBox/C/B=1/F=0/A=MC/N=text5/X=-36.28/Y=-36.36 "\\Z18\\f01SW"
TextBox/C/B=1/F=0/A=MC/N=text6/X=-47.49/Y=0.27 "\\Z18\\f01W"
TextBox/C/B=1/F=0/A=MC/N=text7/X=-34.51/Y=37.27 "\\Z18\\f01NW"
Make/O/N=3 radialLabel
variable/G s_res
radialLabel=(p+1)*(WaveMax(ypolar)+s_res)/4
TextBox/C/N=text8/F=0/X=37.26/Y=47.88 num2str(radialLabel[0])
TextBox/C/N=text9/F=0/X=25.58/Y=47.88 num2str(radialLabel[1])
TextBox/C/N=text10/F=0/X=14.99/Y=47.88 num2str(radialLabel[2])
End Function
March 18, 2015 at 03:16 pm - Permalink