Is there a better way to find the polynomial equation for this curve?

Let $\varepsilon=e^{2\pi i/3}$. Your equation will be $$ \prod_{k=0}^2 \prod_{\ell=0}^2 \Big( \varepsilon^k x^{1/3} + \varepsilon^\ell y^{1/3} + \varepsilon^{k+\ell} (xy)^{1/3} -1 \Big) =0. $$

Let $K = \mathbf{Q}(\zeta_3)$. You could just directly take the norm of the element $x^{1/3} + y^{1/3} + x^{1/3}y^{1/3} - 1$ from $K(x^{1/3}, y^{1/3})$ down to $K(x, y)$. This norm is computed as the product of 9 Galois conjugates of your element and will be an element of $\mathbf{Q}(x, y)$. Should be very easy to do in any computer algebra system which can work with 3rd roots of 1.

But this is just exactly the same as what you suggest!

Two automated approaches. With Gröbner basis: see the first entry of this result. With resultants: first compute this and then this.