Error "Valueerror: operands could not be broadcast Together with shapes"

Question

I’m getting this mistake:

Traceback (Most recent call last): File "geraEstatistica.py", line 148, in cnplot = map.Contour(x,y,result,clevs,Cmap=plt.cm.jet) File "/home/caio/anaconda3/envs/showcast2/lib/python3.7/site-Packages/mpl_toolkits/Basemap/init.py", line 543, in with_transform Return plotfunc(self,x,y,data,*args,**kwargs) File "/home/caio/anaconda3/envs/showcast2/lib/python3.7/site-Packages/mpl_toolkits/Basemap/init.py", line 3569, in Contour Mask = np.logical_or(ma.getmaskarray(data),xymask) Valueerror: operands could not be broadcast Together with shapes (97,153) (2,2)

When I run the following python script:

  # -*- coding: utf-8 -*-
from netCDF4 import Dataset
#from mpl_toolkits.axes_grid1.inset_locator import inset_axes # Add a child inset axes to this existing axes.
from datetime import datetime, timedelta                     # Library to convert julian day to dd-mm-yyyy
#from cpt_convert import loadCPT                              # Import the CPT convert function
#from matplotlib.colors import LinearSegmentedColormap        # Linear interpolation for color maps
import matplotlib.pyplot as plt                              # Plotting library
#import matplotlib.colors                                     # Matplotlib colors
import numpy as np                                           # Scientific computing with Python
#import cartopy, cartopy.crs as ccrs                          # Plot maps
#import cartopy.io.shapereader as shpreader                   # Import shapefiles
#import time as t                                             # Time access and conversion
import math                                                  # Import math
#from matplotlib.image import imread                          # Read an image from a file into an array
import os                                                    # Miscellaneous operating system interfaces
#import sys                                                   # Import the "system specific parameters and functions" module
#from cartopy.feature.nightshade import Nightshade            # Draws a polygon where there is no sunlight for the given datetime.
#from html_update import update                               # Update the HTML animation
from remap import remap                                      # Import the Remap function
from mpl_toolkits.basemap import Basemap 
#from  netCDF4 import Dataset as open_ncfile
import re
global extent_
global lat
global lon
dpi = 100



def interpreta(path): #interpreta os dados da matriz ACMF Clear Sky Mask

    #start = t.time()
    #nc = open_ncfile('netsacmf_jan_2020/'+str(path))

    #print(nc.variables)


    # Image path
    path = "netsacmf_jan_2020/"+str(path)+"_SEC"

    # Remove the composite id
    path_proc = path[:-4]

    # Read the image
    file = Dataset(path_proc)

    # Read the satellite
    satellite = getattr(file, 'platform_ID')

    # Read the main variable
    variable = list(file.variables.keys())[0]

    # Read the product name
    prod_name = getattr(file, 'title')

    # Desired resolution
    resolution = 8

    # Read the resolution
    band_resolution_km = getattr(file, 'spatial_resolution')
    band_resolution_km = float(band_resolution_km[:band_resolution_km.find("km")])

    # Division factor to reduce image size
    f = math.ceil(float(resolution / band_resolution_km))

    #print(resolution)
    #print(band_resolution_km)
    #print(f)

    # Read the central longitude
    longitude = file.variables['goes_imager_projection'].longitude_of_projection_origin

    # Read the semi major axis
    a = file.variables['goes_imager_projection'].semi_major_axis

    # Read the semi minor axis
    b = file.variables['goes_imager_projection'].semi_minor_axis

    # Calculate the image extent
    h = file.variables['goes_imager_projection'].perspective_point_height
    x1 = file.variables['x_image_bounds'][0] * h
    x2 = file.variables['x_image_bounds'][1] * h
    y1 = file.variables['y_image_bounds'][1] * h
    y2 = file.variables['y_image_bounds'][0] * h

    # Reading the file time and date
    add_seconds = int(file.variables['time_bounds'][0])
    date = datetime(2000,1,1,12) + timedelta(seconds=add_seconds)
    date_formated = date.strftime('%Y-%m-%d %H:%M UTC')
    date_file = date.strftime('%Y%m%d%H%M')
    year = date.strftime('%Y')
    month = date.strftime('%m')
    day = date.strftime('%d')
    hour = date.strftime('%H')
    minutes = date.strftime('%M')

    # Choose the visualization extent (min lon, min lat, max lon, max lat)
    extent_ = [float("-54.0"), float("-25.0"), float("-43.0"), float("-18.0")]

    # Call the reprojection funcion
    grid = remap(path_proc, variable, extent_, resolution, h, a, b, longitude, x1, y1, x2, y2)

    # Read the data returned by the function
    data = grid.ReadAsArray()

    # Convert from int16 to uint16
    data = data.astype(np.float64)
    return data


