Rose diagram code

Today I show you how to draw geological rose diagram using Python. All you need to prepare is to download and install Python, Numpy and Matplotlib.

In this code, we have 5 files listed as following:

NorthPolarAxes.py
azimuth.py
dip.py
azimuth_dip.py
configure

These two pictures are azimuth rose diagram and dip rose diagram separately.

NorthPolarAxes.py is the file for drawing rose-diagram-like polar coordinates;
azimuth.py takes charge of drawing rose diagram according to dip direction;
dip.py is responsible for plotting rose diagram base on dip;
azimuth_dip.py simply draws both plots mentioned above together;
configure is the file where you input your parameters for plotting.

##################NorthPolarAxes.py
import numpy as np
from matplotlib.projections import PolarAxes
from matplotlib.transforms import Affine2D, Bbox, IdentityTransform

class NorthPolarAxes(PolarAxes):

name = ‘northpolar’

def format_coord(self, angle, radius):
angle = np.pi * 0.5 – angle
if angle < 0:
angle = angle + 2 * np.pi
angle = angle * 180 / np.pi
return u’\u03b8=%f\u00b0, r=%f’ % (angle, radius)

class NorthPolarTransform(PolarAxes.PolarTransform):
def transform(self, tr):
xy = np.zeros(tr.shape, np.float_)
t = tr[:, 0:1]
r = tr[:, 1:2]
x = xy[:, 0:1]
y = xy[:, 1:2]
x[:] = r * np.sin(t)
y[:] = r * np.cos(t)
return xy

transform_non_affine = transform

def inverted(self):
return NorthPolarAxes.InvertedNorthPolarTransform()

class InvertedNorthPolarTransform(PolarAxes.InvertedPolarTransform):
def transform(self, xy):
x = xy[:, 0:1]
y = xy[:, 1:]
r = np.sqrt(x*x + y*y)
theta = np.arctan2(y, x)
return np.concatenate((theta, r), 1)

def inverted(self):
return NorthPolarAxes.NorthPolarTransform()

def _set_lim_and_transforms(self):
PolarAxes._set_lim_and_transforms(self)
self.transProjection = self.NorthPolarTransform()
self.transData = (
self.transScale +
self.transProjection +
(self.transProjectionAffine + self.transAxes))
self._xaxis_transform = (
self.transProjection +
self.PolarAffine(IdentityTransform(), Bbox.unit()) +
self.transAxes)
self._xaxis_text1_transform = (
self._theta_label1_position +
self._xaxis_transform)
self._yaxis_transform = (
Affine2D().scale(np.pi * 2.0, 1.0) +
self.transData)
self._yaxis_text1_transform = (
self._r_label1_position +
Affine2D().scale(1.0 / 360.0, 1.0) +
self._yaxis_transform)
##################

##################azimuth.py
import numpy as np
import matplotlib.pyplot as pl
from matplotlib.projections import register_projection
from NorthPolarAxes import NorthPolarAxes
import ConfigParser

config = ConfigParser.ConfigParser()
config.read(‘configure’)
filename = config.get(‘configure’, ‘filename’)
angle = int(config.get(‘configure’, ‘angle’))

register_projection(NorthPolarAxes)

data = np.genfromtxt(filename)

nsection = 360 / angle
direction = np.linspace(0, 360, nsection, False) / 180 * np.pi
print direction
frequency = [0] * (nsection)
for i in range(len(data)):
tmp = int((data[i][1] – data[i][1] % angle) / angle)
frequency[tmp] = frequency[tmp] + 1
width = angle / 180.0 * np.pi * np.ones(nsection)

ax = pl.subplot(1, 1, 1, projection = ‘northpolar’)
bars = ax.bar(direction, frequency, width=width, bottom=0.0)
for r,bar in zip(frequency, bars):
bar.set_facecolor(pl.cm.jet(0.8))
bar.set_edgecolor(‘grey’)
bar.set_alpha(0.8)

pl.show()
##################

##################dip.py
import numpy as np
import matplotlib.pyplot as pl
from matplotlib.projections import register_projection
from NorthPolarAxes import NorthPolarAxes
import ConfigParser

config = ConfigParser.ConfigParser()
config.read(‘configure’)
filename = config.get(‘configure’, ‘filename’)
angle = int(config.get(‘configure’, ‘angle’))

register_projection(NorthPolarAxes)

data = np.genfromtxt(filename)

nsection = 360 / angle
direction = np.linspace(0, 360, nsection, False) / 180 * np.pi
frequency = [0] * (nsection)
dipangle = [0] * (nsection)
for i in range(len(data)):
tmp = int((data[i][1] – data[i][1] % angle) / angle)
frequency[tmp] = frequency[tmp] + 1
dipangle[tmp] = dipangle[tmp] + data[i][0]
for i in range(nsection):
if frequency[i] > 0:
dipangle[i] = dipangle[i] / frequency[i]

width = angle / 180.0 * np.pi * np.ones(nsection)

ax = pl.subplot(1, 1, 1, projection = ‘northpolar’)
bars = ax.bar(direction, dipangle, width=width, bottom=0.0)
for r,bar in zip(frequency, bars):
bar.set_facecolor(pl.cm.jet(0.8))
bar.set_edgecolor(‘grey’)
bar.set_alpha(0.8)

pl.show()
##################

##################azimuth_dip.py
import numpy as np
import matplotlib.pyplot as pl
from matplotlib.projections import register_projection
from NorthPolarAxes import NorthPolarAxes
import ConfigParser

