Sunday, August 22, 2010

Drawing polar plots in three ways using R

A polar plot relates r as a function of angle theta, $$\theta$$. The rectangular components can be
computed can be computed separately as $$x = r cos(theta)$$ and $$y = r sin(theta)$$.

Our solver page at http://extreme.adorio-research.org/solvers/rplotpage/polar/ has three different ways to present polar plots, namely
[polar, rectangular and space]. The range of theta can be specified and the number of points can be set up to 15000, to allow fast drawing without overloading the server.

Here are the three images possible for the butterfly equation, exp(cos(theta)) - 2 * cos(4*theta) + sin(theta/12)^5, presented in Venables and Ripley authoritative book: MASS "Modern Applied Statistics using S" book.








The last image was generated using ntheta = 15000 and pch=5.

Here is the page QP/QPY code for the polar plot page.

# pie.qpy
# 2006.09.02   0.0.1  first version
# 2006.09.20   0.0.2  split from pie,pie, dot page.
# 2010.08.23   0.0.3  third version.

__version__ = "0.0.3"
__date__    = "07.08.12"
__author__  = "eadorio@yahoo.com"
__title__   = "Polar plots using R"
__catalog__ = "POLAR-RPLOT-0127"
__url__     = "/solvers/rplotpage/polar/"
__author__  = "E.P. Adorio"


import time
import tempfile
import commands
import os

from   qp.fill.directory  import Directory
from   qp.fill.form       import Form, StringWidget, TextWidget,CheckboxWidget,SingleSelectWidget
from   qp.fill.css   import BASIC_FORM_CSS

from   qp.sites.extreme.lib.tmpfilesmanager import TmpFilesManager
from   qp.sites.extreme.lib.uicommon        import renderheader, renderfooter, processheader, processfooter
from   qp.pub.common                        import page
from   qp.sites.extreme.lib.checkinput      import checkInputs, getFormStrings
from   qp.sites.extreme.lib.qpyutils        import printRlines, showLogo
from   qp.sites.extreme.lib.webutils        import vecRead,GraphicsFile, as_R_cvector, as_R_vector, as_R_matrix, runRcode
from   qp.sites.extreme.lib import config

       

def Solve(fields):
    (gX, gY, theta0, theta1, ntheta, gtype, psize, equation, col, main, sub, xlab, ylab) = fields

    fname1    = GraphicsFile("png")
    barefile1 = fname1.split(str("/")) [-1]

    Rcode = "png('%s', width=%s*72, height=%s*72)\n" % (fname1, gX, gY)
    Rcode += """
theta  <- seq(%s, %s, len=%s)   
radius <- %s
x      <- radius * cos(theta)
y      <- radius * sin(theta)
""" % (theta0, theta1, ntheta, equation)

    if gtype == "polar":
       Rcode += """
plot(y, x, col="%s",type="l",
     main="%s",sub="%s",xlab="%s",ylab="%s", axes=FALSE)
"""    % (col, main, sub, xlab, ylab)
    elif gtype == "rect":
       Rcode += """
vlo = min(radius, x, y)
vhi = max(radius, x, y)
plot(theta, radius,  ylim=c(vlo, vhi), col="black",type="l", main="%s",sub="%s",xlab="%s",ylab="%s")
points(theta, x,  col="red",type="l")
points(theta, y,  col="blue",type="l")
""" % (main, sub, xlab, ylab)
    elif gtype == "3d":
       Rcode += """
library(scatterplot3d)
scatterplot3d(x, y, theta, highlight.3d=TRUE, col.axis="blue",
              col.grid="lightblue", main="Space Curve", pch=16, cex.symbols = %s)
""" % psize
    Rcode += """
graphics.off()
#img %s
""" % (barefile1)
    status, output = runRcode(Rcode)
    return output
    

