sorted_coo = coo[sortidx]
unqID_mask = np.append(True, np.any(np.diff(sorted_coo, axis=0), axis=1))
ID = unqID_mask.cumsum() - 1
unq_coo = sorted_coo[unqID_mask]
average_values = np.bincount(ID, values[sortidx]) / np.bincount(ID)
Numpy: Average of values corresponding to unique coordinate positions