Polar contour plot with irregular 'bins'

So I have sky data in azimuth and elevation (imported from matlab) that is 'binned' in 3 deg x 3 deg samples. So the matrix is 30 x 120 (3x30=90 deg of el, 3x120=360 of az). Each elevation, 0, 3, 6, ..., 87, has 120 azimuth values associated with it...

However, as you get closer to the zenith you really don't want 120 cells (az cols): As you move up in el, 3 deg of az represents a larger and larger portion of the sky. The final row of the data only has 3 non-zero values that essentially map to az 0-120, 120-240, and 240-360 @ an elevation of 87-90 deg.

So the matrix contains more and more 0's as you go from row 0 (which is el 0-3 deg) to row 29 (el 87-90 deg). Data is here.

I want to plot this on a polar plot, but I don't know how to use contour or pcolor to plot these irregular grids (where the color is set by the value in the matrix).

I've tried:

epfd = np.loadtxt("data.txt")
azimuths = np.radians(np.arange(0, 360, 3))
els      = np.arange(0,90,3)

r, theta = np.meshgrid(els, azimuths)

fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))
im = ax.contourf(theta, r,epfd.T,norm=mpl.colors.LogNorm())
ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2.0)
ax.invert_yaxis()

But that results in

Thanks for your ideas.

Here is how I usually do this:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

data = np.loadtxt(r"C:\Users\serge.degossondevare\Downloads\data.txt")

def create_custom_bins(data):
    azimuths = []
    values = []
    # Offset elevations such that bin edge is at tick for that value
    elevations = np.linspace(0, 90, data.shape[0]) + 90/data.shape[0]/2
    num_values = np.count_nonzero(data, axis=1) # Number of nonzero values in each row, array
    for i in range(data.shape[0]):
        # Create rows of bins that map to the shape of the elevation data - 30 rows
        # Each row has a number of non-zero entries... These are how many bins we want
        # To align the az bins on the left edge we need to subtract half the binsize
        az_bins = np.linspace(0, 360, num_values[i], endpoint=False) + 360/num_values[i]/2
        azimuths.append(az_bins) 
        values.append(data[i, :num_values[i]])
    return elevations, azimuths, values

elevations, azimuths, values = create_custom_bins(data)
min = values[values>0].min()
max = values.max()

fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))

for el, az, val in zip(elevations, azimuths, values):
    theta, r = np.meshgrid(np.radians(az), [el, el + 3])
    val = np.tile(val, (2, 1))
    # Important to set min max since each row will have a different scale based on that rows min/max.
    im = ax.pcolormesh(theta, r, val, norm=mpl.colors.LogNorm(vmin=min,vmax=max), shading='auto')

ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2.0)
ax.invert_yaxis()

cbar = fig.colorbar(im, ax=ax, orientation='vertical')
cbar.set_label('PFD [W/m2/s/Hz]')

plt.show()

which gives

If you want a smoother plot, a variation of the offered solution is

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy.interpolate import griddata

data = np.loadtxt(r"C:\Users\serge.degossondevare\Downloads\data.txt")

def create_custom_bins(data):
    azimuths = []
    values = []
    elevations = np.linspace(0, 90, data.shape[0]) + 90/data.shape[0]/2
    num_values = np.count_nonzero(data, axis=1) # Number of nonzero values in each row, array
    for i in range(data.shape[0]):
        # Create rows of bins that map to the shape of the elevation data - 30 rows
        # Each row has a number of non-zero entries... These are how many bins we want
        # To align the az bins on the left edge we need to subtract half the binsize
        az_bins = np.linspace(0, 360, num_values[i], endpoint=False) + 360/num_values[i]/2
        azimuths.append(az_bins) 
        values.append(data[i, :num_values[i]])
    return elevations, azimuths, values

elevations, azimuths, values = create_custom_bins(data)

azimuth_flat = np.hstack([az for az in azimuths])
elevation_flat = np.hstack([[el]*len(az) for el, az in zip(elevations, azimuths)])
value_flat = np.hstack(values)

azimuth_pad = np.hstack([azimuth_flat, azimuth_flat + 360, azimuth_flat - 360])
elevation_pad = np.hstack([elevation_flat, elevation_flat, elevation_flat])
value_pad = np.hstack([value_flat, value_flat, value_flat])

azimuth_fine = np.linspace(0, 360, 361)
elevation_fine = np.linspace(0, 90, 180)
azimuth_grid, elevation_grid = np.meshgrid(azimuth_fine, elevation_fine)

value_grid = griddata((azimuth_pad, elevation_pad), value_pad, (azimuth_grid, elevation_grid), method='linear')

fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))

theta = np.radians(azimuth_grid)
r = elevation_grid

im = ax.pcolormesh(theta, r, value_grid, norm=mpl.colors.LogNorm(), shading='auto', cmap='viridis')

ax.set_theta_direction(-1)
ax.set_theta_offset(np.pi / 2.0)
ax.invert_yaxis()

cbar = fig.colorbar(im, ax=ax, orientation='vertical')
cbar.set_label('PFD [W/m2/s/Hz]')

plt.show()