In order to be able to obtain the (spin weighted) spherical harmonic modes of the strain $h_{lm}$ from two polarization modes $h_{+}$ and $h_{\times}$, you will need to know $h_{+}$ and $h_{\times}$ on all points of the celestial sphere of the sources, i.e. you'd need to observe the event from all possible directions.
This, of course, is not a situation that we would find ourselves in with gravitational waves observed "in the wild". We generally observe a gravitational wave event from only one direction. However, we may find ourselves in exactly this situation when doing numerical (or even semi-analytical) simulations of gravitational wave events.
So suppose we are given polarization modes $h_{+}(t,\theta,\phi)$ and $h_{\times}(t,\theta,\phi)$ on the entirety of future null infinity. (I stress again that $\phi$ and $\theta$ measure all directions from the source, not all directions from us on earth.) We can then recover the spherical harmonic modes $h_{lm}(t)$ as follows. First we construct the so-called complex strain
$$ h(t,\theta,\phi) = h_{+}(t,\theta,\phi) -i h_{\times}(t,\theta,\phi) .$$
As was explained in the linked question, this complex strain is related to the spherical harmonic modes as follows
$$ h(t,\theta,\phi) = \sum_{m=-\infty}^{\infty}\sum_{l=|m|}^{\infty} h_{lm}(t) {_{-2}Y_{lm}}(\theta,\phi),$$
where ${_{-2}Y_{lm}}(\theta,\phi)$ is a spin-weighted spherical harmonic of weight -2. These spin-weighted harmonics satisfy an orthogonality condition
$$ \int_{0}^{2\pi}\int_{0}^{\pi} {_{-2}Y_{l'm'}^{*}}(\theta,\phi){_{-2}Y_{lm}}(\theta,\phi) \sin(\theta)\,d\theta d\phi = \delta_{ll'}\delta_{mm'}.$$
Using this relation we can invert the expression for the mode expansion (similar to an Fourier transform),
$$ h_{lm}(t) = \int_{0}^{2\pi}\int_{0}^{\pi} {_{-2}Y_{lm}^{*}}(\theta,\phi) h(t,\theta,\phi)\,d\theta d\phi.$$
