Hi,
I have a Python script that generates a 3D plot from a Excel .csv file with three columns (Phi, Theta and Power). I'm pretty new with Python. Now I have to modify the code so that Phi have intervalls up to 360 and Theta up to 180. I have tried couple of times now, and cant see the problem. Basically...the logic og Theta and Phi have to be oposite. otherwies the code is working fine. Pretty sure I have to make a change in "deg read" or "def interpolate_x". I have tried to just switch the names og "Theta" and "Phi", but that does not help.
Any help would be strongly appreciated.
This is the code:
I have a Python script that generates a 3D plot from a Excel .csv file with three columns (Phi, Theta and Power). I'm pretty new with Python. Now I have to modify the code so that Phi have intervalls up to 360 and Theta up to 180. I have tried couple of times now, and cant see the problem. Basically...the logic og Theta and Phi have to be oposite. otherwies the code is working fine. Pretty sure I have to make a change in "deg read" or "def interpolate_x". I have tried to just switch the names og "Theta" and "Phi", but that does not help.
Any help would be strongly appreciated.
This is the code:
import pandas as pd
import click
import matplotlib
#matplotlib.use('TkAgg')
# This import registers the 3D projection, but is otherwise unused.
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
import os
import ast
class PythonLiteralOption(click.Option):
def type_cast_value(self, ctx, value):
try:
return ast.literal_eval(value)
except:
raise click.BadParameter(value)
def dbm_to_W(db):
"""
Convert decibell-milliWatts to Watts
:param db:
:return:
"""
return 10**(db/10)/1000
def read(filename):
df = pd.read_csv(filename).sort_values(by = ['Theta','Phi'])
phi = df['Phi'].to_numpy(dtype='float')
theta = df['Theta'].to_numpy(dtype='float')
r = df['Power'].to_numpy(dtype='float')
return [r,theta,phi]
def interpolate_x(x, theta, phi, theta_resolution = 5, phi_resolution = 5, close_gap=True):
"""
Apply bicubic interpolation.
Assumptions: theta and phi are in degrees in limits according to ISO conventions
Theta and phi values cover their respective domain, if
this assumption doesn't hold, extrapolation would occur.
:param x: data values
:param theta: elevation values, must cover the 0-180 range
:param phi: azimuth values, must cover the 0-360 range
:param theta_resolution: output elevation resolution
:param phi_resolution: output azimuth resolution
:param close_gap: boolean, when true sets X[360] to X[0]
:return: [X, THETA, PHI] - interpolated values of X and the THETA-PHI interpolation grid
"""
THETA, PHI = np.mgrid[0.:180.000005:theta_resolution,
0:360.000005:phi_resolution]
# X = interpolator(theta_range, phi_range)
points=np.array([theta.ravel(), phi.ravel()]).T
X = griddata(points, x.ravel(), (THETA, PHI), method='cubic')
if close_gap:
X[PHI == 360.] = X[PHI == 0.]
return [X, THETA, PHI]
def arrange_levels(r,theta,phi):
"""
This generates contour lines to be consumed by matplotlib
surface plot. We rely on data being sorted by theta here.
:param r:
:param theta:
:param phi:
:return:
"""
#TODO: maybe support different number of point in each contour
# using masked arrays? Sub below
#
# R = []
# THETA = []
# PHI = []
# i = 0
# begin = 0
# # Not the most efficient implementation, but makes clear what's happening
# # for i in range(len(theta)):
# n = len(theta)
# while True:
# t = curtheta = theta[i]
# while curtheta==t and i<n:
# t = theta[i]
# i+=1;
# THETA.append(theta[begin:i])
# PHI.append(phi[begin:i])
# R.append(r[begin:i])
# if i==n:
# break
#
# return (np.array(R),
# np.array(THETA),
# np.array(PHI))
u = len(np.unique(theta))
return [x.reshape((u,-1)) for x in (r, theta, phi)]
def sphere2cart(r, theta, phi):
"""
Convert spherical polar coordinates to Cartesian coordinates according to
ISO. See https://en.wikipedia.org/wiki/Spherical_coordinate_system#Conventions
:param r:
:param theta:
:param phi:
:return:
"""
t = np.deg2rad(theta)
p = np.deg2rad(phi)
x = r*np.sin(t)*np.cos(p)
y = r*np.sin(t)*np.sin(p)
z = r*np.cos(t)
return (x,y,z)
def _visuv(r,theta,phi):
"""
A utility I used to get insight into data
:param r:
:param theta:
:param phi:
:return:
"""
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Create the mesh in polar coordinates and compute corresponding Z.
# Plot the surface.
ax.contour(theta,phi,r)
ax.plot_surface(theta,phi,r)#, cmap=plt.cm.YlGnBu_r)
ax.set_xlabel('theta')
ax.set_ylabel('phi')
ax.set_zlabel('r')
plt.show()
def set_view_deg(ax,azim, elev):
ax.view_init(azim=azim, elev = elev)
ax.get_proj()
ax.stale = True
ax.figure.canvas.draw_idle()
def axisEqual3D(ax):
"""
Helper to get correct aspect ratio in 3d plots
https://stackoverflow.com/a/19248731/3791466
:param ax:
:return:
"""
extents = np.array([getattr(ax, 'get_{}lim'.format(dim))() for dim in 'xyz'])
sz = extents[:,1] - extents[:,0]
centers = np.mean(extents, axis=1)
maxsize = max(abs(sz))
r = maxsize/2
# ax.auto_scale_xyz(*np.column_stack((centers - r, centers + r)))
for ctr, dim in zip(centers, 'xyz'):
getattr(ax, 'set_{}lim'.format(dim))(ctr - r, ctr + r)
def add_3d_axes(ax):
x, y, z = np.zeros((3, 3))
u, v, w = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
ax.quiver(x, y, z, u, v, w, arrow_length_ratio=0.1)
for coords, text in zip((u,v,w),"xyz"):
ax.text(*coords, text )# , zdir=None, **kwargs)
def save_standard_TI_notation(path, fig,ax,prefix=''):
standard_rotations = {
"90_0":(0,0),
"0_0":(-90,90),
"90_270":(-90,0),
"90_90":(90,0),
"90_180":(180,0),
"180_0":(90,-90)
}
for r in standard_rotations:
azim = standard_rotations[r][0]
elev = standard_rotations[r][1]
set_view_deg(ax,azim=azim,elev=elev)
pythonFriendlyPath = path.replace('"\"', '"/"')
filename = r+".png"
filename = os.path.join(pythonFriendlyPath, filename)
fig.savefig(filename)
#fig.savefig(r+".png")
def visualize(path, X,Y,Z,V,cmap_limits = (0,0), prefix="", show=True, show_wireframe =False):
if cmap_limits=='auto' or cmap_limits is None or cmap_limits[0]==cmap_limits[1]:
(slope,offset) = get_autoscale(V,0,1)
bounds=(V.min(),V.max())
else:
# cast to float to be safe
bounds = np.array(cmap_limits,dtype='float')
(slope,offset) = get_autoscale(bounds,0,1)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')#,proj_type = 'ortho')
plt.tight_layout()
plt.axis('scaled')
cmap=plt.get_cmap('jet')
patchcolors = cmap(rescale(V,slope,offset), alpha=1)
surf = ax.plot_surface(X, Y, Z, facecolors=patchcolors,rstride=1, cstride=1, linewidth=0.5)#, cmap=plt.cm.YlGnBu_r)
if show_wireframe:
wf = ax.plot_wireframe(X, Y, Z,rstride=1, cstride=1, linewidth=0.5, color='k')#, cmap=plt.cm.YlGnBu_r)
add_3d_axes(ax)
# set_view_deg(ax,azim=30,elev=90)
ax.set_xlim((-1,1))
ax.set_ylim((-1,1))
ax.set_zlim((-1,1))
plt.axis('off')
# cbax = fig.add_subplot(122)
cbax = plt.axes((0.8,0.5, 0.05, 0.45))
norm = matplotlib.colors.Normalize(vmin=bounds[0], vmax=bounds[1])
cb1 = matplotlib.colorbar.ColorbarBase(cbax, cmap=cmap,
norm=norm,
orientation="vertical")
save_standard_TI_notation(path, fig,ax,prefix=prefix)
if show:
plt.show()
# axisEqual3D(ax)
# import vispy.plot as vp
# from vispy import color
#
# fig = vp.Fig()
# ax = fig[0,0]
# ax.surface(Z)
def rescale(x,slope, offset):
return x*slope+offset
def get_autoscale(x, lower, upper):
xmin = x.min()
xmax = x.max()
slope = (upper-lower)/(xmax-xmin)
offset = upper - slope*xmax
return (slope, offset)
@click.command(context_settings=dict(help_option_names=['-h', '--help']))
@click.option('--dbm_to_watts', is_flag=True, help="if present will convert dbm to watts")
@click.option('--disable_interpolation', is_flag=True , help="if present will disable interpolation with 5 degree resulution")
@click.option('--show/--no-show', default=True, help="whether to show the interactive window")
# @click.option('--cmap_limits', cls=PythonLiteralOption, default='"auto"')#,help="preset colormap limits")
@click.option('--cmap_limits', type=(float, float), default=(0,0))#,help="preset colormap limits")
@click.argument('filename', nargs=-1)#,help="paths to csv files to process")
@click.argument('path', nargs=1)#,help="path to where file should be saved")
def main(filename, path, dbm_to_watts = False, disable_interpolation=False, cmap_limits='auto', show=True):
"""
:param filename:
:param path:
:param dbm_to_watts:
:param disable_interpolation:
:param cmap_limits: Limits for the colormap. Use 'auto' or a tuple of lower and upper bounds
in data units.
:return:
"""
for f in filename:
data = read(f)
if dbm_to_watts:
db = data[0] #save dbm to use in the legend adn coloring
data[0] = dbm_to_W(db) #convert to Watts
arranged = arrange_levels(*data)
if disable_interpolation:
interpolated = arranged
else:
interpolated= interpolate_x(*arranged)
slope, offset= get_autoscale(interpolated[0], lower=0, upper=1)
r = rescale(interpolated[0], slope, offset)
xyz = sphere2cart(r,interpolated[1],interpolated[2])
prefix = os.path.splitext(f)
visualize(path, *xyz,interpolated[0],cmap_limits=cmap_limits,prefix=prefix, show=show)
if __name__=='__main__':
main()