Julia, 200 bytes

f(X)=while true
Q=map(x->Tuple(rand(findall(X.==x))),unique(X))
all(map(A->all(map(q->count(A.==q).==1,A)),map(n->getindex.(Q,n),(1,2))))&&!in(2,(sum(abs.(x.-y)) for x in Q for y in Q))&&(return Q)end

Attempt This Online!

This one randomly tries solutions (Q) until they satisfy the conditions. Instead of looking for diagonally adjacent queens, it just looks for queens at a distance of exactly two which includes those and also some orthogonal (but none that are permitted).

JavaScript (ES11), 163 bytes

Expects a square matrix of positive integers. Updates the matrix in-place, putting \$0\$'s at queen positions.

f=(m,N,C,R)=>R+1>>m.length|m.some((r,y)=>r.some((n,x)=>![N>>n,C>>x,R>>y,r[x]=0].some((v,i)=>v&1|m[y+1-i*2%4]?.[x+1-2%~i]<1)&&f(m,N|1<<n,C|1<<x,R|1<<y)||![r[x]=n]))

Attempt This Online! (easy 5×5 test case, instantaneous)

Attempt This Online! (harder 9×9 test case, ~20 seconds)

Commented

f = (                     // f is a recursive function taking:
  m,                      //   m[] = input matrix
  N, C, R                 //   N, C, R = bit-masks for regions, columns, rows
) =>                      //
R + 1                     // force a truthy result if R + 1 is equal to ...
>> m.length |             // ... 2 ** m.length (one queen per row -> success)
m.some((r, y) =>          // for each row r[] at index y in m[]:
  r.some((n, x) =>        //   for each value n at index x in r[]:
    ![                    //     we want to know if there's already a queen in:
      N >> n,             //       the region n
      C >> x,             //       or the column x
      R >> y,             //       or the row y
      r[x] = 0            //     attempt to put a queen at the current position
    ]                     //
    .some((v, i) =>       //     for each value v at index i in the list above:
      v & 1 |             //       test the LSB of v
      m[                  //       test the diagonally touching square at:
        y + 1 - i * 2 % 4 //         y + [ +1, -1, +1, -1 ][i]
      ]?.[                //
        x + 1 - 2 % ~i    //         x + [ +1, +1, -1, -1 ][i]
      ] < 1               //       is it a queen?
    ) &&                  //     end of some() / invert the result
    f(                    //     if true, do a recursive call:
      m,                  //       pass m[]
      N | 1 << n,         //       pass the updated region bit-mask
      C | 1 << x,         //       pass the updated column bit-mask
      R | 1 << y          //       pass the updated row bit-mask
    )                     //     end of recursive call
    ||                    //     unless all of the above is true ...
    ![r[x] = n]           //     ... restore r[x] to n
  )                       //   end of inner some()
)                         // end of outer some()

Jelly,  23  18 bytes

...brute-force method

ŒĠŒpZQ€Ƒ$Ƈạþ`’ḄȦƲƇ

A monadic Link that accepts a list of list of comparable values, e.g. integers (region identifiers) and yields a list of solutions, each as a list of coordinate pairs, [row, column] ([1,1] at top-left).

Try it online! (Third test case, since it completes.)

The first test case runs in just under seven minutes and consumes 12.7 GB of RAM locally (!)

How?

ŒĠŒpZQ€Ƒ$Ƈạþ`’ḄȦƲƇ - Link: list of lists of integers, Grid
ŒĠ                 - group multidimensional indices of Grid by their values
  Œp               - Cartesian product -> all coordinate selections
                                          of one Queen per region
         Ƈ         - keep those for which:
        $          -   last two links as a monad:
    Z              -     transpose
       Ƒ           -     is that invariant under?:
     Q€            -       deduplicate each
                      -> all Queen coordinate selections with one Queen per
                         region, per row, and per column
                 Ƈ - keep those for which:
                Ʋ  -   last four links as a monad:
          ạþ`      -     table of absolute difference with itself
                          -> [[[RowDelta, ColumnDelta], ...], ...]
             ’     -     decrement every delta
              Ḅ    -     convert every pair from binary *
               Ȧ   -     any and all? - i.e. no Queens in the same neighbourhood

* § (sums) would work of course, but it turns out to be slightly slower!

Charcoal, 59 bytes

⊞υＥθＳＦＬθ≔ΣＥυΣＥκＥ⌕ＡμＩ⊕ιＥκ⭆ρ⎇∨∨⁼ςν⁼τξ⁼²⁺↔⁻ςν↔⁻τξ§.Q∧⁼ςν⁼τξσυυ

Attempt This Online! Link is to verbose version of code. Outputs all found solutions. Explanation:

⊞υＥθＳ

Start a search for solutions with the input position.

ＦＬθ

Loop over each digit (ideally the digits would be [0, N) which would save a byte and also allow for N=10).

≔ΣＥυΣＥκＥ⌕ＡμＩ⊕ιＥκ⭆ρ⎇∨∨⁼ςν⁼τξ⁼²⁺↔⁻ςν↔⁻τξ§.Q∧⁼ςν⁼τξσυ

Loop over each position, for each of the remaining rows and columns of the next digit (if any), create a new position with that position replaced by Q and any square that is no longer valid replaced by ., collecting all of the resulting positions.

υ

Output the solutions.