Inferring lat/lon based on known location coordinates

I have a point dataset with a set of coordinates that is in an unknown projection. For some of these points I know the unprojected latitude and longitude. Based on these known data points, would it be possible to infer what the latitude and longitudes are for the unknown points? Is there a formula that can do this?

Here is a sample of the data - I would like to get the Lat/Lon of the first five based on what I know for the last 5 data points:

PROJECTED_X PROJECTED_Y Latitude        Longitude
2635536     1326596
2253740     2804378
2561993     1430129
2919984     1356273
2919931     1355933
2658206     1371867     33.8105204      -117.8349355
2363687     1478454     34.11562481     -118.8053659
2101778     2290419     36.35547224     -119.6548331
2095540     2274891     36.31321472     -119.6764159
2699433     1470349     34.07804599     -117.691484
• Projections are based on complex formulas... I'm afraid you'll have to find out your original projection to get your Lat/Lon. I would be glad to be contradicted though. – ArMoraer Jan 15 '16 at 23:50
• Your data seems almost linear (x as longitude, y as latitude), but not quite. Do you know anything about where this data came from? – barrycarter Jan 16 '16 at 17:50
• Unfortunately no - I am not sure which system this data comes from. – sriramn Jan 17 '16 at 3:05
• In arcgis there is spatial adjustment tool you can use to move points in unknown system to degrees. – FelixIP Jan 17 '16 at 3:54

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.

• This is excellent, thank you! I was thinking along similar lines but was unsure of how to assess the accuracy of the prediction. An error of 1 to 5 m is acceptable for my empirical application. – sriramn Jan 21 '16 at 20:07
• Please note that the 1-5 m error applies only to the fitted data--that is, the points you already know. The prediction errors are what matter for assessing how close you might come to guessing the coordinates of the other points--and those could be off by around half a kilometer. – whuber Jan 21 '16 at 20:08

Yes, it's possible.

Convert the known coordinates from lat-long into various coordinate systems and see what matches.

Looking at a few of your lat-long values, it seems that your data is in California, so my guess would be California State Plane; this may be in meters or feet, but the conversion will tell you quickly.