Generate windrose polar plot from frequency table

295 Views Asked by At

I'm trying to port some code from MATLAB to Python, this snippet should generate two identical windrose plot starting either from raw data of direction and speed for waves of from already calculated frequencies for directions. The first case (with timeseries=1) works as expected, the second one does not. For the first case I used windrose module, but it is my understanding that it cannot generate plots from frequencies, so I'm trying manually using matplotlib but something is wrong. Any hints?

import pandas as pd
import numpy as np
import datetime
from windrose import WindroseAxes
import matplotlib.pyplot as plt
import matplotlib
from scipy.interpolate import interp1d

timeseries = 0
var1_label = 'Hs (m)'  # Choose between Hs (m), Tp (s), ...

# Input
if timeseries == 1:
    filename = '097.00W022.00N.DAT'
    day_start = datetime.datetime(2007, 1, 13)
    day_end = datetime.datetime(2012, 5, 22)

    var1_col = 5
    settoriH = np.arange(1, 4, 0.5)

    theta_col = 7

    nsettori = 12

    A = pd.read_csv(filename, delim_whitespace=True, header=None, comment='*', encoding='cp1252')

    years = A.iloc[:, 0]
    months = A.iloc[:, 1]
    days = A.iloc[:, 2]
    hours = A.iloc[:, 3]
    minutes = A.iloc[:, 4]
    seconds = np.zeros_like(minutes)

    t = pd.to_datetime(dict(year=years, month=months, day=days, hour=hours, minute=minutes, second=seconds))

    Var1 = A.iloc[:, var1_col]
    theta = A.iloc[:, theta_col]

    theta[theta > 1000] = np.nan

    # Crop time series within selected time
    selected_indexes = (t >= day_start) & (t < day_end)
    Var1 = Var1[selected_indexes]
    theta = theta[selected_indexes]
    months = months[selected_indexes]
    t = t[selected_indexes]

    ax = WindroseAxes.from_ax()
    bins = np.arange(0, settoriH[-1] + 1, 0.5)
    ax.bar(theta, Var1, normed=True, opening=0.8, edgecolor='white', bins=np.arange(0, settoriH[-1] + 1, 0.5), nsector=nsettori)
    ax.set_legend()
    plt.show()

elif timeseries == 0:
    filename = 'test_sca.csv'

    A = pd.read_csv(filename, header=None, skiprows=2)
    settoriH0 = A.iloc[0, 1:-1]
    theta_sectors = A.iloc[1:-1, 0]
    theta_sectors = theta_sectors.astype(float)
    nsettori = len(theta_sectors)

    probability = A.iloc[1:-1, 1:-1]
    count2 = np.cumsum(probability, axis=1)
    settoriH = np.arange(1, 4, 0.5)

    count3 = interp1d(settoriH0, count2.T, kind='linear', axis=0)(settoriH)

    count3 = count3.T

    count3 = np.append(count3, count2.iloc[:, -1].to_frame(), axis=1)

    directions = np.arange(0, 360, 30)
    speeds = np.arange(1, 8, 1)

    count3 = pd.DataFrame(count3, index=directions, columns=speeds)

    ax = plt.subplot(111, polar=True)
    cmap = matplotlib.cm.get_cmap("jet", len(count3.columns))  # Color map for different speeds

    for i, speed in enumerate(count3.columns):
        radii = count3[speed]
        theta = np.deg2rad(count3.index -90 )
        width = 2 * np.pi / len(count3.index)
        bottom = np.zeros_like(radii) if i == 0 else count3[count3.columns[i - 1]].values
        bars = ax.bar(theta, radii, bottom=bottom, width=width, color=cmap(i))
    plt.show()

This is the code I'm trying to port in Python, specifically for the frequency part I comment this part from WindRose function:

% %% Bin intensities
% if isempty(vwinds)                                                         % If user did not specify the wind speeds he/she wants to show
%     if ~isempty(WindSpeedRound)                                            % If user did specify the rounding value
%         if isempty(NSpeeds); NSpeeds = 6; end                              % Default value for NSpeeds if not user-specified
%         vmax      = ceil(max(speed)/WindSpeedRound)*WindSpeedRound;        % Max wind speed rounded to the nearest whole multiple of WindSpeedRound (Use round or ceil as desired)
%                     if vmax==0; vmax=WindSpeedRound; end;                  % If max wind speed is 0, make max wind to be WindSpeedRound, so wind speed bins are correctly shown.
%         vwinds    = linspace(0,vmax,NSpeeds);                              % Wind speeds go from 0 to vmax, creating the desired number of wind speed intervals
%     else                                                                   % If user did nor specify the rounding value
%         figure2 = figure('visible','off'); plot(speed);                    % Plot wind speed
%         vwinds = get(gca,'ytick'); delete(figure2);                        % Yaxis will automatically make divisions for us.
%         if ~isempty(NSpeeds)                                               % If a number of speeds are specified
%             vwinds = linspace(min(vwinds),max(vwinds),NSpeeds);            % create a vector with that number of elements, distributed along the plotted windspeeds. 
%         end 
%     end
% end
% 
% %% Histogram in each direction + Draw
% count     = PivotTableCount(N,n,vwinds,speed,dir,NumberElements); 

and call the function like this:

filename='test_sca.csv';

A=readmatrix(filename); settoriH0=A(1,2:end-1); theta_sectors=A(2:end-1,1); nsettori=length(theta_sectors);

probability=A(2:end-1,2:end-1); count2=cumsum(probability,2);

settoriH=1:0.5:3.5;

count3=interp1(settoriH0,count2',settoriH)';

count3(:,end+1)=count2(:,end);


[fig,count,speeds,directions,Table]=WindRose(count3,'FreqLabelAngle',30,'nDirections',nsettori,'vWinds',[0 settoriH],'TitleString',{'Polar plot';' '},'LabLegend',var1_label,'LegendVariable','H_s','AngleNorth',0,'AngleEast',90);
0

There are 0 best solutions below