A Simple Heuristic for Random Error in Possibly Google-able Coordinates

Nov. 21, 2021

For my project mapping income in China during the 1960s with Oliver Kim, we had to decide whether or not to trust the provided coordinates of some of our village ground truth data. We figured villages are pretty sticky in their location over time, so it might make sense to just Google them and use those coordinates. A first pass showed that this isn't so simple and there are lots of villages with the same or similar-enough names in China to foil a naive approach. After supplying additional information, such as the province, we were able to get some results that appeared reasonable but were wary.

To decide whether or not we should stick with the provided or Google-derived coordinates, I devised a simple heuristic based on a test whose null hypethesis is that the Google coordinates are correct and the provided coordinates are noisy in a random direction (that is independent of the distance). The heuristic statistic is

$T=\sum\limits_{i=1}^N [(ProvidedCoords_i \notin County_i) - (\textrm{share of circle with radius }d(ProvidedCoords_i, GoogleCoords_i)\textrm{'s circumference}\notin County_i)]$

and when it is close to zero, we fail to reject the null. The intuition is as follows: if the provided coordinates are noisy in a random direction, then conditional on how far they're displaced they should appear outside of some boundary as often as if they were randomly placed in a different direction (in expectation). That latter number is the share of the circle centered on the Google coordinates with radius of the distance between the pairs of coordinates that lies outside the boundary. I illustrate the test below using Chinese provinces as the boundaries.

import pandas as pd
import geopandas as gpd
import random
import shapely.geometry as shg
from shapely.prepared import prep

random.seed(123)

province_df = gpd.read_file("./berkeley-s7dt3k-shapefile/s7dt3k.shp")
province_df.boundary.plot()

<AxesSubplot:>

Now I'll make some fake Google coordinates of the villages, which will be the true coordinates under the null hypothesis.

# Make Google coordinates, which are true under null hypothesis
minx,miny,maxx,maxy = province_df.geometry.unary_union.bounds

xs = [random.uniform(minx, maxx) for i in range(0,10000)]
ys = [random.uniform(miny, maxy) for i in range(0,10000)]

coords = pd.DataFrame({"google_x":xs,"google_y":ys})
coords['point'] = coords.apply(lambda x: shg.Point(x['google_x'],x['google_y']),axis=1)
coords = gpd.GeoDataFrame(coords,geometry='point')
coords = gpd.sjoin(province_df,coords,how="inner",op="contains")

<ipython-input-1-81c355b66054>:10: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.

Left CRS: EPSG:4326
Right CRS: None

  coords = gpd.sjoin(province_df,coords,how="inner",op="contains")

Next, we have to make the noisy coordinates, which will be displaced by a random direction and distance.

import math

# Choose diplacement distance
dists = [random.uniform(0,1) for i in range(0,len(coords.index))]
coords['dist'] = dists
# Choose displacement angles
angles = [random.uniform(0,2*math.pi) for i in range(0,len(coords.index))]
coords['angles'] = angles

# Make noisy coordinates
coords['noisy_x'] = coords['google_x']+coords['angles'].apply(math.cos)*coords['dist']
coords['noisy_y'] = coords['google_y']+coords['angles'].apply(math.sin)*coords['dist']

# Make sure this worked
coords['dist_check'] = ((coords['noisy_x']-coords['google_x'])**2+(coords['noisy_y']-coords['google_y'])**2)**(1/2)
coords[['dist','dist_check']]

	dist	dist_check
9	0.765766	0.765766
11	0.383786	0.383786
12	0.625888	0.625888
12	0.707815	0.707815
12	0.358853	0.358853
...	...	...
2861	0.068102	0.068102
2861	0.767153	0.767153
2861	0.470232	0.470232
2861	0.874593	0.874593
2861	0.239674	0.239674

4192 rows × 2 columns

Now we need to make the test statistic.

coords['noisy_within'] = coords.apply(lambda i: shg.Point(i['noisy_x'],i['noisy_y']).within(i['geometry']),axis=1)
T1 = coords['noisy_within'].mean()
coords['circle'] = coords.apply(lambda i: shg.Point(i['google_x'],i['google_y']).buffer(i['dist']),axis=1)

def test(circle,poly): # Thanks to our RA Yiyang Chen for this f'n
    circle_points = list(circle.exterior.coords)
    count=0
    for e in circle_points:
        if poly.contains(shg.Point(e)):
            count += 1
    return count/len(circle_points)

coords['circ_within'] = coords.apply(lambda i: test(i['circle'],i['geometry']),axis=1)
T2 = coords['circ_within'].mean()

T = T1-T2
T

0.0008367586611862099

As you can see, this is quite close to 0, suggesting that the Google coordinates are correct. When we applied this test to our true data we also got test statistic very close to zero which helped us decide to just use the Google coordinates.