The first order differential operators aren't exactly what you put. The author of that post didn't mention that his D operator has to be a vector-valued first order differential operator on vector-valued arguments, so it's more like ## D|_{(x,y,z)}[\vec{A}] = \sum_{i,j,k}^{3} c_{ijk}(x,y,z) \vec{e}_{i} \frac{\partial}{\partial x_{j}} A_{k}(x,y,z) ##.
To test for invariance under translations and rotations, you have to first write how each affects the coordinates and basis vectors.
For translations: ## (x,y,z) \to (x+\epsilon, y, z) ## (it doesn't matter which direction the displacement is in, because I haven't specified what direction x is)
##D|_{(x,y,z)} \to D|_{(x+\epsilon,y,z)}## which implies that ##c_{ijk}(x,y,z) = c_{ijk}(x+\epsilon,y,z)## for all values of ##\epsilon##, or in other words that ##c_{ijk}## is a constant for all x. Since you can repeat this for y and z, you know that ##c_{ijk}## have to be constant for all (x,y,z).
For rotations: The general formula for rotations about any axis is a pain to work with, but we can say that if D is invariant under all rotations then D is invariant under the subset of right-handed rotations through 90 degrees. (This is, in fact, equivalent to saying the D is invariant under the action of infinitesimal rotations, but conceptually simpler).
##Case 1: (x,y,z) \to (-y, x, z) = (x',y',z')##
##(\frac{\partial}{\partial x},\frac{\partial}{\partial y}, \frac{\partial}{\partial z}) \to ( - \frac{\partial}{\partial y}, \frac{\partial}{\partial x}, \frac{\partial}{\partial z}) ##
##(A_{x}, A_{y}, A_{z}) \to (-A_{y}, A_{x}, A_{z})##
##D[\vec{A}] = \sum_{i,j.k=1}^3 \vec{e}_{i} \frac{\partial}{\partial x_{j}} A_{k} = \sum_{i,j.k=1}^3 \vec{e}_{i} \frac{\partial}{\partial x'_{j}} A_{k} ##
From this we obtain, ## c_{311} \frac{\partial}{\partial x} A_{x} + c_{312} \frac{\partial}{\partial x}A_{y} + c_{321} \frac{\partial}{\partial y} A_{x} + c_{322} \frac{\partial}{\partial y} A_{y} = (-1)^{2} c_{311} \frac{\partial}{\partial y} A_{y} + (-1) c_{312} \frac{\partial}{\partial y}A_{x} + (-1) c_{321} \frac{\partial}{\partial x} A_{y} + c_{322} \frac{\partial}{\partial x} A_{x} ##
And finally we get
## (c_{311} - c_{322}) (\frac{\partial}{\partial x}A_{x} - \frac{\partial}{\partial y}A_{y}) + (c_{312} + c_{321}) (\frac{\partial}{\partial x}A_{y}) + (c_{321} + c_{312}) (\frac{\partial}{\partial y}A_{x}) = 0##
Since each of the bracketed combinations of the partials of ##\vec{A}## are arbitrary smooth functions, the only way to solve this equation is if all the scalar coefficients (the parenthesized terms with the c's in them) are equal to 0. There are exactly 2 non-trivial normalized solutions for the coefficients:
Solution 1: ##c_{311} = c_{322} = 1##, ##c_{312} = c_{321} = 0##. (The sum of like partial derivatives of the x and y components of A, i.e. the 2D divergence of A restricted to the z-plane)
Solution 2: ##c_{311} = c_{322} = 0##, ##c_{312} = -c_{321} = 1##. (The difference of unlike partial derivatives of the x and y components of A, i.e. the z-component of the curl of A)
You don't really need to compute the conditions on ##c_{ijk}## for invariance under rotations in the x or y planes to see that only Solution 2 is invariant under general 3D rotations. Solution 1 is the ordinary 2D divergence for the values of A on points in the z-plane. It clearly won't be invariant under rotations about the x or y axes. You don't even need to repeat the derivation to see what solution 2 will look like for each set of rotations, you can just use the cyclic substitution property to see that ##c_{ijk}## has to be the Levi-Civita symbol (the coefficients of curl).