There is a straightforward, disciplined way to analyze such data: linear regression. It can work because (a) the range of these points is a relatively small part of the world and (b) to a close approximation over such ranges, projection formulas are nearly linear as a function of latitude and longitude and--with more data--even their departures from linearity can be accommodated.
In this particular case the possibilities for investigation are limited by the small amount of data. We cannot hope to use more than four parameters per coordinate to fit these five points (a fifth parameter would result in a perfect fit and thereby be useless for predicting). We might therefore try models of the form
latitude = intercept1 + b1*X + b2*Y + b3*X*Y + error
longitude = intercept2 + c1*X + c2*Y + c3*X*Y + error
The inclusion of an X*Y term evaluates a likely form of nonlinearity.
Least-squares fits find all eight coefficients to be significant. They result in the following predicted values along with (95%) lower and upper prediction limits for those values:
Latitudes
X Y fit lwr upr
1 2635536 1326596 33.68766 33.68749 33.68784
2 2253740 2804378 37.75574 37.75518 37.75630
3 2561993 1430129 33.97425 33.97414 33.97437
4 2919984 1356273 33.75631 33.75609 33.75653
5 2919931 1355933 33.75538 33.75516 33.75560
Longitudes
X Y fit lwr upr
1 2635536 1326596 -117.9130 -117.9148 -117.9112
2 2253740 2804378 -119.1223 -119.1280 -119.1166
3 2561993 1430129 -118.1500 -118.1512 -118.1488
4 2919984 1356273 -116.9695 -116.9718 -116.9673
5 2919931 1355933 -116.9697 -116.9720 -116.9675
These prediction intervals (from lwr to upr) suggest uncertainty in the second digit after the decimal point--an error of up to a half kilometer or so. There is this much possible error because many of the points are extrapolations: they are relatively far from the points for which lat-lon are known. But if you can live with this much uncertainty, you might not need to guess which of the thousands of possible projections might have been used.
This is the complete R code that performed the analysis. First, read in the data:
D <- data.frame(X=c(2658206, 2363687, 2101778, 2095540, 2699433),
Y=c(1371867, 1478454, 2290419, 2274891, 1470349),
Lat=c(33.8105204, 34.11562481, 36.35547224, 36.31321472, 34.07804599),
Lon=c(-117.834935, -118.8053659, -119.6548331, -119.6764159, -117.691484))
E <- data.frame(X=c(2635536, 2253740, 2561993, 2919984, 2919931),
Y=c(1326596, 2804378, 1430129, 1356273, 1355933))
Next, fit the models:
fit.lat <- lm(Lat ~ X + Y + X*Y, D)
fit.lon <- lm(Lon ~ X + Y + X*Y, D)
You can see their coefficients and assess their significance with the summary command. For instance, the latitude regression produces this output:
summary(fit.lat)
Residuals:
1 2 3 4 5
5.155e-07 1.202e-07 5.533e-06 -5.578e-06 -5.908e-07
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.012e+01 5.640e-04 53404.26 1.19e-05 ***
X -2.445e-08 2.659e-10 -91.96 0.00692 **
Y 2.774e-06 3.894e-10 7124.37 8.94e-05 ***
X:Y -1.451e-14 1.854e-16 -78.22 0.00814 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.896e-06 on 1 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 3.522e+10 on 3 and 1 DF, p-value: 3.917e-06
The residuals, which do not exceed 1e-06, show the fit is good to within one meter for these data points. The summary of the longitude fit is not quite so good, but all residuals are within five meters or so.
Incidentally, the intercept terms are so close to 0 in each case that we ought to suspect that whatever projection was used is for the entire world. Such a projection, when applied to a set of points covering such a small area, can be expected not to be highly accurate. The coefficients themselves tell us these coordinates must be in feet, not meters. The small coefficients of X:Y suggest we might not need the nonlinear "X*Y" term in the model--but including it appreciably narrows the prediction intervals.
Finally, here is the code to output the predictions and the prediction intervals:
cbind(E, predict(fit.lat, newdata=E, interval="prediction"))
cbind(E, predict(fit.lon, newdata=E, interval="prediction"))
Although many GISes can perform the same calculations, they typically do not offer the flexibility of a statistical package nor the additional power of seeing the prediction intervals as well as the predictions themselves.
