Value error: expected input data X have 1 features, but got 2 features during Gaussian Mixture Models

I am currently self learning how to use gaussian mixture models to detect anomalies, but facing some problems in it. I have looked at many blogs but they dont seem to explain what each line signifies. I am trying to figure out which point is an outlier using these 2 variables (r_max, b_max).

This is my dataset (180 rows, 4 columns):

r_max | b_max | SPAD | model
255.0 | 46.0  | 35.1   | Redmi 5A
198.0 | 36.0  | 32.5   | Vivo 1820
237.0 | 145.0 | 27.1   | CPH1920

Note: my r_max and b_max ranges from 0 to 255

Firstly, i filtered out model since it it a string and we can't fit() strings right? I also filtered away SPAD since it is not needed. Then, i changed to a multi-dimensional array to use fit(), seen from this post What is the correct way to fit a gaussian mixture model to single feature data?. Now i will fit() to GaussianMixture(). Not sure how do people know how many n_components to specify?

dataf = df[['r_max', 'b_max']] 
dataf = np.array(dataf).reshape(-1,1) 

gmm = mixture.GaussianMixture()
gmm.fit(dataf)

I get this when i print gmm.fit(dataf):

GaussianMixture()

Now, i will plot the contour plots to see which points are outliers. I specified np.linspace(0, 255) since r_max and b_max varies from 0 to 255.

X, Y = np.meshgrid(np.linspace(0, 255), np.linspace(0, 255))  
XX = np.array([X.ravel(), Y.ravel()]).T

Z = gmm.score_samples(XX)
Z = Z.reshape(X.shape)

CS = plt.contour(X, Y, Z, norm=LogNorm(vmin=1.0, vmax=100.0))
CB = plt.colorbar(CS, shrink=0.8, extend='both')
plt.scatter(dataf['r_max'].values,dataf['b_max'].values)
plt.title('log-likelihood trained using GMM')
plt.axis('tight')
plt.show()

However, i am faced with this error in this line:

--> Z = gmm.score_samples(XX)
ValueError: Expected the input data X have 1 features, but got 2 features

when i print XX i get this:

[[  0.           0.        ]
 [  5.20408163   0.        ]
 [ 10.40816327   0.        ]
 ...
 [244.59183673 255.        ]
 [249.79591837 255.        ]
 [255.         255.        ]]

When i print X and i got this:

[[  0.           5.20408163  10.40816327 ... 244.59183673 249.79591837
  255.        ]
 [  0.           5.20408163  10.40816327 ... 244.59183673 249.79591837
  255.        ]
 [  0.           5.20408163  10.40816327 ... 244.59183673 249.79591837
  255.        ]
 ...
 [  0.           5.20408163  10.40816327 ... 244.59183673 249.79591837
  255.        ]
 [  0.           5.20408163  10.40816327 ... 244.59183673 249.79591837
  255.        ]
 [  0.           5.20408163  10.40816327 ... 244.59183673 249.79591837
  255.        ]]

when i print Y i get this:

[[  0.           0.           0.         ...   0.           0.
    0.        ]
 [  5.20408163   5.20408163   5.20408163 ...   5.20408163   5.20408163
    5.20408163]
 [ 10.40816327  10.40816327  10.40816327 ...  10.40816327  10.40816327
   10.40816327]
 ...
 [244.59183673 244.59183673 244.59183673 ... 244.59183673 244.59183673
  244.59183673]
 [249.79591837 249.79591837 249.79591837 ... 249.79591837 249.79591837
  249.79591837]
 [255.         255.         255.         ... 255.         255.
  255.        ]]

Not sure why i have this error? I am also wondering if its possible to do a contour plot using just 1 variabler_max?

Edits: I updated my codes as per @ronkov's suggestions:

dataf = final_df[['r_max',  'b_max', 'SPAD', 'model']]

dataf = dataf[['r_max', 'b_max']].values

gmm = mixture.GaussianMixture()
gmm.fit(dataf)

X, Y = np.meshgrid(np.linspace(0, 300), np.linspace(0, 300))
XX = np.array([X.ravel(), Y.ravel()]).T
Z = gmm.score_samples(XX)
Z = Z.reshape(X.shape)
# LogNorm only accept positive values, plot -Z
CS = plt.contour(X, Y, -Z, norm=LogNorm(vmax = 300.0), levels=np.logspace(0, 3, 10))
CB = plt.colorbar(CS, shrink=1.0, extend='both')
plt.scatter(dataf[:,0], dataf[:,1], marker = "x", cmap='viridis')
plt.title('log-likelihood trained using GMM')
plt.xlabel('r_max')
plt.ylabel('b_max')
plt.show()

This was what i printed: