Solving the field membership problem using Grobner bases
In theory, Groebner basis software can do this. I don't know whether there is a standard package for this purpose, or what practical issues may arise. Here is the theory:
Your goal is the following. You are given rational functions $r_1$, $r_2$, ..., $r_n$ in $x$ and $y$. You would like know whether there are polynomials $f(t_1, \ldots, t_n)$ and $g(t_1,\ldots, t_n)$ such that $f(r_1, \ldots, r_n) - g(r_1, \ldots, r_n) x =0$, and similarly for $y$ in place of $x$.
Our strategy will be the following: Map $k[t_1, \ldots, t_n, u] \to k(x,y)$ by $t_i \mapsto r_i$, $u \mapsto x$. Let $I$ be the kernel of this map.
We will test whether $I$ contains an element whose degree in $u$ is $1$. Actually we need to do slightly more than that: we need to find an element $f u - g$ where $f$ and $g$ are not in $I$, but the Grobner technology will do that for us.
The problem breaks into two parts
Getting the ideal $I$ If the $r_i$ are themselves polynomials in $x$ and $y$, you can see a worked example of doing this with Grobner bases in the Macualay documentation here. Macaulay also provides a function kernel which does this; I don't know what the advantages are of using kernel or doing it by hand.
If $r_i=p_i/q_i$, my first idea is to define the ring $S = k[x,y,s_1, \ldots, s_n]/(q_i(x,y) s_i -1)$, and map $k[t_1, \ldots, t_n]$ to $S$ by $t_i \mapsto p_i(x,y) s_i$, and adapt the Grobner basis methods to this case. If you are working in Macaulay specifically, then Macaualy has a class FractionField, with Ring as an ancestor class, meaning you should be able to directly talk about maps $k[t_1, \ldots, t_n,u] \to k(x,y)$, and not just about maps to quotients of polynomial rings. However, the documentation of FractionField warns that using it is a performance hit; so you might still be better off only inverting the elements you need to invert.
Testing whether $I$ contains an element of degree $1$ in $u$ The way to do this in theory is to compute a Grobner basis of $I$ in lex order, putting the variable $u$ first. If the answer contains an element $fu-g$ of degree $1$, then this is your answer; if not, then $x$ cannot be written as a rational expression in the $r_i$. Note that $f$ and $g$ will not be in $I$ as, if they were, $fu-g$ would be eliminated from the Grobner basis as redundant.
But lex term orders are notoriously slow, so there might be a better way to do this in practice. Hopefully, someone will show up soon who understands these issues.
Between David's answer and much staring at the paper I linked to above, I think I've mostly figured out what's going on there so I thought I'd post a summary of their method here as well in case anyone else had a similar problem. I'm focusing here on $how$ to use their algorithm, not $why$ it works: I found it was a lot easier to follow the theory once I knew where I was going!
Their setup is: given a field $k(x_1,\dots, x_n) = k(\mathbf{x})$ and a set of elements $\mathbf{g} = \{g_1, \dots, g_m\}\subset k(\mathbf{x})$, we can form the field $k(\mathbf{g})$ generated by the elements of $\mathbf{g}$ and then ask: given $f \in k(\mathbf{x})$, do we have $f \in k(\mathbf{g})$ and if so, how can we write $f$ in terms of the generators? The authors' aim is to answer this without resorting to using the lex ordering, which as mentioned above can be very slow.
So let $g_i = n_i/d_i$, and introduce a new variable $T_i$ for each $g_i$. Collect together all of the prime divisors of each of the denominators $d_i$ into a set $\{p_1, \dots, p_r\}$ (eliminate any repeats) and introduce a new variable $Y_i$ for each $p_i$.
In the polynomial ring $k[Y_1, \dots, Y_r, x_1, \dots, x_n, T_1, \dots, T_m]$ define the ideal $I=\langle n_1-T_1d_1, \dots, n_m-T_md_m, Y_1p_1-1,\dots,Y_rp_r-1\rangle$ and compute a Groebner basis of $I$. The only restriction on ordering is that you should have $\textbf{Y} >> \textbf{x} >> \textbf{T}$, within each block they can be ordered however you like.
Eliminate all terms involving any $Y_i$ from the Groebner basis; call the resulting set $G$. Let $H$ be the subset of $G$ involving terms only in the $T_i$: this may be empty. (This step seems to deal with syzygies between the generators, giving you the simplest result at the end and avoiding dividing by zero at any point.)
Let $S=\langle H \rangle$ and define $L = Q(k[\mathbf{T}]/S)$, so if $H$ was empty this is just $k(\mathbf{T})$. Let $\pi$ be the natural embedding $k[\mathbf{T}][\mathbf{x}] \hookrightarrow L[\mathbf{x}]$ and $G' := \pi(G)$.
If $f=n_f/d_f \in k(\mathbf{x})$ is the element of interest then introduce a new variable $A$ (so we're now in the ring $L[\mathbf{x},A]$) and reduce the polynomial $n_f-Ad_f$ to normal form wrt the Groebner basis $G'$: you should obtain an expression $N-AD$.
Optionally at this point, you can reduce $N$ and $D$ to their simplest form by reducing their numerators wrt $\pi(H)$; recall that $N,D \in L[\mathbf{x}]$ where $L$ is a localization of $k[\mathbf{T}]$ so they could be rational functions in the $T_i$.
If $D = 0$ then $f \not\in k(\mathbf{g})$. If $D \neq 0$ then solve $N-AD = 0$, i.e. $A=N/D$. If $N/D \not\in k(\mathbf{T})$ then $f \not\in k(\mathbf{g})$, otherwise $f=N(g_1, \dots, g_m)/D(g_1,\dots,g_m)$.
Whew! I hope that's fairly clear if you just want to apply the algorithm; if you want to understand why each step works, the paper I linked to in the question has all the technical details. I suspect it's essentially the same approach as David's answer, with a few tweaks to avoid necessarily using the lex ordering and to deal with relations between the set of generators.