Этот проект на Python автоматизирует создание гексогональных карт из любых геопространственных полигонов *.
Пример:
* Этот код не был тщательно протестирован. Пожалуйста, не стесняйтесь вносить свой вклад, предлагать улучшения или дайте мне знать, где в вашем применении этого кода есть ошибки, чтобы я мог сделать его лучше..
Проект на github
Алгоритм работы следующий:
- Считываем полигональный слой используя GeoPandas, создаем GeoDataFrame (т.е. GDF).
- Вычисляем центр для каждого полигона в GDF.
- Обобщаем геометрию полигонов, создав выпуклые оболочки.
- Создаем другую GDF с размытой геометрией краев. Это создаси единый полигональный объект, представляющий обобщенный контур входного объекта. Позже это будет использовано при создании шестиугольников.
- Добавьте геометрию центроида к обобщенному контуру.
- Вычислите минимальный ограничивающий прямоугольник, чтобы определить начальную точку создания шестиугольника.
- Создаем базовый шестиугольник в левом нижнем углу.
- Копируем этот шестиугольник, смещая его в направлениях X и Y так, чтобы он образовывал бесшовную шестиугольную сетку.
- Сохраняйте только те шестиугольники, центр тяжести которых находится в пределах обобщенного контура объектов.
- Если количество созданных шестиугольников не равно количеству полигонов во входных данных, измените размер шестиугольников и повторите попытку.
- Конечный продукт может быть выведен в любой геопространственный формат (Shapefile, GeoJSON, и т.д.).
Содержание
- Первоначальная настройка проекта
- Импортируем библиотеки
- Определите переменные для входных данных, имена полей и т.д.
- Первая группа — это имена полей (т.е. атрибуты), которые будут добавлены к входным данным (т.е. GeoPandas GeoDataFrames).
- Вторая группа специфична для входных данных. Это:
- Убедитесь, что путь к файлу является правильным
- Определите функцию для считывания геопространственного файла в GeoDataFrame
- Считываем геопространственный файл в GeoDataFrame, сортировка по admin_name/code, удаляем ненужные поля
- Определите полезные функции и запустите их.
- Создавайте геометрию для последующего пространственного сравнения
- Создайте обобщенную ограничивающую геометрию, чтобы позже убедиться, что шестиугольники соответствуют форме входных данных
- Отображение административных районов и центроидов
- Отображение выпуклых оболочек и ограничивающей геометрии
- Функции для получения некоторой базовой информации об объеме данных
- Функции, которые фактически создают шестиугольники
- Функция: displace_vertex
- Функция: build hexagon
- Функция: tesselate_hexagons_translate
- Функция: iterate_size
- Функция: verify_count
- Вызов функций для создания шестиугольников
- Установите входные параметры
- Отображение двух версий шестиугольных карт
- Снова отобразите две версии карт шестиугольников с перечисленными параметрами и заполненными шестиугольниками для проверки наличия пустот
- Отобразите первую шестиугольную карту с номерами плиток
- Функции для вычисления расстояния между центроидом (например, центроидом страны) и центроидами шестиугольников и административных районов
- Сопоставьте шестиугольные плитки и имена районов
- Отображение имен шестиугольников
- Другая забавная карта, которую можно распечатать и присвоить названия от руки
- Вывод шестиугольников в формат ГИС
Первоначальная настройка проекта
Импортируем библиотеки
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
import os import sys import math import operator import geopandas as gpd import pandas as pd import matplotlib.pyplot as plt import config as c from IPython.display import display, HTML from tabulate import tabulate from shapely.geometry import Polygon, Point from math import sqrt from shapely.affinity import translate |
Определите переменные для входных данных, имена полей и т.д.
Первая группа — это имена полей (т.е. атрибуты), которые будут добавлены к входным данным (т.е. GeoPandas GeoDataFrames).
- Их не нужно менять
Вторая группа специфична для входных данных. Это:
- Данные из файлов (shapefile, geojson, и т.д.).
- Имя атрибута, содержащее административные имена или коды.
- Имена атрибутов, от которых не следует избавляться во входных данных (для минимизации размера данных).
- Код EPSG, если он у вас есть. Если вы его не укажете, данные будут считываться с любой существующей пространственной привязки.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
# Имена атрибутов, добавляемые в фрейм геоданных во время обработки tile_id = 'tile_id' admin_id = 'admin_id' admin_id_admin = 'admin_id_admin' admin_id_tile = 'admin_id_tile' jo_tile = 'join_order_tile' jo_admin = 'join_order_admin' # Имена переменных для входных данных # Tatarstan polys = c.adm1 admin_name_col = 'HASC_1' keep_flds = ['NAME_1', admin_id] epsg = 4326 |
Убедитесь, что путь к файлу является правильным
|
1 |
os.path.isfile(polys) |
Определите функцию для считывания геопространственного файла в GeoDataFrame
- При этом будет предпринята попытка перепроектировать файл в указанную вами пространственную привязку EPSG
- Если повторное проектирование завершится неудачей, по умолчанию будет использоваться пространственная привязка, определенная во входных данных
|
1 2 3 4 5 6 7 8 9 |
def read_file_and_attempt_projection(file, epsg_code): try: gdf = gpd.read_file(file).to_crs({'init': 'epsg:{}'.format(epsg_code)}) print('Data projected to EPSG code: {}'.format(epsg_code)) except: gdf = gpd.read_file(file) print('*** Failed to project to EPSG code: {} \ \n*** Data read in with this CRS: {}'.format(epsg_code, gdf.crs)) return gdf |
Считываем геопространственный файл в GeoDataFrame, сортировка по admin_name/code, удаляем ненужные поля
|
1 2 3 4 5 6 7 8 9 10 11 12 |
admin = read_file_and_attempt_projection(polys, epsg) admin[admin_id] = admin['district'].str[0:6] #admin[admin_id] = admin[admin_name_col].str.split('.').str[1] #For USA #admin[admin_id] = admin[admin_name_col].str.split(' ').str[0] #For India #admin[admin_id] = admin[admin_name_col] #For Kenya #admin = admin.sort_values(by=admin_name_col).reset_index(drop=True) admin = admin.rename(columns={'geometry': 'geom_poly'}).set_geometry('geom_poly') #keep_flds.append(admin.geometry.name) #del_flds = [x for x in list(admin) if x not in keep_flds] #admin.drop(columns=del_flds, inplace=True) print(display(HTML(admin.head(5).to_html()))) |
Определите полезные функции и запустите их.
Создавайте геометрию для последующего пространственного сравнения
- Добавьте геометрию центроида к каждому многоугольнику
Создайте обобщенную ограничивающую геометрию, чтобы позже убедиться, что шестиугольники соответствуют форме входных данных
- Создайте геометрию выпуклой оболочки (т.е. минимальную ограничивающую геометрию)
- Сведите геометрию выпуклой оболочки к одному многоугольнику
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
def display_point_labels(gdf, label_col, fs=10, c='blue', angle=0): texts = [] for x, y, label in zip(gdf.geometry.x, gdf.geometry.y, gdf[label_col]): texts.append(plt.text(x, y, label, fontsize=fs, color=c, rotation=angle, horizontalalignment='center', verticalalignment='center')) def add_centroid_geom(gdf): gdf['geom_centroid'] = gdf.centroid return gdf def add_convexhull_geom(gdf): gdf['geom_convexhull'] = gdf.convex_hull return gdf def dissolve_polygons(gdf): gdf['dissolve_col'] = 1 gdf_dissolve = gdf.dissolve(by='dissolve_col') #gdf_dissolve['geom_exterior'] = Polygon(gdf_dissolve.exterior[1]) #gdf_dissolve = gdf_dissolve.set_geometry('geom_exterior') gdf_dissolve['geom_centroid'] = gdf_dissolve.centroid return gdf_dissolve.reset_index(drop=True) admin = add_centroid_geom(admin) admin = add_convexhull_geom(admin) admin_hulls_dissolved = dissolve_polygons(admin.set_geometry('geom_convexhull')) print(display(HTML(admin_hulls_dissolved.to_html()))) |
Отображение административных районов и центроидов
|
1 2 3 4 5 6 7 |
fig, ax = plt.subplots(figsize=(10, 10)) # Plot on axis 1 (left) ax.set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax, linewidth=0.5, edgecolor='lightgrey') admin.set_geometry('geom_centroid').plot(ax=ax, color='red') plt.show() |
Отображение выпуклых оболочек и ограничивающей геометрии
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
fig, (ax1,ax2) = plt.subplots(nrows=1, ncols=2, figsize=(20, 20), sharex=True, sharey=True) # Plot on axis 1 (left) ax1.set_axis_off() admin.set_geometry('geom_convexhull').plot(ax=ax1, linewidth=0.5, edgecolor='black', color='lightgrey') admin_hulls_dissolved.set_geometry('geom_centroid').plot(ax=ax1, color='red') # Plot on axis 2 (right) ax2.set_axis_off() admin_hulls_dissolved.plot(ax=ax2, linewidth=0.5, edgecolor='black', color='lightgrey') admin_hulls_dissolved.set_geometry('geom_centroid').plot(ax=ax2, color='red') #admin.loc[admin[admin_id] == 'ME'].set_geometry('geom_centroid').plot(ax=ax2, color='black') # Only for USA plt.show() |
Функции для получения некоторой базовой информации об объеме данных
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 |
def round_dict_vals(dictionary, places=10): for k, v in dictionary.items(): dictionary[k] = float(round(dictionary[k], places)) return dictionary def get_bounding_coordinates(gpd, round_nums=True): tb = gpd.total_bounds total_bounds = {'minx': tb[0], 'miny': tb[1], 'maxx': tb[2], 'maxy': tb[3], 'rangex': abs(tb[0] - tb[2]), 'rangey': abs(tb[1] - tb[3])} if round_nums: total_bounds = round_dict_vals(total_bounds) return total_bounds def build_bounding_polygon(extent, crs): BL = (extent['minx'], extent['miny']) BR = (extent['maxx'], extent['miny']) TR = (extent['maxx'], extent['maxy']) TL = (extent['minx'], extent['maxy']) corners = [BL, BR, TR, TL] rectangle = gpd.GeoDataFrame(pd.DataFrame({'geom_poly': [Polygon(corners)]}), geometry='geom_poly', crs=crs) return rectangle bbox = get_bounding_coordinates(admin.set_geometry('geom_poly')) bbox_gdf = build_bounding_polygon(bbox, admin.crs) bbox |
Функции, которые фактически создают шестиугольники
Функция: displace_vertex
- Эта функция перемещает точку на заданное расстояние и угол.
- 6 точек создают вершины для шестиугольника.
Функция: build hexagon
- Эта функция вызывает displace_vertex 5 раз, чтобы создать 5 точек в дополнение к начальной точке
- Затем соединяет эти вершины, чтобы создать многоугольник
Функция: tesselate_hexagons_translate
- Эта функция использует отдельный шестиугольник, созданный ранее, копирует его и переводит (т.е. сдвигает) в направлении x, создавая строку.
- Когда строка завершена, она начинается с исходного шестиугольника и перемещает его в направлениях x и y, чтобы создать первый шестиугольник в следующем ряду.
- Это повторяется до тех пор, пока оба измерения x и y не будут заполнены шестиугольниками.
Функция: iterate_size
- Эта функция итеративно изменяет (уменьшает) размер шестиугольников до тех пор, пока количество созданных шестиугольников не будет соответствовать или станет меньше количества административных областей во входных данных
Функция: verify_count
- Это простая функция для вывода инструкции, указывающей, является ли количество созданных шестиугольников правильным или их слишком мало.
- Если их слишком мало, пользователю предлагается изменить некоторые параметры
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 |
def displace_vertex(x, y, length, angle): ''' calculates new point at a given distance away from original point. All values should be provided in meters point = (x,y) ''' #The number pi PI = 3.1415926535 #Convert the random angle from degrees to radians angle_radian = (angle) * (PI/180) #Generate the offset by applying trig formulas (law of cosines) #using the distance as the hypotenuse solving for the other sides xOffset = math.sin(angle_radian) * length yOffset = math.cos(angle_radian) * length #Add the offset to the orginal coordinate new_x = x + xOffset new_y = y + yOffset return (new_x, new_y) def build_hexagon(x, y, length, bearing): pts = [(x,y)] for pt in range(0,5): x, y = displace_vertex(x, y, length, bearing) bearing += 60 pts.append((x,y)) hex_poly = Polygon(pts) return hex_poly def tesselate_hexagons_translate(bb, hull, length, x_offset_factor, y_offset_factor, flat, keep): x_row = 100 y_col = 0 hexes = [] x_offset_seed = x_offset_factor * length y_offset_seed = y_offset_factor * length # flat='up' means flat sides are up and down # flat='side' means flat sides are on the sides if flat=='up': bearing = 30 next_start_pt_idx = 3 elif flat=='side': bearing = 0 next_start_pt_idx = 4 else: raise Exception('Acceptable parameters for flat: \'up\' and \'side\'') hexagon_seed = build_hexagon(bb['minx'] + x_offset_seed, bb['miny'] + y_offset_seed, length, bearing) hexagon = hexagon_seed row=1 shift = 1 xy_inside_bbox = True while xy_inside_bbox: y_col += 1 hex_gdf = gpd.GeoDataFrame(pd.DataFrame({'geom_poly': [hexagon]}), geometry='geom_poly') hex_gdf[tile_id] = (x_row + y_col) centroid = hexagon.centroid x_inside = ((centroid.x >= bb['minx'] - abs(x_offset_seed)) and (centroid.x <= bb['maxx'] + abs(x_offset_seed))) y_inside = ((centroid.y >= bb['miny'] - abs(y_offset_seed)) and (centroid.y <= bb['maxy'] + abs(y_offset_seed))) if keep=='center_in' and hull.intersects(centroid): hexes += [hex_gdf] elif keep=='touches' and hull.intersects(hexagon): hexes += [hex_gdf] if x_inside and y_inside: if flat == 'up': move_x = 3 * length elif flat == 'side': move_x = sqrt(3) * length move_y = 0 elif not x_inside and y_inside: hexagon = hexagon_seed if flat == 'up': if row%2 == 0: #Every other row move_y = (sqrt(3) * length * shift) + y_offset_seed move_x = 0 + x_offset_seed shift+=1 else: move_y = (((sqrt(3) * length) * shift) - (0.5 * sqrt(3) * length)) + y_offset_seed move_x = (length * 1.5) + x_offset_seed elif flat == 'side': if row%2 == 0: #Every other row move_y = 1.5 * length * row move_x = 0# + x_offset_seed else: move_y = 1.5 * length * row move_x = -((0.5 * sqrt(3) * length))# + x_offset_seed) x_row += 100 y_col = 0 row+=1 else: xy_inside_bbox = False hexagon = translate(hexagon, xoff=move_x, yoff=move_y) hexes_gdf = gpd.GeoDataFrame(pd.concat(hexes, ignore_index=True, sort=False), crs=hexes[0].crs, geometry='geom_poly') hexes_gdf = add_centroid_geom(hexes_gdf) return hexes_gdf def iterate_size(adm, bb, flat, x_offset_factor, y_offset_factor, len_factor, iter_factor, keep='touches'): length = bb['rangex'] / len_factor print('\nFlat side: {}'.format(flat)) hex_count = len(adm) + 1 i = 1 while hex_count > len(adm): hexagons = tesselate_hexagons_translate(bb=bbox, hull=admin_hulls_dissolved['geom_convexhull'][0], length=length, x_offset_factor=x_offset_factor, y_offset_factor=y_offset_factor, flat=flat, keep=keep) hex_count = len(hexagons) print('i{} h{} |'.format(i, hex_count), end =' ') length = length * iter_factor i+=1 hexagons.crs = adm.crs return hexagons def verify_count(hexes, num): # hexes is passed in as a list of geodataframes print('\n') i=1 for h in hexes: if len(h) == num: print('Cool, we got the correct number of hexagons, {}, for item #{}'.format(num, i)) else: print('*** Dang, we got {} hexagons but needed {} for item #{}. Try increasing lfac variable, or decreasing ifac variable.'.format(len(h), num, i)) i+=1 |
Вызов функций для создания шестиугольников
Установите входные параметры
- flat_side*: Определяет ориентацию шестиугольников.
- ‘up’ означает, что шестиугольники будут иметь плоские стороны вверх и вниз
- ‘side’ означает, что шестиугольники будут иметь плоские стороны справа и слева
- xofac*: коэффициент смещения x — определяет смещение в направлении x для начального шестиугольника
- yofac*: коэффициент смещения по оси y — определяет смещение в направлении y для начального шестиугольника
- lfac: коэффициент длины — определяет размер первого шестиугольника путем деления ширины экстента на lfac
- ifac: коэффициент масштабирования — значение, умноженное на длину шестиугольника для итеративного увеличения / уменьшения его размера до тех пор, пока не будет создано нужное количество шестиугольников.
*В этом примере создаются две версии шестиугольных карт. Вы можете создать больше, создав дополнительные *flat_side, xofac, и yofac переменные, затем дополнительно вызовите iterate_size **
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
# Hexagon creation variables flat_side1 = 'up' # Hexagon orientation, valid values: ;'up', and 'side' x_off1 = -0.1 # X offset factor 1 - moves seed hexagon in the X direction. Ideally between -1 and 1 y_off1 = 0.1 # Y offset factor 1 - moves seed hexagon in the Y direction. Ideally between -1 and 1 lfac1 = 27.6 # Length factor - width / lfac - determines the initial number of hexagons # - larger value creates more (i.e. smaller) initial hexagons ifac1 = 1.0005 # Iteration factor - multiplied to hexagon length to create larger hexagons each iteration # - Value > 1 keep_hex1 = 'center_in' # Variable to determine which hexagons to keep. Valid values are: 'center_in', and 'touches' # If creating a second hexagon configuration set variables here. (refer to the above for explanations) flat_side2 = 'side' x_off2 = -0.5 y_off2 = 0.9 lfac2 = 27.6 ifac2 = 1.0005 keep_hex2 = 'center_in' print('We need {} hexagons.'.format(len(admin))) hexagons1 = iterate_size(admin, bbox, flat=flat_side1, x_offset_factor=x_off1, y_offset_factor=y_off1, len_factor=lfac1, iter_factor=ifac1, keep=keep_hex1) hexagons2 = iterate_size(admin, bbox, flat=flat_side2, x_offset_factor=x_off2, y_offset_factor=y_off2, len_factor=lfac2, iter_factor=ifac2, keep=keep_hex2) verify_count([hexagons1, hexagons2], len(admin)) #hexagons1.head(3) |
Отображение двух версий шестиугольных карт
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
fig, (ax1,ax2) = plt.subplots(nrows=1, ncols=2, figsize=(25, 20), sharex=True, sharey=True) # Plot on axis 1 (left) ax1.set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax1, linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') hexagons1.set_geometry('geom_centroid').plot(ax=ax1) hexagons1.plot(ax=ax1, edgecolor='black', facecolor="none") bbox_gdf.plot(ax=ax1, facecolor="none", edgecolor='lightgrey', linestyle=':') admin_hulls_dissolved.plot(ax=ax1, linewidth=1, edgecolor='red', facecolor='none', linestyle=':') #pt.plot(ax=ax1, color='red') # Plot on axis 2 (left) ax2.set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax2, linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') hexagons2.set_geometry('geom_centroid').plot(ax=ax2) hexagons2.plot(ax=ax2, edgecolor='black', facecolor="none") bbox_gdf.plot(ax=ax2, facecolor="none", edgecolor='lightgrey', linestyle=':') admin_hulls_dissolved.plot(ax=ax2, linewidth=1, edgecolor='red', facecolor='none', linestyle=':') plt.show() |
Снова отобразите две версии карт шестиугольников с перечисленными параметрами и заполненными шестиугольниками для проверки наличия пустот
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 |
fig, ax = plt.subplots(nrows=2, ncols=2, constrained_layout=True, figsize=(15, 10), sharex=True, sharey=True) title = '''tile count: {}\n x offset factor: {}\n y offset factor: {}\n length factor: {}\n iteration factor: {}''' # Plot on axis 1 (top left) ax[0,0].set_title(title.format(len(hexagons1), x_off1, y_off1, lfac1, ifac1), fontdict={'fontsize': '15', 'fontweight' : '500'}) ax[0,0].set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax[0,0], linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') bbox_gdf.plot(ax=ax[0,0], facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons1.plot(ax=ax[0,0], edgecolor='black', facecolor="none") admin_hulls_dissolved.set_geometry('geom_centroid').plot(ax=ax[0,0], color='red') # Plot on axis 2 (top right) ax[0,1].set_title(title.format(len(hexagons2), x_off2, y_off2, lfac2, ifac2), fontdict={'fontsize': '15', 'fontweight' : '500'}) ax[0,1].set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax[0,1], linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') bbox_gdf.plot(ax=ax[0,1], facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons2.plot(ax=ax[0,1], edgecolor='black', facecolor="none") admin_hulls_dissolved.set_geometry('geom_centroid').plot(ax=ax[0,1], color='red') # Plot on axis 3 (bottom left) ax[1,0].set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax[1,0], linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') bbox_gdf.plot(ax=ax[1,0], facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons1.plot(ax=ax[1,0], edgecolor='black', facecolor="skyblue") admin_hulls_dissolved.set_geometry('geom_centroid').plot(ax=ax[1,0], color='red') # Plot on axis 4 (bottom right) ax[1,1].set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax[1,1], linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') bbox_gdf.plot(ax=ax[1,1], facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons2.plot(ax=ax[1,1], edgecolor='black', facecolor="skyblue") admin_hulls_dissolved.set_geometry('geom_centroid').plot(ax=ax[1,1], color='red') plt.show() |
Отобразите первую шестиугольную карту с номерами плиток
- Последние два числа = номер столбца
- Первые 1-2 числа = номер строки
|
1 2 3 4 5 6 7 8 9 10 11 12 |
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 10)) ax.set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax, linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') bbox_gdf.plot(ax=ax, facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons1.plot(ax=ax, edgecolor='black', facecolor="none") admin_hulls_dissolved.set_geometry('geom_centroid').plot(ax=ax, color='red') #hexagons1.set_geometry('geom_centroid').plot(ax=ax, color='lightblue') #display_point_labels(admin.set_geometry('geom_centroid'), admin_id, c='black') display_point_labels(hexagons1.set_geometry('geom_centroid'), tile_id) #admin.set_geometry('geom_centroid').plot(ax=ax, color='green') |
Функции для вычисления расстояния между центроидом (например, центроидом страны) и центроидами шестиугольников и административных районов
- Это позволяет нам ранжировать и упорядочивать шестиугольные плитки и административные полигоны, чтобы присваивать имена районов шестиугольникам
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 |
def dist_to_centroid_fld(centers, centroid, idx_col, order_col, sort_order): distances = {} print(type(centroid['geom_centroid'])) centroid_pt = centroid['geom_centroid'][0] #Get Shapely point geom from geodataframe try: centers.drop(columns=['dist_to_centroid', order_col], inplace=True) except: pass for i, row in centers.iterrows(): pt = row['geom_centroid'] idx = row[idx_col] distance = pt.distance(centroid_pt) distances[idx] = distance #print('{} {} {}'.format(i, idx, distance)) if sort_order == 'farthest_first': rev = True elif sort_order == 'closest_first': rev = False distances_dict = dict(sorted(distances.items(), key=operator.itemgetter(1), reverse=rev)) distances_df = pd.DataFrame(list(distances_dict.items()), columns=[idx_col, 'dist_to_centroid']) distances_df[order_col] = distances_df.index centers = centers.join(distances_df.set_index(idx_col), on=idx_col).sort_values(by=[order_col]).reset_index(drop=True) return centers def build_admin_tile_distance_df(admin_centers, tile_centers): cols = [admin_id, tile_id, 'distance'] df_dist = pd.DataFrame(columns=cols, dtype=int) for idx1, admin_row in admin_centers.iterrows(): admin_pt = admin_row['geom_centroid'] admin_idx = admin_row[admin_id] for idx2, tile_row in tile_centers.iterrows(): tile_pt = tile_row['geom_centroid'] tile_idx = tile_row[tile_id] distance = tile_pt.distance(admin_pt) #print(distance) df_tmp = pd.DataFrame([[admin_idx, tile_idx, distance]], columns=cols) df_dist = pd.concat([df_dist, pd.DataFrame(df_tmp, index=[0])]) #df_dist.append(df_tmp, ignore_index=True) return df_dist admin = dist_to_centroid_fld(admin, admin_hulls_dissolved, admin_id, jo_admin, 'farthest_first') hexagons1 = dist_to_centroid_fld(hexagons1, admin_hulls_dissolved, tile_id, jo_tile, 'farthest_first') distance_df = build_admin_tile_distance_df(admin, hexagons1) display(HTML(admin[[admin_id, jo_admin]].head(10).to_html())) display(HTML(hexagons1[[tile_id, jo_tile]].head(10).to_html())) #display(HTML(distance_df.head(50).to_html())) # display(HTML(admin.head(3).to_html())) # display(HTML(hexagons1.head(3).to_html())) |
Сопоставьте шестиугольные плитки и имена районов
Используются два метода
- admin-distance-to-centroid
- tile-distance-to-centroid prioritization
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
def match_one_to_another(admin_centers, tile_centers, dist_df, order_defined_by): if order_defined_by == 'admin': cols = [admin_id_admin, tile_id] df_match = pd.DataFrame(columns=cols, dtype=int) for i, row in admin_centers[[admin_id, jo_admin]].iterrows(): admin_idx = row[admin_id] df = dist_df.loc[dist_df[admin_id] == admin_idx] tile_idx = df.loc[df.groupby(admin_id)['distance'].idxmin()][tile_id].values[0] dist_df = dist_df[dist_df[tile_id] != tile_idx] df_tmp = pd.DataFrame([[admin_idx, tile_idx]], columns=cols) df_match = pd.concat([df_match, pd.DataFrame(df_tmp, index=[0])]) #df_match = df_match.append(df_tmp, ignore_index=True) if order_defined_by == 'tile': cols = [admin_id_tile, tile_id] df_match = pd.DataFrame(columns=cols, dtype=int) for i, row in tile_centers[[tile_id, jo_tile]].iterrows(): tile_idx = row[tile_id] df = dist_df.loc[dist_df[tile_id] == tile_idx] admin_idx = df.loc[df.groupby(tile_id)['distance'].idxmin()][admin_id].values[0] dist_df = dist_df[dist_df[admin_id] != admin_idx] df_tmp = pd.DataFrame([[admin_idx, tile_idx]], columns=cols) df_match = pd.concat([df_match, pd.DataFrame(df_tmp, index=[0])]) #df_match = df_match.append(df_tmp, ignore_index=True) return df_match farthest_states_match = match_one_to_another(admin, hexagons1, distance_df, 'admin') farthest_tile_match = match_one_to_another(admin, hexagons1, distance_df, 'tile') display(HTML(farthest_states_match.head(5).to_html())) display(HTML(farthest_tile_match.head(5).to_html())) hexagons1 = hexagons1.merge(farthest_states_match, on=tile_id) hexagons1 = hexagons1.merge(farthest_tile_match, on=tile_id) display(HTML(hexagons1.head(5).to_html())) |
Отображение имен шестиугольников
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
fig, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(11, 10)) ax1.set_title('Farthest state from the US center is assigned to the nearest hexagon.\n(State location prioritization)', fontdict={'fontsize': '15', 'fontweight' : '500'}) ax1.set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax1, linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') #bbox_gdf.plot(ax=ax1, facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons1.plot(ax=ax1, edgecolor='black', facecolor="none") hexagons1[tile_id] = hexagons1.index display_point_labels(hexagons1.set_geometry('geom_centroid'), admin_id_admin, c='black', angle=20) fig, ax2 = plt.subplots(nrows=1, ncols=1, figsize=(12, 10)) ax2.set_title('Farthest hexagon from the US center is assigned value of its nearest state.\n(Hexagon location prioritization)', fontdict={'fontsize': '15', 'fontweight' : '500'}) ax2.set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax2, linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') #bbox_gdf.plot(ax=ax2, facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons1.plot(ax=ax2, edgecolor='black', facecolor="none") hexagons1[tile_id] = hexagons1.index display_point_labels(hexagons1.set_geometry('geom_centroid'), admin_id_tile, c='black', angle=20) plt.show() #fig2.savefig('USA_hexes_labeled.png', bbox_inches='tight', dpi=300) # For saving the figure as a PNG |
Другая забавная карта, которую можно распечатать и присвоить названия от руки
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 10)) #ax.set_title('States', # fontdict={'fontsize': '15', 'fontweight' : '500'}) ax.set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax, linewidth=0.5, edgecolor='black', color='whitesmoke') #bbox_gdf.plot(ax=ax, facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons1.plot(ax=ax, edgecolor='lightgrey', facecolor="none") hexagons1[tile_id] = hexagons1.index display_point_labels(admin.set_geometry('geom_centroid'), admin_id, c='black') fig, ax2 = plt.subplots(nrows=1, ncols=1, figsize=(12, 10)) #ax2.set_title('Blank', # fontdict={'fontsize': '15', 'fontweight' : '500'}) ax2.set_axis_off() admin.set_geometry('geom_poly').plot(ax=ax2, linewidth=0.5, edgecolor='lightgrey', color='whitesmoke') #bbox_gdf.plot(ax=ax2, facecolor="none", edgecolor='lightgrey', linestyle=':') hexagons1.plot(ax=ax2, edgecolor='black', facecolor="none") hexagons1[tile_id] = hexagons1.index |
Вывод шестиугольников в формат ГИС
|
1 |
hexagons1.drop(columns=['geom_centroid']).to_file('output/USA_hexagons9.geojson', driver='GeoJSON') |
