PySAL

import numpy as np
# Global Moran's I on a 6x6 grid with rook contiguity
rng=np.random.default_rng(0); g=np.zeros((6,6)); g[:3,:3]=10; g+=rng.normal(0,1,(6,6)); x=g.ravel(); n=x.size
W=np.zeros((n,n))
for i in range(6):
    for j in range(6):
        for di,dj in [(1,0),(-1,0),(0,1),(0,-1)]:
            a,b=i+di,j+dj
            if 0<=a<6 and 0<=b<6: W[i*6+j,a*6+b]=1
z=x-x.mean(); I=n/W.sum()*(z@W@z)/(z@z)
perm=[(lambda s:n/W.sum()*(s@W@s)/(s@s))(rng.permutation(z)) for _ in range(499)]
print("Moran's I =",round(I,3)," pseudo p =",round((np.sum(np.array(perm)>=I)+1)/500,3))