This is more a collection of really simple comments than an actual answer, but you mention you want to know where to start to look, so maybe this helps a bit and may guide you to ask more specific questions in the future if you want to continue with it. I'm by no means an expert.
Going to assume the problem statement is minimize $F(q) = \int_L d(p,q) dp$, where $L$ is a subset of the sphere $\mathbb{S}^2$. If $L$ is a curve then I'm going to interpret the integral as the integral over the curve, and if it's a collection of points it's going to be the sum of the distance over all points.
There are some obstacles to this problem:
The land $L$ is what changes the whole computation. It's trivial to do the calculation for some examples of $L$, if we take $L=\{p_0\}$ then it's clear $q$ should be $p_0$. If $L=\mathbb{S}^2$ then all the points are valid as $q$ because $F$ is constant. But for your actual statement, you want $L$ to be somewhat similar to the surface land on Earth. That's difficult.
The distance $d$ is the great-circle distance. It's important to note that this distance has some weird behaviour that will be translated to your problem in some weird scenarios. With weird I mean lack of differentiability over antipodal points. A good thing is that we can compute it, check the link to wikipedia for some formula examples, which are already in spherical coordinate system. This is a benefit of working in a sphere instead of working over other manifold, which means the general problem is actually even harder.
Some starting points I would do:
Theory: You want to analyze this problem mathematically, for academic purposes or just because. Then define $F$ over spherical coordinates as $F(\phi, \lambda)$. Notice that $d((\phi_1, \lambda_1), (\phi_2, \lambda_2))$ has some differentiability problems when we reach the antipodal point. It would be nice to have something like 'Every $(\phi_0, \lambda_0)$ that minimizes $F$ must hold that $\frac{\partial F}{\partial \phi} (\phi_0, \lambda_0) = 0$', like we have in real functions and other settings. In mathematics we always try to connect extrema points with derivatives because it gives us more information. But $d$ having that differentiability problem doesn't let us continue. In the real line, the distance which is the absolute value of the difference also has this problem. The way we can fix that is by minimizing $|x-y|^2$. The natural question here is: does this squaring method work for the distance over the sphere? It wouldn't be the expectation anymore though. Can we fix it by other methods and arrive to a formula that $q$ should hold? Notice that if we can fix it we could apply the differentiation under integral sign to continue.
Practice: If you want a specific solution to your specific question (actual Earth, actual expectation, no need to study it mathematically), then I would go to numerical methods to find the actual value. Find the Earth projection you like the most with an actual map with land of one color and water of another. Find the formula of said projection to the sphere, and program a little function with that image that given a point on Earth it gives you "Land" if it's land or "Water" if it's water. You will also need to discretize the sphere, so that the integral can be approximated as the average over points that are in land multiplied by area of land. After that, you have the function $F$ that you can compute with your numerical methods and formula from wikipedia, and a mesh of points where you can calculate it, you just need to find the $q$ that minimizes it. If you don't care about accuracy, then discretizing in $100$ subintervals both axis (with ~$400$km between each point) and just greedily find the minimum would be enough. If you have more subintervals you'd need to be smarter about how to find this minimum.
