The direct calculation of the derivatives isn't that hard. But one can also quickly see the values of the Riemann tensor for a sphere – and similar simple shapes – by using the definition of the Riemann tensor via the parallel transport of vectors.
$$\delta V^\alpha = R_{\alpha\beta\gamma\delta}V^\beta d\Sigma^{\gamma\delta}$$
Around a point of the sphere $S^d$, the transport around an area given by $d\Sigma^{\gamma\delta}$ for fixed values of the indices (locally orthonormal basis) allows you to see that all of this is happening on an $S^2$ only. The other dimensions are unaffected. That's why you get $R_{\alpha\beta\gamma\delta}$ equal to $(1/a^2)$ times $g_{\alpha\gamma}g_{\beta\delta}-g_{\alpha\delta}g_{\beta\gamma}$. Effectively, the antisymmetrized pair of indices $\alpha\beta$ has to be the same as the pair $\gamma\delta$. I didn't assume anything special about the point; all points on a sphere are equally good by a symmetry. So the Ansatz for the Riemann tensor has to hold everywhere.
Note that it's multiplied by $1/a^2$, the inverse squared radius of the sphere. Many of your formulae omit it; moreover, you are using a confusing symbol $R$ for the radius which looks like the Ricci scalar – a different thing.
In $d=2$, the Riemann tensor only has one independent component and the formula for the Riemann tensor in terms of the metric tensor above actually holds for any surface if $1/a^2$ is replaced by $R/2$. Note that the two-sphere has (Ricci scalar) $R=2/a^2$. Also, the Ricci tensor is $R_{ij}=Rg_{ij}/2$ in $d=2$ so that the vacuum Einstein equations are obeyed identically.
For $d=3$, the Riemann tensor has 3 independent components, just like the Ricci tensor, so the Riemann tensor may be written in terms of the Ricci tensor. That's not true for higher dimensions, either.
In practice, given a stress-energy tensor $T_{\mu\nu}$, we may attempt to find solutions to the Einstein field equations using perturbation theory. The basic idea is to expand around a known solution $g_{\mu\nu}$ by a perturbation $h_{\mu\nu}$. In the case of a flat background,
$$\delta G_{\mu\nu} = 8\pi G \delta T_{\mu\nu} = \partial_\mu \partial_\nu h - \partial_\mu \partial_\alpha h^\alpha_\nu -\partial_\nu \partial_\alpha h^\alpha_\mu + \partial_\alpha \partial^\alpha h_{\mu\nu} - \eta_{\mu\nu} \partial_\alpha \partial^\alpha h + \eta_{\mu\nu} \partial_\alpha \partial_\beta h^{\alpha \beta}.$$
In some cases, one may solve the equations exactly, or more typically, employ numerical methods. Another alternative approach when faced with a stress-energy tensor is to try to determine the symmetries the metric may have, and then plug in an ansatz for the metric to yield a set of differential equations which may be more tractable, analytically and numerically.
Yet another alternative is to generate new solutions from old ones through the introduction of pseudopotentials, a method due to Harrison, Eastbrook and Wahlquist, called the method of prolongation structures.
There are several other methods such as those that rely on Lie point symmetries of differential equations. There is also a Backlund transformation, which relies on identifying a simpler differential equation whose solution satisfies a condition involving the solution to the harder problem.
These methods are too involved to present here and require a significant background. They are explained in Exact Solutions to the Einstein Field Equations by H. Stephani et al.
Addressing your other question, if given a metric $g_{\mu\nu}$, of course one can compute $T_{\mu\nu}$ through the Einstein field equations, you just plug it in and tediously compute all the curvature tensors. There are strictly speaking probably some requirements on the functions in $g_{\mu\nu}$, but you can get away with most things, even distributions, such as a delta function, which may lead to a stress-energy tensor describing a brane.
Best Answer
The key to proving item 2 is to express the metric in Riemann normal coordinates, which is usually what is meant when you say you are working in a locally inertial coordinate system. In these coordinates, the metric is equal to the Minkowski metric at a point, the first derivatives of the metric at the point vanish, and the second derivatives of the metric are given by the Riemann tensor at that point. The explicit form of the metric components are then (see e.g. this document)
$$g_{\mu\nu} = \eta_{\mu\nu}-\frac13R_{\mu\alpha\nu\beta}x^\alpha x^\beta + ...$$
where the dots represent higher order corrections in the coordinate distance from the origin, $x^\alpha = 0$.
We need to compute the volume of a sphere of coordinate radius $r$. For this we need the spatial metric, which is $h_{\alpha \beta} \equiv g_{\alpha\beta}+u_\alpha u_\beta$, and $u^\alpha$ is tangent to the inertial observer, so points in the time direction. The spatial volume element comes from the determinant of $h_{ij}$ as a spatial tensor ($i,j$ are only spatial indices). We have
$$h_{ij} = \delta_{ij} -\frac13R_{i\mu j\nu}x^\mu x^\nu+...$$
and the first order correction to the determinant just adds the trace of this tensor,
$$\sqrt{h} = 1 + \frac12\left(-\frac13\delta^{kl}R_{kilj}x^i x^j\right). $$
It will be useful to work with spacetime indices in a moment, where the background spatial metric is given by $\delta_{\mu\nu} = \eta_{\mu\nu} +u_\mu u_\nu$, and its inverse is $\delta^{\mu\nu} = \eta^{\mu\nu}+u^\mu u^\nu$. Now the volume of the sphere is simply
$$V = \int d^3x \sqrt{h} = \int d^3x\left(1-\frac16 \delta^{\mu\nu}x^\alpha x^\beta R_{\mu\alpha\nu\beta}\right).$$
(the limit of integration is over a coordinate sphere centered at the origin).
The first term will give the flat space volume of the sphere, so we need to compute the second term to get what Feynman is calling the spatial curvature of space. Remember that the Riemann tensor is taken to be constant since it is evaluated at the origin. Also, when integrated over a spherical region, only the trace of $x^\alpha x^\beta$ contributes, the other parts canceling out, so we can replace $x^\alpha x^\beta \rightarrow \frac13 r^2 \delta^{\alpha\beta}$. So the integral we are computing becomes
$$\Delta V = -\frac16\frac{4\pi}{3}\delta^{\mu\nu}\delta^{\alpha\beta}R_{\mu\alpha\nu\beta}\int_0^{r_s} r^4 dr = -\frac{2\pi}{45} r^5\delta^{\mu\nu}\delta^{\alpha\beta}R_{\mu\alpha\nu\beta}.$$
The numerical coefficient is not important, we only care about the dependence on the Riemann tensor. Re-writing the $\delta$'s in terms of the background metric $\eta^{\mu\nu}$ and $u^\alpha$, we get
$$\delta^{\mu\nu}\delta^{\alpha\beta}R_{\mu\alpha\nu\beta}=(\eta^{\mu\nu}+u^\mu u^\nu)(\eta^{\alpha\beta}+u^\alpha u^\beta)R_{\mu\alpha\nu\beta} = R + 2 R_{\mu\nu}u^\mu u^\nu,$$
where $R$ is the Ricci scalar at the origin. Now we can easily check that this is proportional to the $uu$-component of the Einstein tensor (using $u^\mu u^\nu g_{\mu\nu} = -1$),
$$G_{\mu\nu}u^\mu u^\nu = \frac12(2 R_{\mu\nu}u^\mu u^\nu + R) \checkmark$$
Then from the rest of your arguments, we arrive at Feynman's conclusion: the energy density is proportional to the spatial curvature in all locally inertial frames.