1 Introduction

2 Notation

We need to find out two functions, f_x(x,y) and f_y(x,y), to describe the deviation. x and y are the true longitude and latitude respectively. f_x(x,y) and f_y(x,y) are the deviation delta at point (x,y). In other words, if (x,y) is the true coordinate, then (x + f_x(x,y),y + f_y(x,y)) is the deviated coordinate. Since |f_x(x + d,y) -f_x(x,y)| is much smaller than d (and same for y) we can also use the two functions to calculate the true coordinate from a deviated coordinate. Let (x,y) be the deviated coordinate, (x - f_x(x,y),y - f_y(x,y)) is the true coordinate.

In the following analysis, we will use octave. All the data used in the analysis has already been posted in my previous article¹ .

3 The four 1D cases

We first solve the functions in one dimension by setting one parameter to be constant. I.e. instead of solving f_x(x,y) directly, we first solve f_x(x,c) and f_x(c,y), where c is a constant. For example, we will look at f_x(x,39.36) and f_x(110.16,y). Similarly, we will solve f_y(x,c) and f_y(c,y) before f_y(x,y). In total, there are four 1D cases.

3.1 f_x(const,y)

f_x(const,y) looks like to be the simplest, so we solve it first.

> data=load("-ascii", "iny-f-110.txt");
> y=data(:,2);
> dx=data(:,3);
> plot(y,dx) Figure 1(left)

Figure 1(left) shows the plot, f_x(110.16,y). The curve looks like a second order curve, but if we go that direction, it’s going to be very tough. It is a lot simpler if we measure the deviation in meters rather than degrees. Note that earth’s radius is 6378137 meters in Google Map. I didn’t know this until the experiment was all finished. However, this shouldn’t cause any error even if I set the radius to 1 meter, because in the end, meters is converted back to degrees and the error cancels.

> plot(y,dx.*cos(y./180.*pi).*6371000.*pi./180) Figure 1(right)

Figure 1(right) shows the plot where the deviation is measured in meters. It’s a linear function, thus we can do a simple linear regression.

The standard error, 0.14664 meter, is acceptable, because the input data’s resolution is about 0.6 meter² on average. The final function is

(1)

After trying out other constant longitudes, we find out that the two coefficients are different for different longitude. We have

(2)

where g₁(x) and g₂(x) are not yet known.

In the following analysis, we will measure x and y in degrees, while measure f_x(x,y) and f_y(x,y) in meters. It makes our life a lot easier.

3.2 f_y(x,const)

f_y(x,const) seems to be the second easiest one. It’s quite obvious that there are sinusoidal and linear components in the function. To find out the frequency of the sinusoidal components, we do a Fourier transform.

> plot(x,dy) Figure 2(left)
> plot(abs(fft(dy))) Figure 2(right)

Figure 2(left) shows the data and Figure 2(right) shows the data in frequency domain. The two peaks in the frequency domains are at 47 and 138, thus the frequencies are 46/45.575=1.0093 and 137/45.575=3.0060, where 45.575=119.23-73.655 is the range of x. Since the range of x is not multiple of the periods, we don’t get exact results, but after cutting the range, the exact results are 1 and 3. We then try to model it using b + kx + a₁sin(2πx) + a₂cos(2πx) + a₃cos(6πx) + a₄cos(6πx).

> plot(x,dy,x,dy-mx*(pinv(mx)*dy)) Figure 3

The standard error, 0.38114, is acceptable. We can also see this in Figure 3. The regression error (difference between the data and our function) is close to zero. Let’s see the coefficients:

We can combine sin and cos of the same frequency.

So the function is

After trying out other constant latitudes, we find out that the first two coefficients to be changing but the last 4 coefficients to be fixed. Thus, we know

3.3 f_x(x,const)

Similar to previous one, we use Fourier transform to find out the sinusoidal components first.

> plot(x,dx) Figure 4(left)
> plot(abs(fft(dx))) Figure 4(right)

Figure 4(left) shows the curve that we are going to fit. Figure 4(right) shows the frequency domain. The peaks are at 9, 24, 47 and 138, so the frequencies are 8/45.575=0.17553, 23/45.575=0.50466, 46/45.575=1.0093 and 137/45.575=3.0060. The exact frequencies are 1/6, 1/2, 1 and 3. Now, we try to fit the curve with the 4 sinusoidal components and a linear component.

The standard error 119.96 meters is not acceptable. To see it, we plot the difference in Figure 5. It seems that there are some low frequency components. We can deal with low frequency components with a bunch of polynomials.

It fits well with a 9-order polynomial. Note that we divide x by mean(x) to make sure the values in the matrix have low dynamic range so as to reduce numeric error. Mathematically, they don’t make a difference. To see how our polynomial works, we can plot the polynomial and compare with the data (see Figure 6).

Although the 9-order polynomial fits the data well (error is 0.79 meters), the function is ugly. After lots of tries³ , I figured out that the low frequency component actually consists of a quadratic curve and two sinusoidal wave of period 60 and 24. Using this model, the regression error is 0.39 meters and the function is simple and beautiful.

Now let’s merge sin and cos.

Finally, we have

After trying out other constant latitude, we found that only the first two coefficients are affected by latitude. The rest of the coefficients do not change. Thus, we have

3.4 f_y(const,y)

The process of finding f_y(const,y) is similar to f_x(x,const), except that the two high frequency sinusoidal components do not appear in f_y(const,y).

We have

4 The final 2D case: f_x(x,y) and f_y(x,y)

From Equation 6 and 2 we know that f_x(x,y) must be in the form of Equation 9. Similarly, from Equation 4 and 8 we know f_y(x,y) must be in the form of Equation 10. What’s left now is to determine the parameters.

Applying the previous linear regression method on the data file google-grid40.tsv, we can determine all the parameters easily. The octave commands are skipped. The resulting equations are shown below and their standard errors are 0.4373 and 0.43998 meters respectively, which are larger than the four 1D cases, but still below the input data’s resolution.

5 Evaluation

To check their correctness, we can compute a deviated coordinate and check whether it aligns with the deviated map. We now compute the deviated coordinate of the Monument to the People’s Heroes using its true coordinate (116.391667,39.903056).⁴ f_x(116.391667,39.903056) = 531.91 meters and f_y(116.391667,39.903056) = 155.39 meters. This means that the deviated coordinate is 531.91 meters to the east and 155.39 to the north. We then convert meters into degrees and adds them to the true coordinate.

The green arrow in Figure 8 shows the true coordinate in the Google Maps’ satellite image. Figure 9 shows the deviated coordinate in the Google Maps’ street map. The red balloon, which is the nearest landmark, should be ignored. The green arrows in both maps are slightly to the southeast of the monument center, but as long as they offset the same amount, it is as good.