Well basically title is what I'm trying to do except I would also like to

sum up (integrate) the projected values in bins.

In other words it could be described as a 2D heatmap of a 3D point set.

In another other words, imagine a normal distribution projected from the z

axis to the xy plane. What you see is basically a 2D normal distribution.

Moving an azimuthal angle theta (like in spherical coordinates) away from

the z axis when you project the same normal distribution what you would get

is not the "normal" normal dist but a skewed normal distribution.

I thought I can get the same results by rotating the normal distribution,

projecting it on xy plane and then doing a heatmap, hexbin, histogram or

something similar. Here's the rotation an projection via the contourf

method, but I can't reproduce the projections myself.

I create an x, y and z as 2D arrays (xy meshgird of coordinates, and z as a

Cartesian product of 1D normal dist.). I rotate them by multiplying the

rotation matrix with the column vector of every point on the grid just like

in basic lin alg. No issues.

Rotation matrix is created in the rotation_matrix_weave but for the people

that have too modern scipy's/mpl's the rotation_matrix_numpy does the same

thing it just takes more time.

The points are rotated in the RotateFunc (where the meshgrid and Cartesian

product is done too).

I plot the rotated normal dist. and use the contourf functions to show the

projections on xy and yz planes. The kind of projection that is visible in

the xy plane is the one I'd like to get. In all fairness because these are

projections of points if you would just rotate an empty xy plane what you

would get as a projection back on to the xy plane is a gradient. So I

rotate an "empty" plane and subtract it from the rotated normal projection

to get rid of the gradient and just keep the skewed gaussian.

But this is where all the problems start. I can't figure out how to get

that projection. I tried finding where the countourf functions are defined,

so I could try to figure out how it does it but couldn't make out the code

(it's fairly big) and it doesn't really do what I need it to.

Easier approach is to just construct the projections over pcolor or

pcolormesh (much faster) as I've added at the bottom of the script. But the

online wisdom on how to reconstruct the arrays won't work because

pcolormesh/pcolor understands that "indices aren't coordinates". That is to

say, pcolormesh doesn't simply just plot rz, but it plots the value rz[i,j]

at the coordinates ( rx[i,j], ry[i,j] ) in this new 2D array that is the

final image.

I've been fishing for help online and in IRC for 2-3 days now and I've been

lost for over a week, by this point honestly, whatever works. I've resorted

to saving the graph as an image and reading that in with OpenCv. But the

problem is they are all "aspected" so from the original (N, N) array you

end up with a (1000, 800) image and I really don't like the possibility

that the original data can change like that. I mean contourf, pcolor and

pcolormesh can do it, there's got to be a way to extract the 2D array or at

least some kind of bin-value pairs or something from them.

So to reiterate:

1) can someone show me how to get the end graph shown by pcolormesh as an

array?

2) if that's not possible is pcolor any different?

3) if it's not, does anyone have insight on how pcolormesh/pcolor can end

up drawing the graph correctly and can I get pointers on how to do what

they do myself?

4) is there perhaps some other method/function I'm not looking at (i.e

hexbin) that can be used to get the result?

5) if that's not possible at all, could it be possible to change from the

current representation to i.e. [ [x, y, z]....] of the points and use that

to get a histogram/hexbin? Those are the xyz coordinates of points, that is

[ [ rx[0,0] , ry[0,0] , rz[0,0] ], ...... [ rx[i,j] , ry[i,j] , rz[i,j]

], ..... [ rx[-1,-1] , ry[-1,-1] , rz[-1,-1] ] ].

6) am I stupid?

Sorry for the long mail, thanks for all the help.

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from scipy import stats

from matplotlib import cm

from math import cos, sin, sqrt

from scipy import weave

def rotation_matrix_numpy(axis, theta):

mat = np.eye(3,3)

axis = axis/sqrt(np.dot(axis, axis))

a = cos(theta/2.)

b, c, d = -axis*sin(theta/2.)

return np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)],

[2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)],

[2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]])

def rotation_matrix_weave(theta, axiss,

mat = [1., 0., 0., 0., 1., 0., 0., 0., 1.]):

support = "#include <math.h>"

code = """

double axis [3] = {axiss[0], axiss[1], axiss[2]};

double x = sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] *

axis[2]);

double a = cos(theta / 2.0);

double b = -(axis[0] / x) * sin(theta / 2.0);

double c = -(axis[1] / x) * sin(theta / 2.0);

double d = -(axis[2] / x) * sin(theta / 2.0);

mat[0] = a*a + b*b - c*c - d*d;

mat[1] = 2 * (b*c - a*d);

mat[2] = 2 * (b*d + a*c);

mat[3*1 + 0] = 2*(b*c+a*d);

mat[3*1 + 1] = a*a+c*c-b*b-d*d;

mat[3*1 + 2] = 2*(c*d-a*b);

mat[3*2 + 0] = 2*(b*d-a*c);

mat[3*2 + 1] = 2*(c*d+a*b);

mat[3*2 + 2] = a*a+d*d-b*b-c*c;

"""

weave.inline(code, arg_names=["mat", "axiss", "theta"],

support_code = support)

return np.reshape(mat, (3,3))

def RotateFunc(X, Y, Z, theta, axis=[1.0, 0.0, 0.0]):

x, y = np.meshgrid(X, Y)

z = np.outer(Z, Z)

R = rotation_matrix_weave(theta, axis)

rx, ry, rz = np.zeros(x.shape), np.zeros(y.shape), np.zeros(z.shape)

for i in range(x.shape[0]):

for j in range(x.shape[1]):

tmpx, tmpy, tmpz = np.matmul(R, [x[i, j], y[i, j], z[i, j]])

rx[i,j]+=tmpx

ry[i,j]+=tmpy

rz[i,j]+=tmpz

return rx, ry, rz

## ···

###############################################

### The plot

###############################################

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

xmin, xmax = -5, 5

ymin, ymax = -5, 5

theta = 0.3

N = 1000

x = np.linspace(xmin, xmax, N)

y = np.linspace(ymin, ymax, N)

z = 5.*stats.multivariate_normal().pdf(x)

xp = np.linspace(xmin, xmax, N)

yp = np.linspace(xmin, xmax, N)

zp = np.zeros(x.shape)

rx, ry, rz = RotateFunc(x, y, z, theta)

rx2, ry2, rz2 = RotateFunc(xp, yp, zp, theta)

ax.plot_surface(rx, ry, rz, color='red',alpha=0.65, linewidth=1)

#cset = ax.contourf(rx[:500], ry[:500], rz[:500], zdir='x', offset=-4,

cmap=cm.coolwarm)

cset = ax.contourf(rx, ry-ry2, rz, zdir='y', offset=-4, cmap=cm.coolwarm)

cset = ax.contourf(rx, ry, rz-rz2, zdir='z', offset=-4, cmap=cm.coolwarm)

plt.show()

plt.pcolormesh(rx, ry, rz-rz2)

plt.show()

plt.pcolormesh(rx, rz, ry-ry2)

plt.show()

-------------- next part --------------

An HTML attachment was scrubbed...

URL: <http://mail.python.org/pipermail/matplotlib-users/attachments/20170125/e4bd2b9d/attachment-0001.html>