How to count the sum of at least two consecutive days with same the values, from xarray dataset in python

An easy logic is to do a runsum on the binary field using the rolling sum function with a window size equal to your event length N, identify the Ns with a logical expression and then just add those up. Thus you can do it in these two lines:

event_len=3 # this is N
ds_sum = ds.rolling(time=event_len, center=False).sum()
events = (ds_sum == event_len).astype(int).sum(dim="time")

If this logic is unclear, I actually have a youtube video posted on this on my channel for further clarification.

However, as pointed out in the comment by @Patrick, while this works fine for continuous data, it has a minor error if your data is only for specific months, e.g. July, as you may get a heatwave that consists of the last day of July in year X and then the first two days of July in year X+1 for example!!! To avoid this, we want to effectively turn the 3D [lon,lat,time] array into a 4D [lon,lat,year,day] array and apply the rolling function on the last "day" time dimension for each year separately. Luckily for us, there is the xarray function groupby which is ideal for this and we then just use apply like this.

events_per_year = ds['data'].groupby(ds.time.dt.year).apply(
    lambda x: (x.rolling(time=event_len, center=False).sum() == event_len).astype(int).sum(dim="time")
)

I made some dummy July data for 4 years up to illustrate this for a small 3x3 grid (so this is reproducible after the OPs netcdf link disappears in the future):

import xarray as xr
import numpy as np
import pandas as pd

# make up some dummy July data for a few year on a small grid:
time = pd.date_range(start='2010-07-01', end='2013-07-31', freq='D')
time = time[time.month == 7]
lon = np.arange(0, 3, 1)
lat = np.arange(0, 3, 1)
data = np.random.randint(0, 2, size=(len(time), len(lat), len(lon)))
ds = xr.Dataset(
    {"data": (("time", "lat", "lon"), data)},
    coords={"time": time, "lon": lon, "lat": lat},)

### now we define events of this length:

event_len=3

# this function wraps up 3 steps:
# 1. Calculate the rolling sum for 'event_len' days for each year
# 2. Check where this rolling sum equals 'event_len' which I put 3 days here as per the OP def 
# 3. Sum the true/false 1/0 in time dim
events_per_year = ds['data'].groupby(ds.time.dt.year).apply(
    lambda x: (x.rolling(time=event_len, center=False).sum() == event_len).astype(int).sum(dim="time")
)

print(events_per_year)

This gives the following output

<xarray.DataArray (year: 4, lat: 3, lon: 3)>
array([[[3, 2, 2],
        [0, 6, 8],
        [3, 5, 4]],

       [[8, 2, 8],
        [4, 5, 6],
        [7, 4, 4]],

       [[9, 3, 1],
        [6, 4, 3],
        [2, 2, 7]],

       [[1, 3, 3],
        [2, 6, 8],
        [3, 2, 0]]])
Coordinates:
  * lat      (lat) int64 0 1 2
  * lon      (lon) int64 0 1 2
  * year     (year) int64 2010 2011 2012 2013

I can propose following solution:

#!/usr/bin/env ipython
# ---------------------
import numpy as np
import xarray as xr
from pylab import pcolormesh,show,colorbar,plot,title,legend,subplot,savefig
# -------------------
fin = 'EC7_heatjulZA26_30.nc' # filename in...

dfin = xr.open_dataset(fin) # let us use xarray to read data ...
vin = dfin['tasmax'].values # get only the values ...
ntime,ny,nx = np.shape(vin) # get the dimensions...
# --------------------------------------
dvin = np.diff(vin,axis=0) # let us take the diff along time axis...
# --------------------------------------
counts = np.sum(vin,axis=0) # let us count days with temp over threshold...
pcolormesh(counts);colorbar();title('Number of days over threshold');show() # let us view the map...
# --------------------------------------
dvin[dvin<0] = 0.e0; # let us remove the -1 when the temperature has dropped below the treshold...
icounts = np.sum(dvin,axis=0)
pcolormesh(icounts);colorbar();title('Number of instances over threshold (simple)');savefig('simple.png',bbox_inches='tight');show() # let us view the map...
# let us check:
plot(vin[:,1,0]);title('Number of instances found '+str(icounts[1,0]));show() # if you check, the number of instances is too high -- 9 instead of 6
# ---------------------------------------
# let us calculate correct one:
ntime,ny,nx = np.shape(vin) # get the dimensions...
dvin_org = np.diff(vin,axis=0); # the diff...

dvin = np.concatenate((np.zeros((1,ny,nx)),dvin_org),axis=0); # make the diff and original data same size...
dvin_n = np.concatenate((np.zeros((2,ny,nx)),dvin_org[:-1,:,:]),axis=0); # shift the diff +1 index in time
# ------------------------------------------------------------------------------------------------------
dvin[np.where((dvin_n==1)&(vin==1))] = 1.0 # in original diff, add one to the location, where the derivate does not change -- these are the time moments we actually wanted...
# -------------------------------------------------------
icounts = np.sum(dvin,axis=0) # calculate the correct number of instances over treshold
pcolormesh(icounts);colorbar();title('Number of instances over threshold (complex)');savefig('corrected.png',bbox_inches='tight');show() # let us view the map...
# let us check:
plot(vin[:,2,2]);title('Number of instances found '+str(icounts[2,2]));show()

So, the original calculation gives figure like this:

where the occurrences goes up to 9 times in one grid cell. But, this is overestimated as the time-series look like this (9 is taken from the simple solution):

Problem is that we counted also couple of events where the temperature went over threshold only for one day. Therefore, made additional check with using the shifted derivate, see the code up.

In any case, the final solution comes then like this:

Hope this helps!