config = ConfigParser.ConfigParser()
config.read(‘configure’)
filename = config.get(‘configure’, ‘filename’)
angle = int(config.get(‘configure’, ‘angle’))

register_projection(NorthPolarAxes)

data = np.genfromtxt(filename)

nsection = 360 / angle
direction = np.linspace(0, 360, nsection, False) / 180 * np.pi
frequency = [0] * (nsection)
dipangle = [0] * (nsection)
for i in range(len(data)):
tmp = int((data[i][1] – data[i][1] % angle) / angle)
frequency[tmp] = frequency[tmp] + 1
dipangle[tmp] = dipangle[tmp] + data[i][0]
for i in range(nsection):
if frequency[i] > 0:
dipangle[i] = dipangle[i] / frequency[i]

width = angle / 180.0 * np.pi * np.ones(nsection)

ax = pl.subplot(1, 2, 1, projection = ‘northpolar’)
bars = ax.bar(direction, frequency, width=width, bottom=0.0)
for r,bar in zip(frequency, bars):
bar.set_facecolor(pl.cm.jet(0.8))
bar.set_edgecolor(‘grey’)
bar.set_alpha(0.8)

ax = pl.subplot(1, 2, 2, projection = ‘northpolar’)
bars = ax.bar(direction, dipangle, width=width, bottom=0.0)
for r,bar in zip(frequency, bars):
bar.set_facecolor(pl.cm.jet(0.8))
bar.set_edgecolor(‘grey’)
bar.set_alpha(0.8)

pl.show()
##################

##################configure
[configure]
# filename is the name of the file which contains dip and dip direction
# file contains two columns: 1st column is dip; 2nd column dip direction
filename = dip_azim
# angle is a degree, based on which a whole circle is going to be divided into 360 / angle parts
angle = 5
##################

You can download the source and example and test it by yourself.

This entry was posted in code and tagged , , , , , . Bookmark the permalink.

13 Responses to Rose diagram code

  1. yonggeng says:

    I hope it’s not difficult to use. If you have any questions about it, please leave me a message.

  2. mikepsn says:

    Hi,
    I used your north polar projection code and it works great except for one minor problem.
    When I use pyplot to show the plot and then hover the mouse over the plot to get the coordinates I get error messages saying that the global name “InvertedNorthPolarTransform” is not defined.

    File “C:\Python27\lib\lib-tk\Tkinter.py”, line 1410, in __call__
    return self.func(*args)
    File “C:\Python27\lib\site-packages\matplotlib\backends\backend_tkagg.py”, lin
    e 279, in motion_notify_event
    FigureCanvasBase.motion_notify_event(self, x, y, guiEvent=event)
    File “C:\Python27\lib\site-packages\matplotlib\backend_bases.py”, line 1624, i
    n motion_notify_event
    guiEvent=guiEvent)
    File “C:\Python27\lib\site-packages\matplotlib\backend_bases.py”, line 1262, i
    n __init__
    LocationEvent.__init__(self, name, canvas, x, y, guiEvent=guiEvent)
    File “C:\Python27\lib\site-packages\matplotlib\backend_bases.py”, line 1182, i
    n __init__
    xdata, ydata = self.inaxes.transData.inverted().transform_point((x, y))
    File “C:\Python27\lib\site-packages\matplotlib\transforms.py”, line 1953, in i
    nverted
    return CompositeGenericTransform(self._b.inverted(), self._a.inverted())
    File “C:\Python27\lib\site-packages\matplotlib\transforms.py”, line 1953, in i
    nverted
    return CompositeGenericTransform(self._b.inverted(), self._a.inverted())
    File “C:\Documents and Settings\papasimm\Desktop\rcs\northpolar.py”, line 28,
    in inverted
    return InvertedNorthPolarTransform()
    NameError: global name ‘InvertedNorthPolarTransform’ is not defined

    • yonggeng says:

      I think you forgot the overall class name when you used InvertedNorthPolarTransform.
      In my file it is not a global name. You can check line 30 in file NorthPolarAxes.py.

  3. lrs says:

    Hi, thank you for your piece of code. I wanted to use your north polar projection and since I didn’t get it running, I tried your example azimuth_dip.py. Also here I receive the following error message:

    File “NorthPolarAxes.py”, line 46, in _set_lim_and_transforms
    self.transProjection = self.NorthPolarTransform()
    TypeError: __init__() takes exactly 2 arguments (1 given)

    Do you have any idea why this happens?

    • yonggeng says:

      Hello Irs, I tested with Python 2.6 and 2.7. I works fine.
      Which version do you use? or did you change anything?

      • lrs says:

        Hi, thanks for your quick reply. I use Python 2.6 and the only thing which I had to change is the reading of the data in azimuth_dip.py because the input with config also didn’t work properly. But this shouldn’t affect the problem stated above in NorthPolarAxes.py as far as I can see. I’ll try to install Python 2.7, maybe it’ll work then.

  4. yonggeng says:

    Good luck. If it still doesn’t work, you can send me your code and let me take a look if you want.

    • lrs says:

      Which matplotlib version are you using? I use 0.98.5.2 on an ubuntu Linux machine. I didn’t get the code running, I’ll try now to write an own northpolar version starting from the polar.py on my machine and using yours as guideline…

    • lrs says:

      It works now! I guess the problem was due to different matplotlib versions!

  5. dinesh says:

    Hi,

    I am trying to download the sourse code, but the weblink is broken. Could let me know from where i can down load the sourse code of the rose diagram.

  6. Nate T. says:

    I’m grateful you’re sharing this code, but it could use some comments for us mere mortals to understand it better.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s