Standard deviation of several measurements with uncertainties I have two 2 hours of GPS data with a sampling rate of 1 Hz (7200 measurements). The data are given in the form $(X, X_\sigma, Y, Y_\sigma, Z, Z_\sigma)$, where $N_\sigma$ is the measurement uncertainty.
When I take the mean of all measurements (e.g. the average Z value of those two hours), what is its standard deviation? I can of course calculate the standard deviation from the Z values, but then I neglect the fact that there are known measurement uncertainties...
Edit: The data is all from the same station, and all coordinates are remeasured every second. Due to satellite constellations etc., every measurement has a different uncertainity. The purpose of my analysis is to find the displacement due to the an external event, (i.e. an earthquake). I would like to take the mean for 7200 measurements (2h) before the earthquake and another mean for 2h after the earthquake, and then calculate the resulting difference (in height for example). In order to specifiy the standard deviation of this difference, I need to know the standard deviation of the two means.
 A: I suspect that the previous responses to this question may be a bit off the mark.  It seems to me that what the original poster is really asking here could be rephrased as, "given a series of vector measurements: $$\vec{\theta}_{i} = \left( \begin{array}{c} X_{i} \\ Y_{i} \\ Z_{i} \end{array}\right)$$ with $i=1, 2, 3,...,7200$, and measurement covariance: $$C_{i} = \left( \begin{array}{ccc} X_{\sigma,i}^{2} & 0 & 0 \\ 0 & Y_{\sigma,i}^{2} & 0 \\ 0 & 0 & Z_{\sigma,i}^{2} \end{array} \right)$$ how would I correctly calculate the covariance-weighted mean for this series of vector measurements, and afterward, how would I correctly calculate its standard deviation?" The answer to this question can be found in a lot of textbooks specializing in statistics for the physical sciences.  One example that I like in particular is Frederick James, "Statistical Methods in Experimental Physics", 2nd edition, World Scientific, 2006, Section 11.5.2, "Combining independent estimates", pg. 323-324.  Another very good, but more introductory-level text, which describes the variance-weighted mean calculation for scalar values (as opposed to full vector quantities as presented above) is Philip R. Bevington and D. Keith Robinson, "Data Reduction and Error Analysis for the Physical Sciences", 3rd edition, McGraw-Hill, 2003, Section 4.1.x, "Weighting the Data--Nonuniform Uncertainties".  Because the poster's question happened to have a diagonalized covariance matrix in this case (i.e., all of the off-diagonal elements are zero), the problem is actually separable into three individual (i.e., X, Y, Z) scalar weighted mean problems, so the Bevington and Robinson analysis applies equally well here too.
In general, when responding to stackexchange.com questions, I don't normally find it useful to repackage long derivations that have already been presented before in numerous textbooks--if you want to truly understand the material, and understand why the answers look the way they do, then you really should just go and read the explanations which have already been published by the textbook authors.  With that in mind, I'll simply jump directly to re-stating the answers that others have already provided.  From Frederick James, setting $N=7200$, the weighted mean is: $$\vec{\theta}_{mean} = \left( \sum_{i=1}^{N} C_{i}^{-1} \right)^{-1} \left( \sum_{i=1}^{N} C_{i}^{-1} \vec{\theta_{i}} \right) $$ and the covariance of the weighted mean is: $$ C_{mean} = \left( \sum_{i=1}^{N} C_{i}^{-1} \right)^{-1} $$  This answer is completely general, and will be valid no matter what the form of the $C_{i}$, even for non-diagonal measurement covariance matrices.
Since it so happens that the measurement covariances are diagonal in this particular case, the Bevington and Robinson analysis can also be used to calculate variance-weighted means for the individual $X_{i}$, $Y_{i}$, and $Z_{i}$.  The form of the scalar answer is similar the form of the vector answer: $$ X_{mean} = \frac{\sum_{i=1}^{N} \frac{X_{i}}{X_{\sigma,i}^{2}}}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}} $$ and the variance is $$ X_{\sigma,mean}^{2} = \frac{1}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}} $$ or equivalently, $$ X_{\sigma,mean} = \sqrt{\frac{1}{\sum_{i=1}^{N} \frac{1}{X_{\sigma,i}^{2}}}} $$  and similarly for $Y_{mean}, Y_{\sigma, mean}$ and $Z_{mean}, Z_{\sigma, mean}$.  A brief wikipedia entry which also arrives at this same answer for the scalar-valued case is available here. 
A: This should be easily solved using bayesian inference. You know the measurement properties of the individual points with respect to their true value and want to infer the population mean and SD that generated the true values. This is a hierarchical model.
Rephrasing the problem (Bayes basics)
Note that whereas orthodox statistics give you a single mean, in the bayesian framework you get a distribution of credible values of the mean. E.g. the observations (1, 2, 3) with SDs (2, 2, 3) could have been generated by the Maximum Likelihood Estimate of 2 but also by a mean of 2.1 or 1.8, though slightly less likely (given the data) than the MLE. So in addition to the SD, we also infer the mean.
Another conceptual difference is that you have to define your knowledge state before making the observations. We call this priors. You might know in advance that a certain area was scanned and in a certain height range. The complete absence of knowledge would be to have uniform(-90, 90) degrees as the prior in X and Y and maybe uniform(0, 10000) meters on height (above the ocean, below the highest point on earth). You have to define priors distributions for all parameters that you want to estimate, i.e. get posterior distributions for. This is true for the standard deviation as well.
So rephrasing your problem, I assume that you want to infer credible values for three means (X.mean, Y.mean, X.mean) and three standard deviations (X.sd, Y.sd, X.sd) which could have generated your data.
The model
Using standard BUGS syntax (use WinBUGS, OpenBUGS, JAGS, stan or other packages to run this), your model would look something like this:
  model {
    # Set priors on population parameters
    X.mean ~ dunif(-90, 90)
    Y.mean ~ dunif(-90, 90)
    Z.mean ~ dunif(0, 10000)
    X.sd ~ dunif(0, 10)  # use something with better properties, i.e. Jeffreys prior.
    Y.sd ~ dunif(0, 10)
    Z.sd ~ dunif(0, 100)
    
    # Loop through data (or: set up plates)
    # assuming observed(x, sd(x), y, sd(y) z, sd(z)) = d[i, 1:6]
    for(i in 1:n.obs) {
      # The true value was generated from population parameters
      X[i] ~ dnorm(X.mean, X.sd^-2)  #^-2 converts from SD to precision
      Y[i] ~ dnorm(Y.mean, Y.sd^-2)
      Z[i] ~ dnorm(Z.mean, Z.sd^-2)

      # The observation was generated from the true value and a known measurement error
      d[i, 1] ~ dnorm(X[i], d[i, 2]^-2)  #^-2 converts from SD to precision
      d[i, 3] ~ dnorm(Y[i], d[i, 4]^-2)
      d[i, 5] ~ dnorm(Z[i], d[i, 6]^-2)
    }
  }

