how to change the particuler elements of an array

The below proposed solution uses skimage.util.view_as_windows in order to get all of the possible 3x3 views onto the larger array and then looping using them with help of Python loops to set the found hits back to ones after zeroing all the ones in the first step.

It seems that OpenCV, scimage, numpy don't provide a method able to label areas in which a 3x3 square can move around having all the values beneath set to 1's, but I am not sure here. So if you read this and know much about dilations, convolutions, etc. please leave a note pointing me in the right direction allowing to label the area of image by a methond implemented in C-code for better speed on huge arrays.

# https://stackoverflow.com/questions/73649612/how-to-change-the-particuler-elements-of-an-array
# Example array: 
import numpy as np
D_in =np.array([[4,4,4,0,1,1,1,1,1,1,0,0,0,0,1],
                [3,0,1,0,1,1,1,1,1,1,1,1,1,1,0],
                [6,0,0,0,1,1,1,1,1,1,1,1,1,1,0],
                [2,0,0,0,1,1,1,0,1,0,1,1,1,0,0],
                [2,0,1,0,1,1,1,0,1,0,1,0,1,0,0]])
# Requirements: If in the D_in array at least one 3x3 square with 1's 
# is found, replace in D_in all 1's with 0's except the square size and more than square size no changes. 
# Otherwise.
# Expected output :
D_tgt=np.array([[4,4,4,0,1,1,1,1,1,1,0,0,0,0,0],
                [3,0,0,0,1,1,1,1,1,1,1,1,1,0,0],
                [6,0,0,0,1,1,1,1,1,1,1,1,1,0,0],
                [2,0,0,0,1,1,1,0,0,0,1,1,1,0,0],
                [2,0,0,0,1,1,1,0,0,0,0,0,0,0,0]])
# ----------------------------------------------------------------------
import numpy as np
D_out  = np.copy(D_in) # D_out for restoring hit 3x3 ONEs
D_out[D_out==1] = 0     # set all ONEs  to zero 
needle = np.ones((3,3)) # not necessary here, except documentation
# Create all possible 3x3 needle views on larger D_in heystack: 
from skimage.util import view_as_windows as winview
# vvvvvvvvvvvvvvvvvvvvvvvvvvv
all3x3 = winview(D_in,(3,3)) # the CORE of the algorithm
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^
all3x3_rows, all3x3_cols, needle_rows, needle_cols = all3x3.shape
print(f'{all3x3_rows=}, {all3x3_cols=}, {needle_rows=}, {needle_cols=}')
noOfHits = 0 # used also to decide about the output if not found Hits
for row in range(all3x3_rows):
    for col in range(all3x3_cols):
        if np.all(all3x3[row,col,:,:]):
            noOfHits += 1
            # print(f'3x3 Ones at: {row=}, {col=}')
            D_out[row:row+3, col:col+3] = 1
print('--------------------------')
print(D_in)
print('--------------------------')
if noOfHits > 0:
    print(D_out)
    assert D_out.all() == D_tgt.all() # make sure the result is correct
else: 
    print(D_in)

gives on output:

all3x3_rows=3, all3x3_cols=13, needle_rows=3, needle_cols=3
--------------------------
[[4 4 4 0 1 1 1 1 1 1 0 0 0 0 1]
 [3 0 1 0 1 1 1 1 1 1 1 1 1 1 0]
 [6 0 0 0 1 1 1 1 1 1 1 1 1 1 0]
 [2 0 0 0 1 1 1 0 1 0 1 1 1 0 0]
 [2 0 1 0 1 1 1 0 1 0 1 0 1 0 0]]
--------------------------
[[4 4 4 0 1 1 1 1 1 1 0 0 0 0 0]
 [3 0 0 0 1 1 1 1 1 1 1 1 1 0 0]
 [6 0 0 0 1 1 1 1 1 1 1 1 1 0 0]
 [2 0 0 0 1 1 1 0 0 0 1 1 1 0 0]
 [2 0 0 0 1 1 1 0 0 0 0 0 0 0 0]]

The same output can be achieved with the one-liner from the answer by mozway. Here some preliminary code:

import numpy as np
data =np.array([[4,4,4,0,1,1,1,1,1,1,0,0,0,0,1],
                [3,0,1,0,1,1,1,1,1,1,1,1,1,1,0],
                [6,0,0,0,1,1,1,1,1,1,1,1,1,1,0],
                [2,0,0,0,1,1,1,0,1,0,1,1,1,0,0],
                [2,0,1,0,1,1,1,0,1,0,1,0,1,0,0]])
print(data)
print("  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^")
from scipy.signal  import convolve2d
from scipy.ndimage import binary_dilation
b = binary_dilation
c = convolve2d
d = data
k = (3,3)
m = data==1 # ; print( m.astype(int) )
n = np.ones
s = 'same'
w = np.where

than

out=w(m,b(c(d==1,n(k),mode=s)==9,n(k)),d)

and finally

print(out)

You can use a 2D convolution on the 1s with a 3x3 kernel of 1s to identify the centers of the 3x3 squares, then dilate them and restore the non 1 numbers

from scipy.signal import convolve2d
from scipy.ndimage import binary_dilation

# get 1s (as boolean)
m = data==1

kernel = np.ones((3, 3))

# get centers
conv = convolve2d(m, kernel, mode='same')

# dilate and restore the other numbers
out = np.where(m, binary_dilation(conv == 9, kernel).astype(int), data)

print(out)

Alternative, erosion and dilation, similarly to my previous answer:

from scipy.ndimage import binary_dilation, binary_erosion

m = data==1
kernel = np.ones((3, 3))

out = np.where(m, binary_dilation(binary_erosion(m, kernel), kernel), data)

Output:

[[4 4 4 0 1 1 1 1 1 1 0 0 0 0 0]
 [3 0 0 0 1 1 1 1 1 1 1 1 1 0 0]
 [6 0 0 0 1 1 1 1 1 1 1 1 1 0 0]
 [2 0 0 0 1 1 1 0 0 0 1 1 1 0 0]
 [2 0 0 0 1 1 1 0 0 0 0 0 0 0 0]]