class PolarplotPage(Directory):
    def get_exports(self):
        yield ('', 'index', 'PolarPlot', '')
           

    def index[html](self):
        form  = Form(enctype="multipart/form-data")  # enctype for file upload
        form.add(StringWidget,  name="gX",  title = "gX",  value="6", size=3)
        form.add(StringWidget,  name="gY", title = "gY", value="6",  size=3)

        form.add(StringWidget,  name="theta0",  title = "theta0", value = "0", size=5)
        form.add(StringWidget,  name="theta1",  title = "theta1", value = "24*pi", size=5)
        form.add(StringWidget,  name="ntheta",  title = "ntheta", value = "2000", size = 5)
        form.add(SingleSelectWidget,  name="gtype",   title = "Type",
                 value = "polar", options = [("polar", "polar"), ("rect", "rect"),("3d", "space")])
        form.add(StringWidget,  name="psize",  title = "Pt size", value = "0.2", size = 3)

        form.add(StringWidget,  name="equation",  title="Equation in theta r = f(theta)", 
                 value= "exp(cos(theta)) - 2 * cos(4*theta) + sin(theta/12)^5", size = 70)
        form.add(StringWidget,  name="col", title="Color", value="red", size=10)
        form.add(StringWidget,  name="main",  title="Main", value="Polar Plot", size=15)
        form.add(StringWidget,  name="sub",   title="Sub",  value="Example",  size = 15)
        form.add(StringWidget,  name="xlab",  title="xlab", value="X",     size=15)
        form.add(StringWidget,  name="ylab",  title="ylab", value="Y",   size=15)

        form.add_hidden("time",   value = time.time())
        form.add_submit("submit", "submit")

        def render [html] ():
            renderheader(__title__)

            """\
    
%s %s %s %s %s %s %s %s
""" % (form.get_widget("gX").render(), form.get_widget("gY").render(), form.get_widget("theta0").render(), form.get_widget("theta1").render(), form.get_widget("ntheta").render(), form.get_widget("gtype").render(), form.get_widget("col").render(), form.get_widget("psize").render(), ) """\
%s
""" % ( form.get_widget("equation").render(), ) """
%s %s %s %s
""" % (form.get_widget("main").render(), form.get_widget("sub").render(), form.get_widget("xlab").render(), form.get_widget("ylab").render(), ) """ Draws 2D polar plots r vs theta or rectangular plots (r,x,y vs theta) or 3D space curves (x,y,theta) for equation involving theta. The size of the plot is approximately gX inch by gY inch (limited to 10). The point in a space curve will be drawn using filled circles with size determined by Pt size. The following bounds apply to the Open Solver version:
gX float 3, 10
gY float 3, 10
ntheta int 2, 15000
psize float 0.2, 10

Space curves are drawn using the R scatterplot3d package of Uwe Ligges. The butterfly polar function example is from Venables and Ripley, MASS.
""" renderfooter(form, __version__, __catalog__, __author__) if not form.is_submitted(): return page('polarplotpage', render(), style= BASIC_FORM_CSS) def process [html] (): processheader(__title__) calctime_start = time.time() # Get the form input values def tmp1(): return getFormStrings(form, [ "gX", "gY", "theta0", "theta1", "ntheta", "gtype", "psize", "equation", "col", "main", "sub", "xlab", "ylab"]) (gX, gY, theta0, theta1, ntheta, gtype, psize, equation, col, main, sub, xlab, ylab) = tmp1() def tmp2(): inflag = checkInputs( [("gX", gX, ("float", 3, 10)), ("gY", gY, ("float", 3, 10)), ("theta0", theta0, ("eqn", "")), ("theta1", theta1, ("eqn", "")), ("ntheta", ntheta, ("int", 2, 15000)), ("psize", psize, ("float", 0.2, 10)), ("equation", equation, ("eqn",""))]) return inflag inflag = tmp2() if inflag[0]: "

"
               inflag[1]
               "
" else: output = Solve((gX, gY, theta0, theta1, ntheta, gtype, psize, equation, col, main, sub, xlab, ylab)) "
"
               printRlines(output)
               "
" showLogo("Rlogo.jpg") processfooter(form, calctime_start, "./", __url__) process()

Constructive criticisms from our reader are very welcome. Email the page author at ernesto.adorio@gmail.com or alternately at eadorio@yahoo.com

No comments:

Post a Comment