Given a set of curved geodesics in 2D space, is there a way to determine the geometrical shape of said space?

The quantity you would like to obtain is the metric, which completely characterizes the shape of a space (Riemannian manifold, to be more precise). Unfortunately, knowing even every geodesic on a Riemannian manifold does not uniquely determine the metric. See, for example, this discussion: https://mathoverflow.net/questions/132244/can-one-recover-a-metric-from-geodesics.

As explained there, $\mathbb{R}^4$ with either Euclidean or Minkowski metric produces the same set of geodesics, namely straight lines, so in general it cannot be possible to obtain the metric from just the geodesics. Other notions of shape, such as the curvature tensor and sectional curvature, are defined in terms of this metric.

This problem is a case of obtaining local properties from global ones. The metric is a local quantity; it depends only on the chosen point of the manifold and a local neighborhood around it. Geodesics, on the other hand, connect points on the manifold which may not even lie on the same coordinate chart, appealing to a larger global structure of the manifold.

There are a variety of theorems dealing with the relationship between local and global properties in Riemannian geometry, a prolific example of which is the Hopf-Rinow theorem. This tells you that a Riemannian manifold is a complete metric space if and only if it is geodesically complete - that is, if at any point you can extend a geodesic infinitely far in any direction, the metric on your space is such that the manifold is complete. We can deduce a few interesting properties like this, but we cannot fully determine the metric from information about the geodesics.

From the set of geodesic, as the previous answers, the shape of the space-time can be partially determined. However knowing a bit more of information, the metric can be fully delimited in a neighborhood of every point. I think is a useful calculation, because represents how we can, as observers inside of the spacetime determine the nature of it through experiments.

The argument is an sketch of the one given in Sec 3.2 in "Large Scale Structure of the Space-Time" by Hawking and Ellis"

Given the null-vectors of the spacetime, the functional form of the metric is determined by local causality and the material content.

• Part 1: Geodesics and null-vector relation

Consider an observer in the spacetime at a point $p$. The observer can throw test particles that will move under non-spacelike geodesics. The tangent vector of the geodesic is an element of $T_p$. Throwing enough test particles following different geodesics (which is equivalent to know all the time-like geodesics that pass through $p$) we can determine the null cone.

In simple words, throwing particles from $p$ and seeing which points of the manifold can be reached, the null cone can be determined as the boundary of such hypersurface.

• Part 2: The null cone determine the functional form of the metric up to a conformal factor.

Consider known all the vectors of the null cone, as well as the timelike vectors (i.e. we can distinguish which geodesic are causal in our spacetime). Then every vector of the spacetime that is not null nor timelike, must be spacelike.

Let $\mathbf{T}$ be a timelike vector, and $\mathbf{S}$ a spacelike vector. Then there exist two values of $\lambda\in\mathbb{R}-\{0\}$ for which $\mathbf{T}+\lambda\mathbf{S}$ is null, so

$$ 0=\mathbf{g}(\mathbf{T}+\lambda\mathbf{S},\mathbf{T}+\lambda\mathbf{S})=\mathbf{g}(\mathbf{T},\mathbf{T})+2\lambda\mathbf{g}(\mathbf{T},\mathbf{S})+\lambda^2\mathbf{g}(\mathbf{S},\mathbf{S}) $$

This is a polynomial on $\lambda$ for which the roots $\lambda_1,\lambda_2$ are known (since we know all the vectors, and their character, we can determine for a given pair $\mathbf{T},\mathbf{S}$ which $\lambda$ makes $\mathbf{T}+\lambda\mathbf{S}$ null), then is true that:

$$ \mathbf{g}(\mathbf{T},\mathbf{T})+2\lambda\mathbf{g}(\mathbf{T},\mathbf{S})+\lambda^2\mathbf{g}(\mathbf{S},\mathbf{S})=\mathbf{g}(\mathbf{S},\mathbf{S})(\lambda-\lambda_1)(\lambda-\lambda_2)\Rightarrow \lambda_1\lambda_2=\frac{\mathbf{g}(\mathbf{T},\mathbf{T})}{\mathbf{g}(\mathbf{S},\mathbf{S})} $$

So the ratio between the norm of a timelike and spacelike vector can be found knowing the null cone.

Now let $\mathbf{W},\mathbf{Z}$ be two non-null vectors, then

$$ \mathbf{g}(\mathbf{W},\mathbf{Z})=\frac{1}{2}\left(\mathbf{g}(\mathbf{W},\mathbf{W})+\mathbf{g}(\mathbf{Z},\mathbf{Z})-\mathbf{g}(\mathbf{W+Z},\mathbf{W+Z})\right) $$

Each of the terms on the RHS can be connected to $\mathbf{g}(\mathbf{S},\mathbf{S})$ using different values of $\lambda_1,\lambda_2$. Sor for every pair $\mathbf{W},\mathbf{Z}$, the value of $\mathbf{g}(\mathbf{W},\mathbf{Z})$ is known up to a factor $\mathbf{g}(\mathbf{S},\mathbf{S})$.

• Part 3: The material content determine the conformal factor up to a measuring gauge.

For now we have that $\mathbf{\hat{g}}=\Omega^2\mathbf{g}$ where $\mathbf{g}$ is known.

Let the energy momentum tensor for the material fields be $T^{ab}$, satisfying $\nabla_aT^{ab}=0$. Since the spacetime must be locally Minkowsky (equivalent to take normal coordinates), there is a neighbourhood of $p$ in which we can define "almost killing vectors", taking the killing vectors of the minkowsky spacetime $K_a$. Since $K_aT^{ab}$ is a conserved current in Minkowsky, it will be almost conserved in our neighbourhood, in the sense that the first approximation vanishes. In particular that means that energy and momentum conservation hold approximatelly in the neighbourhood of $p$.

Given the timelike geodesic with respect to the metric $\mathbf{g}$ trajectory of a particle $\gamma(t)$ with tangent vector $\mathbf{K}=\partial_t$, the geodesic equation reads:

$$ K^{\left[b\right.}\frac{\mathbf{\hat{D}}}{\partial_t}K^{\left.a\right]}=K^{\left[b\right.}\frac{\mathbf{D}}{\partial_t}K^{\left.a\right]}-(K^cK^d\hat{g}_{cd})K^{\left[b\right.}g^{\left.a\right]e}\partial_e(\log\Omega) $$

Since $\gamma$ is a geodesic with respect $\mathbf{g}$, the first term vanishes. By considering another curve $\gamma^\prime$ whose tangent vector is not paralell to $\mathbf{K}$ $\Omega$ can be found up to a constant factor. This constant factor correspond to an arbitrary normalization (i.e. choosing a measure of time).

Not sure if this is exactly what you're looking for, but Synge's Theorem is a classic result in Riemannian geometry relating curvature to topology. In the proof of the theorem, one essentially performs a stability analysis of closed geodesics, which relates the properties of geodesics to the curvature of the manifold. This involves taking the second variation of the arc length functional, which unsurprisingly comes out with a term proportional to the curvature tensor. By examining how the closed geodesics behave under small variations, one can draw conclusions about the global properties of the space (i.e. its topology). For an example of how this works related to your immediate question on surfaces, consider the closed geodesics of $S^2$. These are great circles. Any of these geodesics can be shrunk by a small variation, for example by shifting them slightly towards one of the poles. By repeating this process, one can then shrink the geodesic to a point. This is clearly related to the fact that $S^2$ is simply connected. Hope this helped.