Xarray DataArray georeferenced plot only displays the contour of the data on the map

Question

I am trying to plot my satellite data (a .zarr file loaded as an xarray DataSet and subsetted to a DataArray to display only one band data) as a georeferenced plot. When the plot loaded I get the correct map, and only a contour of my data in the correct place, without the coloured data which should fill the plot.

Here is what I get:

Here is the code:

# Open the zarr store
zarr_store = zarr.open(product_path)
zarr_store.tree()

# Get measurements dictionary
def get_measurements_dict(product_store):
    measurement_keys = list(product_store['measurements'].keys())
    dict = {}
    datasets = []
    for key in measurement_keys:
        dataset = xr.open_dataset(f'{product_path}/measurements/{key}', engine='zarr', consolidated=False)
        datasets.append(dataset)
        dict[key] = dataset
    return dict
# datasets
dataset_dict = get_measurements_dict(zarr_store)
dataset_dict

# Select the band you want to plot

class BandSelector:
    def __init__(self, dataset_dict):
        self.bands = list(dataset_dict.data_vars.keys())  # Get a list of band names from the dataset
        self.selected_band = None

    def select_band(self, band):
        self.selected_band = band
        
        # Store the selected band within a variable
        global selected_band_name
        selected_band_name = self.selected_band

band_selector = BandSelector(dataset_dict[selected_dataset_name])
selected_dataset_name
@interact(band=widgets.Dropdown(options=band_selector.bands))
def select_band(band):
    band_selector.select_band(band)
    print("Selected band:", selected_band_name)

# Check selected band data
band[selected_band_name].data

# Use the compute method to put data into memory
data = band[selected_band_name].data.compute()

# Print result
print("data")
print(data)
print("\ndata.shape:")
print(data.shape)


# Impelement dynamic selection
y_dim, x_dim = data.dims
print(y_dim, x_dim)

# Get data from selected band
x_comp = data[x_dim]
y_comp = data[y_dim]
band_choice = band[selected_band_name].data.compute()

# Get bbox
geometry_from_product = np.squeeze(dt.attrs['stac_discovery']['geometry']['coordinates'])

# Define constant for plotting
L1C_PROJECTION = ccrs.epsg(32633)
DESIRED_PROJECTION = ccrs.PlateCarree()
FIGSIZE: tuple[int, int] = (12, 8)
RESOLUTION_CARTOPY: str = '110m'
GEOGRAPHICAL_LIMITS: tuple[int, int, int, int] = (6, 19, 36, 47) # Coords for Italy

# Speed up plot by sampling data every SKIP_EVERY pixels
SKIP_EVERY: int = 10
# Define plotting arguments for Polygon around the area of interest
POLYGON_THICKNESS: int = 1
POLYGON_COLOR: str = 'r'


# Start plotting
_, ax = plt.subplots(subplot_kw={"projection": DESIRED_PROJECTION},
                     figsize=FIGSIZE)

# Plot cartopy geographic information
ax.coastlines(resolution=RESOLUTION_CARTOPY)
ax.add_feature(cf.BORDERS)
ax.add_feature(cf.OCEAN)
ax.add_feature(cf.LAND)
gl = ax.gridlines(draw_labels=True, 
                  crs=DESIRED_PROJECTION)

plt.contourf(band_choice[::SKIP_EVERY, ::SKIP_EVERY],
             transform = L1C_PROJECTION)
poly = mpatches.Polygon(geometry_from_product, 
                        closed=True, 
                        ec=POLYGON_COLOR, 
                        fill=False, 
                        lw=POLYGON_THICKNESS, 
                        transform=DESIRED_PROJECTION)
ax.add_patch(poly)
ax.set_extent(GEOGRAPHICAL_LIMITS, crs=DESIRED_PROJECTION)
cbar = plt.colorbar(orientation="horizontal")
cbar.set_label('b07_20m')
plt.tight_layout()

I expected a plot similar to this: