I think your confusion is most concisely caught in this part of your question:
If I insist on using these coordinates to describe my space, I don't see how there can be any remaining gauge freedom.
Without getting into the argument that's emerged in comments to your question about whether your particular choice of coordinates makes sense, it is generally true that once you choose a set of coordinates, you do not have any gauge freedom left. The gauge freedom in general relativity is exactly the freedom to choose your coordinates. No more and no less.
I think that's most naturally viewed through an ADM decomposition to 3+1 dimensions. (The part of your question about gravitational waves also gets computed in this formalism, or a variant of it.) Assuming for simplicity at this point a vacuum solution, if $g_{\mu\nu}$ is the metric of your spacetime, Arnowitt, Deser, and Misner showed in the 1950s that if you take a foliation of spacelike hypersurfaces, you can decompose this into an induced 3-metric on the surfaces $\gamma_{ij}$ plus and "lapse" function $\alpha$ and a shift vector $\beta^i$. These are related by
$$ \left( \begin{array}{cc} g_{00} & g_{0j} \\ g_{i0} & g_{ij} \end{array} \right) = \left( \begin{array}{cc} \beta_k \beta^k - \alpha^2 & \beta_j \\ \beta_i & \gamma_{ij} \end{array} \right) $$
where the index on the shift is lowered by the 3-metric, $\beta_k = \gamma_{ik} \beta^i$. There's also a conjugate field $\pi^{ij}$ in this formalism that allow (coordinate) time derivatives $\partial_t \gamma_{ij}$ to be expressed as first-order-in-time PDEs.
I'm not going to write out all of the equations here because they are readily available, but making this split divides the Einstein equations into two types of equations: There are two "evolution" equations for $\gamma_{ij}$ and for $\pi^{ij}$, and there are two constraint equations. The constraints are a scalar constraint called the Hamiltonian constraint and a 3-vector constraint, known as the momentum constraint. (This is analogous, respectively, in EM to the two Maxwell equations that have time derivatives and the two that do not.)
In both the 3+1 formulation of GR and the Maxwell equations, one can prove that the "constraints propagate". That means that if the constraints are satisfied at any given time, they are satisfied at all times. See, for example, the appendix in Wald, which treats this in great detail.
Now note that the lapse and the shift are your "gauge" fields in the this formalism. They have no physical meaning. In particular, they don't appear in the Hamiltonian constraint or in the momentum constraint, which is analogous to the fact that the EM gauge fields don't appear in the constraint equations $\nabla \cdot E = 0$ or $\nabla \cdot B = 0$ (again assuming vacuum). Not appearing in the constraint equations means that if you have any ($\gamma_{ij}$, $\pi^{ij}$) pair that satisfy the constraint equations at any given time, i.e. they are "physical" at any time, you have can construct a solution to the full Einstein equations for any choice of lapse and shift by solving the evolution equations as an initial value problem, due to the fact that the constraints propagate, as mentioned above.
Secondly, note that the choice of lapse and shift are essentially equivalent to a choice of coordinates. You can find figures of this in many references, including MTW and Wald, but basically you can show that between two "nearby" surfaces $\Sigma_0$ and $\Sigma_1$, a point $x$ on $\Sigma_0$ will have the same spatial coordinate as a point $y$ on $\Sigma_1$ if $y \approx x + \alpha n + \beta$, where $n$ is normal to $\Sigma_0$ at $x$. You also have that the lapse is proportional to the proper time (rather than coordinate time) between slices.
So, in some sense you are right in your complaint. Once you fix the coordinates, you don't have gauge freedom left. When you construct your "personal LIGO" ala the end of your question, you will have effectively made your coordinate choice since the "natural" coordinates for your configuration will be determined by how you layout your detector's arms, giving you a rectangular coordinate system on your "spatial slices". In the weak-field region of Earth, that will be completely compatible with unit lapse and zero shift, i.e. the Minkowski space that you would have linearized around to compute your wave solution at the detector.
Even with this understanding, you can still find the effects of gauge. When you eventually detect a wave on your personal LIGO, it will have a polarization, which you'll probably decompose into "plus" and "cross" polarizations. If you rotate your detector, which is effectively a different gauge choice, you'll mix those differently for a given signal. This probably the simplest of multiple factors that you need to "match" between your theoretical calculation of the wave - where you have to fix the gauge - and your actual configuration of the detector on the ground. This is the other end of your mistake / confusion. When you do the theory, if you want to match to your detector, you not only need to choose a gauge, you need to choose a gauge that's compatible with how you're going to read the results from your detector. You cannot not just choose any gauge that you like and try to use the results from the calculation directly with your detector.