Search code examples
pythonmatplotlibpolar-coordinates

Stereographic Sun Diagram matplotlib polar plot python


I am trying to create a simple stereographic sun path diagram similar to these: http://wiki.naturalfrequency.com/wiki/Sun-Path_Diagram

I am able to rotate a polar plot and set the scale to 90. How do I go about reversing the y-axis? Currently the axis goes from 0>90, how do I reverse the axis to 90>0 to represent the azimuth?

I have tried:

ax.invert_yaxis()
ax.yaxis_inverted()

Further, how would I go about creating a stereographic projection as opposed to a equidistant?

My code:

import matplotlib.pylab as plt
testFig = plt.figure(1, figsize=(8,8))
rect = [0.1,0.1,0.8,0.8]
testAx = testFig.add_axes(rect,polar=True)
testAx.invert_yaxis()
testAx.set_theta_zero_location('N')
testAx.set_theta_direction(-1)

Azi = [90,180,270]
Alt= [0,42,0]
testAx.plot(Azi,Alt)
plt.show()

Currently my code doesn't seem to even plot the lines correctly, do I need need to convert the angle or degrees into something else?

Any help is greatly appreciated.


Solution

  • I finally had time to play around with matplotlib. After much searching, the correct way as Joe Kington points out is to subclass the Axes. I found a much quicker way utilising the excellent basemap module.

    Below is some code I have adapted for stackoverflow. The sun altitude and azimuth were calculated with Pysolar with a set of timeseries stamps created in pandas.

    import matplotlib.pylab as plt
    from mpl_toolkits.basemap import Basemap
    import numpy as np
    
    winterAzi = datafomPySolarAzi
    winterAlt = datafromPySolarAlt
    
    # create instance of basemap, note we want a south polar projection to 90 = E
    myMap = Basemap(projection='spstere',boundinglat=0,lon_0=180,resolution='l',round=True,suppress_ticks=True)
    # set the grid up
    gridX,gridY = 10.0,15.0
    parallelGrid = np.arange(-90.0,90.0,gridX)
    meridianGrid = np.arange(-180.0,180.0,gridY)
    
    # draw parallel and meridian grid, not labels are off. We have to manually create these.
    myMap.drawparallels(parallelGrid,labels=[False,False,False,False])
    myMap.drawmeridians(meridianGrid,labels=[False,False,False,False],labelstyle='+/-',fmt='%i')
    
    # we have to send our values through basemap to convert coordinates, note -winterAlt
    winterX,winterY = myMap(winterAzi,-winterAlt)
    
    # plot azimuth labels, with a North label.
    ax = plt.gca()
    ax.text(0.5,1.025,'N',transform=ax.transAxes,horizontalalignment='center',verticalalignment='bottom',size=25)
    for para in np.arange(gridY,360,gridY):
        x= (1.1*0.5*np.sin(np.deg2rad(para)))+0.5
        y= (1.1*0.5*np.cos(np.deg2rad(para)))+0.5
        ax.text(x,y,u'%i\N{DEGREE SIGN}'%para,transform=ax.transAxes,horizontalalignment='center',verticalalignment='center')
    
    
    # plot the winter values
    myMap.plot(winterX,winterY ,'bo')
    

    Note that currently I am only plotting points, you will have to make sure that line points have a point at alt 0 at sunrise/sunset.

    stereographic plot, note I have plotted Winter/Summer Solstice and Autumn Equinox