Naturally, you monitor the .mean and .sd parameters and use their posteriors for inference.
Simulation
I simulated some data like this:
# Simulate 500 data points
x = rnorm(500, -10, 5)  # mean -10, sd 5
y = rnorm(500, 20, 5)  # mean 20, sd 4
z = rnorm(500, 2000, 10)  # mean 2000, sd 10
d = cbind(x, 0.1, y, 0.1, z, 3)  # added constant measurement errors of 0.1 deg, 0.1 deg and 3 meters
n.obs = dim(d)[1]

Then ran the model using JAGS for 2000 iterations after a burnin of 500 iterations. Here's the result for X.sd.

The blue range indicates the 95% Highest Posterior Density or Credible interval (where you believe the parameter is after having observed the data. notice that an orthodox confidence interval does not give you this).
The red vertical line is the MLE estimate of the raw data. It is usually the case that the most likely parameter in Bayesian estimation is also the most likely (maximum likelihood) parameter in orthodox stats. But you should not care too much about the top of the posterior. The mean or median is better if you want to boil it down to a single number.
Notice that MLE/top is not at 5 because the data were randomly generated, not because of wrong stats.
Limitiations
This is a simple model which has several flaws currently.

*

*It doesn't handle the identity of -90 and 90 degrees. This can be done, however, by making some intermediate variable which shifts extreme values of estimated parameters into the (-90, 90) range.

*X, Y and Z are currently modeled as independent though they are probably correlated and this should be taken into account to get the most out of the data. It depends on whether the measurement device was moving (serial correlation and joint distribution of X, Y and Z will give you a lot of information) or standing still (independence is ok). I can expand the answer to approach this, if requested.

I should mention that there's a lot of literature on spatial Bayesian models which I am not knowledgeable about.
A: I first introduce some notation and set up the problem using the simple approach you mentioned. Then go further. I will use $\mathbf{z}$ to refer to the vector Z you gave.
Consider the following model, which lacks the explicit mention measurement error:
$\bar{Z} = \frac{\sum_{i=1}^n\mu_{Z} + \epsilon_i}{n}$, where $\bar{Z}$ is the estimated average value of $\mathbf{z}$, and $\mu_Z$ is the true average value of Z. Here, $\mathbf{\epsilon}$ is a vector of the errors in your data, and you expect that if your sample is large $\bar{Z}$ will converge to $\mu_Z$. If you simply take the observed $Z$ values and average them, you get $\bar{Z}$ and if you compute the sample standard deviation you get $\hat{\sigma}$, the estimate of the true population standard deviation $\sigma$. What if you would like to make use of some knowledge about the measurement error?
First, note that we can reformulate the initial model as:
$\mathbf{z} = \mathbf{1}\beta + \epsilon$, where $\mathbf{1}$ is a vector of ones, and $\beta$ will end up being $\bar{Z}$. Now this really looks like regression, but we are still basically just getting an estimate of $\mu_Z$. If we perform a regression like this, we will also get an estimate for the standard error of $\epsilon$, which is almost what we want - this is nothing but the standard error of $\mathbf{z}$ (but we still want to account for measurement error). 
We can augment our initial model to obtain a mixed effects model.
$\mathbf{z} = \mathbf{1}\beta + \mathbf{Q}u +\epsilon$, where $u$ is a vector of random effects, and $\mathbf{Q}$ is the regressor relating $\mathbf{z}$ to $u$.
As with any random effect, you will need to make an assumption about the distribution of $u$. Is it correct that $Z_\sigma$ is the distribution of the measurement error for $\mathbf{z}$? If yes, this can be used to provide the distribution of the random effects. Typically, software to perform basic mixed effects modeling will assume the random effects have a normal distribution (with mean 0...) and estimate the variance for you. Perhaps you can try this to test the concept. If you wish to use your prior information about the distribution of the measurement error, a Bayesian mixed effects model is in order. You can use R2OpenBUGS.
After estimating this model, the standard error you get for the residuals $\epsilon$ is the standard error you express interest in. 
Intuitively, the random effects component of the model is soaking up some of the variation that you can explain because you know there is measurement error. This allows you to obtain a more relevant estimate of the variation of $\epsilon$
See this paper for a deeper discussion on this approach of random effects to account for measurement error. Your situation is similar to the one the authors introduce for $\mathbf{D}$ and its measurement error corrupted version $\mathbf{W}$. The example in Section 4 may offer some insights into your situation.
As mentioned by whuber, you may wish to account for autocorrelation in your data. Using random effects will not solve that problem.
