trigdefs : [
    cos(phi) = x/r_2,
    sin(phi) = y/r_2,
    cos(theta) = z/r,
    sin(theta) = r_2/r
];

[x,y,z] = [r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)];
jacobian(rhs(%), [r,theta,phi]);
invert(%)$trigsimp(%)$ ratsimp(%);
subst(trigdefs,%)$
ratsimp(%);
drhp_def : transpose(matrix([dr, dtheta, dphi])) = % . transpose(matrix([dx,dy,dz]));

ds_def : (1 - r_s / r)*dr^2 + r^2 * dtheta^2 + r^2 * sin(theta)^2 * dphi^2;
subst(trigdefs, %)$
subst(makelist(lhs(drhp_def)[i,1]=rhs(drhp_def)[i,1],i,1,3), %);
subst([r_2 = sqrt(x^2 + y^2), r = sqrt(x^2 + y^2 + z^2)], %);
ratsimp(%);