Недавно я работал над анализом на основе местоположения с клиентом, где у нас было несколько записей с недействительными почтовыми кодами. Однако у нас были широта и Долгота для каждой из этих записей. Это заставило меня задуматься о процессе обратного геокодирования, чтобы мы могли обновить эти плохие молнии.
На ум пришли два подхода. Один из них использует 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. В конце концов, я думаю, что у обоих есть некоторые достоинства.