Geographical binning with HISTA not creating squares

Hello everybody!

My data contains latitude, longitute, time and multiple other variables. I need to bin this data for my needs.

>> head(w_files{1})

             time               lon       lat      pci_lte    sinr_lte    cinr_lte    rsrp_lte    prb_dl_lte    prb_ul_lte    bler_dl_lte    bler_ul_lte    mod_dl_0_lte    mod_dl_1_lte    mod_ul_lte    earfnc_dl_lte    earfnc_ul_lte    tp_mac_dl_lte    tp_mac_ul_lte    tp_pdsch_lte    tp_pusch_lte    tp_pdcp_dl_lte    tp_pdcp_ul_lte        tm            tx_rx        bw_dl_lte    bw_ul_lte       band          ca_dl         ca_ul                ca_comb              
    _______________________    ______    ______    _______    ________    ________    ________    __________    __________    ___________    ___________    ____________    ____________    __________    _____________    _____________    _____________    _____________    ____________    ____________    ______________    ______________    __________    _____________    _________    _________    ___________    __________    __________    _________________________    

    2020-12-08 11:48:13.000    14.563    50.035      NaN        NaN          0          NaN           NaN           NaN           NaN            NaN         {0×0 char}      {0×0 char}     {0×0 char}         NaN               NaN              NaN              NaN              NaN            NaN             NaN               NaN          {0×0 char}    {0×0 char   }       NaN          NaN       {0×0 char }    {0×0 char}    {0×0 char}    {0×0 char               }    
    2020-12-08 11:48:14.000    14.563    50.035      374        NaN          0          -61         41.49          0.74             9              5         {'256QAM'}      {0×0 char}     {0×0 char}        6200             24200           238.43            0.002           238.43          0.327               0                 0          {'TM4'   }    {'MIMO(4x2)'}        10           10       {'Band 20'}    {'3CA'   }    {'NonCA' }    {'Band 20+Band 7+Band 3'}     
    2020-12-08 11:48:15.000    14.563    50.035      374         11          0          -61          43.3           0.7             7              2         {'256QAM'}      {0×0 char}     {0×0 char}        6200             24200           256.08            0.002           256.08          0.314               0                 0          {'TM4'   }    {'MIMO(4x2)'}        10           10       {'Band 20'}    {'3CA'   }    {'NonCA' }    {'Band 20+Band 7+Band 3'}   
    2020-12-08 11:48:16.000    14.563    50.035      374         10          0          -61         44.31          0.69             6              2         {'256QAM'}      {0×0 char}     {0×0 char}        6200             24200            257.2            0.001            257.2          0.324               0                 0          {'TM4'   }    {'MIMO(4x2)'}        10           10       {'Band 20'}    {'3CA'   }    {'NonCA' }    {'Band 20+Band 7+Band 3'}    
    2020-12-08 11:48:17.000    14.563    50.035      374         10          0          -61         43.57          0.81             6              0         {'256QAM'}      {0×0 char}     {0×0 char}        6200             24200           252.94            0.001           252.94          0.357               0                 0          {'TM4'   }    {'MIMO(4x2)'}        10           10       {'Band 20'}    {'3CA'   }    {'NonCA' }    {'Band 20+Band 7+Band 3'}    
    2020-12-08 11:48:18.000    14.563    50.035      374          9          0          -63         44.93          0.79             7              0         {'64QAM' }      {0×0 char}     {0×0 char}        6200             24200           226.24            0.002           226.24           0.35               0                 0          {'TM4'   }    {'MIMO(4x2)'}        10           10       {'Band 20'}    {'3CA'   }    {'NonCA' }    {'Band 20+Band 7+Band 3'}    
    2020-12-08 11:48:19.000    14.563    50.035      374          7          0          -62         44.45           0.8            11              0         {'16QAM' }      {0×0 char}     {0×0 char}        6200             24200           191.17            0.002           191.17          0.359               0                 0          {'TM4'   }    {'MIMO(4x2)'}        10           10       {'Band 20'}    {'3CA'   }    {'NonCA' }    {'Band 20+Band 7+Band 3'}   
    2020-12-08 11:48:20.000    14.563    50.035      374          2          0          -66         44.18          0.82            12              2         {'16QAM' }      {0×0 char}     {0×0 char}        6200             24200            177.6            0.002            177.6          0.372               0                 0          {'TM4'   }    {'MIMO(4x2)'}        10           10       {'Band 20'}    {'3CA'   }    {'NonCA' }    {'Band 20+Band 7+Band 3'} 

I am using HISTA to find a 20x20 meters grid and applying it on the data sets.

[latbin, lonbin] = hista(w_files{ii}.lat,w_files{ii}.lon,0.0004);% bin at 20x20m,hista binning + computing binsID
[w_files{ii}.latEq, w_files{ii}.lonEq] = grn2eqa(w_files{ii}.lat,w_files{ii}.lon);% Convert coordinates to equidistant cartesian coordinates
[latbinEq, lonbinEq] = grn2eqa(latbin, lonbin);% Convert coordinates to equidistant cartesian coordinates
dist = pdist2([w_files{ii}.lonEq,w_files{ii}.latEq],[lonbinEq, latbinEq]);% Compute distance between each coordinate and each bin-center
[~, w_files{ii}.bin20] = min(dist,[],2);% Add bin ID numbers to table

When I run a figure with plot_google_map I am not getting squares but rectangles. Right now, I am not sure if I am just displaying it wrong or I have made some mistake in the binning. I have added one sample file for you guys to try. (You will need plot_google_map function to have the map picture but this issue is clearly visible without it as well).

   figure()
   hold on
   scatter(lonbin, latbin, 200, 1:numel(lonbin), 'Marker','*','LineWidth',2) % bin centers
   scatter(w_files{ii}.lon,w_files{ii}.lat, 50,w_files{ii}.bin20,'filled')
   cmap = colorcube(255);
  colormap(cmap(1:end-10,:))
   xlabel('longitude'); ylabel('latitude')
   latbinUnq = unique(latbin); 
   lonbinUnq = unique(lonbin);
   set(gca, 'xtick', lonbinUnq(2:end)-diff(lonbinUnq)/2,'ytick',latbinUnq(2:end)-diff(latbinUnq)/2)
   grid on
  plot_google_map('Scale',2,'resize',2,'ShowLabels',0)%google_plot settings

What do you think? If you think you know the answer, please be very specific because I am a beginner and I might be dumb.