1

I am plotting a 3D shape in spherical coordinates. In order to rotate it, I am shifting the phi values by 30 deg as phi_lin and phi_rot show in in the following code. I would expect the result in panel 4 to have the same distribution of panel 2, but rigidly shifted to the right by 30 degrees.

I guess, the problem is that plotting function countorf cannot deal with the phi_rot input vector since it is non-monotonic. It is possible to see in panel 3 the discontinuity du the shifting. How can I overcome this problem?

Here a working code:

import glob
import math
import numpy as np
import pandas as pd

import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import LightSource
%matplotlib inline

import itertools

def ellips(THETA,PHI):
    """
    #Definiton of the ellipsoid
    # from https://arxiv.org/pdf/1104.5145.pdf
    """
    a=1; b=2; c=3
    R = (a*b*c) / np.sqrt(b**2*c**2*np.cos(THETA)**2 + c**2*a**2*np.sin(THETA)**2*np.cos(PHI)**2 + a**2*b**2*np.sin(THETA)**2*np.sin(PHI)**2)
    return np.array(R)


nth=13
theta = np.linspace(0, np.pi, nth)

#length = 13
phi_lin=[-180,-150,-120,-90,-60,-30,0,30,60,90,120,150,180]
phi_rot=[-150,-120,-90,-60,-30,0,30,60,90,120,150,180,-180]

THETA_lin, PHI_lin = np.meshgrid(theta, phi_lin)
THETA_rot, PHI_rot = np.meshgrid(theta, phi_rot)

THETA_deg_lin=[el*180/np.pi for el in THETA_lin]
THETA_deg_rot=[el*180/np.pi for el in THETA_rot]
PHI_deg_lin=[el for el in PHI_lin]
PHI_deg_rot=[el for el in PHI_rot]

fig1, ax = plt.subplots(2,2, figsize=(15,15), constrained_layout=True)

ax[0,0].plot(PHI_deg_lin, "o")
ax[0,0].set_xlabel("# element")
ax[0,0].set_ylabel('phi [DEG]')
ax[0,0].set_title("initial coordinates")

ax[0,1].contourf(PHI_deg_lin, THETA_deg_lin, ellips(THETA_deg_lin,PHI_deg_lin).reshape(len(phi_lin),nth))
ax[0,1].set_xlabel('phi [DEG]')
ax[0,1].set_ylabel('theta [DEG]')
ax[0,1].set_title("Original ellipsoind in spherical coordinates")

ax[1,0].plot(PHI_deg_rot, "o")
ax[1,0].set_xlabel("# element")
ax[1,0].set_ylabel('phi [DEG]')
ax[1,0].set_title("shifted coordinates")

ax[1,1].contourf(PHI_deg_rot, THETA_deg_rot, ellips(THETA_deg_rot,PHI_deg_rot).reshape(len(phi_rot),nth))
ax[1,1].set_xlabel('phi [DEG]')
ax[1,1].set_ylabel('theta [DEG]')
ax[1,1].set_title("Original ellipsoind in spherical coordinates")

and the output:

output_edit

UPDATE: I tried to create an interpolation function z=f(x,y) with the rotated coordinates and to plot the new z:

from scipy import interpolate
    
i2d = interpolate.interp2d(theta, phi_rot, ellips(THETA_deg_rot,PHI_deg_rot))
    znew = i2d(theta,phi_lin)
    
    
ax[1,1].contourf(PHI_deg_rot, THETA_deg_rot,znew.reshape(len(phi_rot),nth))

the shifting occurs as you can see in the following output, but the non linearly-spaced x axis prevents to have a smooth contour:

interpolated

any idea how to fix it?

giammi56
  • 154
  • 1
  • 15
  • the main problem is how to deal with the non-monotonic behaviour of the array describing the x-axis – giammi56 Dec 08 '20 at 09:18

1 Answers1

0

The solution has been inspired by this post.

Since contourf doesn´t accept non-linearly-spaced axis, it is necessary to interpolate the rotated data

from scipy import interpolate

i2d = interpolate.interp2d(theta, phi_rot, ellips(THETA_deg_rot,PHI_deg_rot))

evaluate it on the same axis (lin or rot doesn´t matter at this point)

znew = i2d(theta,phi_lin)

and plotting it using the tricontourf with a suitable numner of levels

ax[1,1].tricontourf(np.array(PHI_deg_rot).reshape(-1), np.array(THETA_deg_rot).reshape(-1),znew.reshape(-1),10)

the output is the expected one:

rotated interpolated output

giammi56
  • 154
  • 1
  • 15
  • probably would work using just `ax[1,1].tricontourf(PHI_deg_rot, THETA_deg_rot, ellips(THETA_deg_rot,PHI_deg_rot).reshape(len(phi_rot),nth))` provided that `x` and `y` are `meshgrid` – giammi56 Dec 23 '20 at 13:23