print('- Mapa de janeiro 2020 - cobertura de nuvem - \n Script iniciado! \n gerando matrizes binárias... \n calculando soma de matrizes... ')
global arquivos;
global soma;
soma = np.array([], dtype = float)
for arquivo in os.walk('netsacmf_jan_2020/'): #lê o diretório que contém os nets e transforma em listas com os nomes dos arquivos
    arquivos = list(arquivo[2])

for arq in arquivos:
    if not re.search('aux.xml', arq):
        net_receptus = interpreta(arq)  #manda para a função interpreta para realizar a leitura do nets e transformar na matriz binária net_receptus
        if (len(soma) == 0): 
            soma = net_receptus
        else:
            soma = np.add(soma,net_receptus)


print("calculando estatística final de janeiro...")

qttd = len(arquivo[2])

resultado = (soma / qttd)

lat = np.array([90, -90], dtype = float)
lon = np.array([0, 357.5], dtype = float)
fig = plt.figure(figsize=(1100/dpi, 1100/dpi), dpi=dpi)
ax  = fig.add_axes([0.1,0.1,0.8,0.9])
map = Basemap(projection='cyl',llcrnrlat= -90.,urcrnrlat= 90.,
              resolution='c',  llcrnrlon=-180.,urcrnrlon=180.)
map.drawcoastlines()
map.drawstates()
map.drawcountries()

map.drawparallels(np.arange( -90., 90.,30.),labels=[1,0,0,0],fontsize=10)
map.drawmeridians(np.arange(-180.,180.,30.),labels=[0,0,0,1],fontsize=10)

x, y = map(*np.meshgrid(lon,lat))
clevs = np.arange(0,100,1)
cnplot = map.contour(x,y,resultado,clevs,cmap=plt.cm.jet)
plt.title('Mapa de indice de combertura de nuvens - Estado de São Paulo - Janeiro 2020')

plt.show()
plt.savefig('mapa_de_indice_de_cobertura_de_nuvens.png', bbox_inches='tight', dpi=dpi)

The purpose of this script is to plot the results matrix, which is the sum matrix divided by the amount of netCDF’s files with cloud coverage data read. The sum matrix is the sum of all binary matrices that are generated by netCDF’s files.

[[0,1,0,1,1,0]] + [[0,1,0,0,1,1]] + ....... = [[[20.40,0,20,60,10]] / number of binary matrices or number of netCDFs files = result matrix with percentages.

I’m trying to plot this result matrix with the percentages on a map that looks like the image below (only what I’m doing for the map of the state of SP):

for this I am using the countr() and countrf() function of matplotlib. but I keep receiving this traceback.

I am using netCDF’s files with ACMF with the cloud cover mask made available through Amazon’s AWS.

The following are several files from: OR_ABI-L2-ACMF-M6_g16_s20200010100216_e20200010109524_c20200010110376.nc each generate a binary matrix. 0 = cloudless and 1 = cloudless.

here is a folder with netCDFs files to be able to run the script: Folder link in drive