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.
Using the standard rules for error propagation, as mentioned on wikipedia, for each element $i$, the sum can be split into two parts: $a = y_i$ and $b = \sum y - y_i$. Then, ignoring cross-correlation between the elements of $y$:
$$
f(a,b) = a / (a+b) \\
\sigma^2_f = |\frac{df(a,b)}{d a}|^2 \sigma_a^2 + |\frac{df(a,b)}{d b}|^2 \sigma_b^2 \\
\sigma^2_f = (\frac{1}{a+b} - \frac{a}{(a+b)^2})^2 \sigma_a^2 + (\frac{a}{(a+b)^2})^2 \sigma_b^2
$$
Best Answer
This answer will assume your errors are standard deviations.
If you have a data set $x_1,\ldots,x_n$, then we can define the discrete mean and variance as $$\langle{x}\rangle\equiv\frac{1}{n}\sum_ix_i \,,\,\hat{\sigma}^2\equiv\langle{x^2}\rangle-\langle{x}\rangle^2$$ which means $$\langle{x}\rangle{n}=\sum_ix_i \,,\, \langle{x^2}\rangle{n}=\sum_ix_i^2$$
Now imagine your data is really giving $$X_k\pm{\Delta}X_k=\langle{x}\rangle_k\pm\hat{\sigma}_k$$ where the "subsample size" $n_k$ is unknown.
Then we have $$\langle{x}\rangle_k=X_k \,,\, \langle{x^2}\rangle_k={\Delta}X_k^2-\langle{x}\rangle_k^2$$ And the aggregate statistics can be computed via $$\langle{X}\rangle=\frac{1}{N}\sum_kX_kn_k \,,\, \langle{X^2}\rangle=\frac{1}{N}\sum_k(X_k^2+\Delta{X}_k^2)n_k \,,\, N=\sum_kn_k$$ from which you can compute $$\hat{\sigma}_X^2=\langle{X^2}\rangle-\langle{X}\rangle^2$$
For simplicity you could take $n_k=1$.