Let's break down the operations required:
First, solve for theta and r in terms of x and y:
theta_r = solve([x == x0 + r*cos(theta), y == y0 + r*sin(theta)], [theta, r]);
theta_r will be a struct with possible substitutions for theta and r as its properties.
Second, choose the substitutions (in this case, the first ones):
theta_r_subs = [theta_r.theta(1), theta_r.r(1)];
Third, perform the substitution:
fxy = subs(f, [theta, r], theta_r_subs);
df_x = diff(fxy, x);
df_y = diff(fxy, y);
Further solutions can be obtained by choosing different substitutions.