Рубрики
Без рубрики

Обратное Геокодирование

Недавно я работал над анализом на основе местоположения с клиентом, где у нас было несколько записей с недействительными почтовыми кодами. Однако у нас были широта и долгота для каждой из этих записей. Это попало…

Автор оригинала: Bob Haffner.

Недавно я работал над анализом на основе местоположения с клиентом, где у нас было несколько записей с недействительными почтовыми кодами. Однако у нас были широта и Долгота для каждой из этих записей. Это заставило меня задуматься о процессе обратного геокодирования, чтобы мы могли обновить эти плохие молнии.

На ум пришли два подхода. Один из них использует KDTree для поиска известного набора координат и почтовых индексов для ближайшего совпадения. Второй подход будет ГеоПандой и стройной комбинацией для выполнения некоторых пространственных операций.

Оба подхода будут использовать последние данные о районах табуляции почтовых индексов переписи населения США

Давайте начнем с одной точки для поиска

import numpy as np
omaha_point = np.array((-95.995102, 41.257160))

Использование подхода к дереву KD с фреймом данных Pandas

import pandas as pd
from sklearn.neighbors import KDTree

url = 'http://www2.census.gov/geo/docs/maps-data/data/gazetteer/2016_Gazetteer/2016_Gaz_zcta_national.zip'df_locations = pd.read_csv(url, dtype={'GEOID' : 'str'},sep='\t')
#some column cleanupdf_locations.columns = df_locations.columns.str.strip()
print (len(df_locations))print(df_locations.head())

33144
GEOID INTPTLAT INTPTLONG
00601 18.180555 -66.749961
00602 18.361945 -67.175597 
00603 18.455183 -67.119887 
00606 18.158345 -66.932911 
00610 18.295366 -67.12513

Я делаю предположение, что координаты представляют центроиды их соответствующих молний в этом наборе данных.

Теперь мы построим наше простое KDTree, используя реализацию scikit-learn. Принимая значения по умолчанию здесь. Смотрите здесь для получения дополнительной информации и примеров.

kdt = KDTree(df_locations[['INTPTLONG', 'INTPTLAT']])

И все, наше дерево построено и готово к использованию. Вы должны запросить дерево с помощью двумерного массива, поэтому сначала мы обратимся к нему с помощью expand_dims numpy. Затем мы выполняем запрос, потому что мне просто нужен ближайший центроид zip к моей точке поиска в Омахе. Я использую результат запроса для индексирования в наш фрейм данных и получения действительного почтового индекса 68132.

omaha_point_kdt = np.expand_dims(omaha_point, axis=0)

nearest_point_index = kdt.query(omaha_point_kdt, k=1, return_distance=False)
print(df_locations.loc[nearest_point_index[0], 'GEOID'])

23609 68132
Name: GEOID, dtype: object

К геопандам и стройному подходу. Здесь мы используем довольно большой файл shp из Переписи в качестве справочного набора данных. Этот файл содержит те же молнии, что и в нашем первом подходе, но в нем есть полигоны, которые представляют границы молний, а не просто центроиды.

import geopandas as gpd
from shapely.geometry import Point

gdf_locations = gpd.read_file('data/tl_2016_us_zcta510/tl_2016_us_zcta510.shp')print(len(gdf_locations))
print(gdf_locations[['GEOID10', 'geometry']].head())

33144

GEOID10 geometry
43451 POLYGON ((-83.674464 41.331119, -83.6744449999...
43452 POLYGON ((-83.067745 41.537718, -83.067729 41....
43456 (POLYGON ((-82.8566 41.681222, -82.856831 41.6...
43457 POLYGON ((-83.467474 41.268186, -83.4676039999...
43458 POLYGON ((-83.222292 41.531025, -83.2222819999...

Этот файл shp занял 2,5 минуты, чтобы прочитать его на моей машине. Довольно многие молнии являются мультиполигонами, которые, я думаю, способствовали длительному времени преобразования SHP-файла в геодатафрам.

Затем мы изменим геокодирование, но сначала нам нужно преобразовать нашу точку omaha_point в фигурную точку. Затем мы используем одну из моих любимых вещей в GeoPandas, которая заключается в возможности комбинировать логическую индексацию с пространственными операциями Shapely. Таким образом, наш почтовый полигон 68132 содержит нашу точку omaha_point. Довольно круто, я бы сказал.

omaha_point_shp = Point(omaha_point)

filter = gdf_locations['geometry'].contains(omaha_point_shp)
print(gdf_locations.loc[filter, 'GEOID10'])

24842 68132
Name: GEOID10, dtype: object

Два разных подхода дают один и тот же результат (68132). Давайте посмотрим на некоторые тайминги.

Подход к дереву KD

%%timeit
nearest_point_index = kdt.query(omaha_point_kdt, k=1, return_distance=False)
df_locations.loc[nearest_point_index[0], 'GEOID']

494 µs ± 24.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Геопанды/Стройный подход

%%timeit
filter = gdf_locations['geometry'].contains(omaha_point_shp)
gdf_locations.loc[filter, 'GEOID10']

243 ms ± 15.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Решение дерева KD в этом случае значительно быстрее. Хотя я предпочитаю геопанды и стройное решение, потому что оно более точное. Я мог видеть случаи, когда Lat Long может быть ближе к центроиду zips, но находится в пределах границы другого zip. В конце концов, я думаю, что у обоих есть некоторые достоинства.