Как выполнять ГИС-вычисления в базах данных, отличных от ГИС

Введение:

Иногда вы можете столкнуться с необходимостью выполнения геопространственных функций в вашем приложении, таких как картографирование местоположений пользователей или анализ географических данных. Для этих задач доступно множество специфичных для языка библиотек, таких как GDAL, Shapely и Geopandas для Python.

В качестве альтернативы геопространственные функции могут быть реализованы через базы данных; например, вы можете использовать расширение PostGIS с реляционной базой данных, такой как PostgreSQL, или использовать встроенную поддержку пространственных типов данных в распределенной базе данных, такой как Azure CosmosDB.

Однако если ваше серверное хранилище, такое как Redis или Google Spanner, не поддерживает пространственные запросы изначально и вам необходимо обрабатывать крупномасштабные геопространственные запросы, эта статья предназначена для вас.

Каковы мои варианты?

Вы всегда можете создать еще один микросервис для обработки пространственных данных, но этот вариант часто влечет за собой затраты на поддержку дополнительного приложения. Другой подход — использовать библиотеки геоиндексации, такие как S2 и H3. S2, разработанный Google, основан на кривой Гильберта, а H3, разработанный Uber, основан на геодезической дискретной глобальной сетке. S2 и H3 имеют много общего: оба делят данную область на ячейки и используют 64-битные целые числа для индексации этих ячеек.

Однако основное отличие заключается в форме ячеек; В S2 используются ячейки квадратной формы, тогда как в H3 используются ячейки шестиугольной формы. Для некоторых приложений H3 может обеспечить более высокую производительность. Однако в целом любой библиотеки должно быть достаточно. В этой статье мы будем использовать S2, но аналогичные функции можно выполнять и с помощью H3.

Основные понятия библиотеки Google S2

• Ячейки: S2 делит сферу на ячейки, каждая из которых имеет уникальный 64-битный идентификатор.

• Уровни ячеек. Иерархия допускает различные уровни детализации: от больших регионов до небольших точных областей. Каждый уровень представляет собой разное разрешение:

  
  

  • Уровень 0: Самые крупные клетки, покрывающие значительную часть поверхности Земли.

    
    

  • Более высокие уровни: клетки постепенно подразделяются на более мелкие квадранты. Например, каждая ячейка уровня 1 делится на четыре ячейки уровня 2 и так далее.

    
    

  • Разрешение и площадь: более высокие уровни соответствуют более высокому разрешению и меньшим площадям ячеек. Эта иерархия обеспечивает точную индексацию и запросы на различных уровнях детализации.

    
    

В таблице ниже показаны различные уровни ячеек и соответствующие им области.

Из представленной таблицы видно, что с помощью S2 можно добиться точности отображения до 0,44 см^2. Внутри каждого квадрата ячейки S2 существует дочерняя ячейка, имеющая одного и того же родителя, что указывает на иерархическую структуру. Уровень ячейки может быть статическим значением (один и тот же уровень применяется ко всем ячейкам) или может быть динамическим, когда S2 решает, какое разрешение работает лучше всего.

Вычисление ближайших соседей

Начнем с примера. Предположим, мы пишем приложение, которое предоставляет функции, подобные сервису близости, для района Сиэтла. Мы хотим вернуть список кофеен в заданном районе. Для выполнения этих операций разделим эту задачу на 4 подзадачи:

• Загрузка карты Сиэтла
• Визуализация ячеек S2 на карте Сиэтла
• Сохраните несколько мест кофеен в базе данных.
• Запрос ближайших кофеен

Загрузка карты Сиэтла

Чтобы загрузить карту Google, мы будем использовать библиотеку gmplot. Для загрузки этой библиотеки требуется ключ API Google Maps. Чтобы сгенерировать ключ API, следуйте инструкциям здесь .

 import gmplot import const # plot seattle with zoom level 13 gmap = gmplot.GoogleMapPlotter(47.6061, -122.3328, 13, apikey=const.API_KEY) # Draw the map to an HTML file: gmap.draw('map.html')

Приведенный выше код создает файл map.html, как показано ниже:

Визуализируйте ячейки S2 на карте Сиэтла

Теперь, когда у нас есть карта, давайте нарисуем несколько ячеек S2 для карт:

 from s2 import * import gmplot # plot seattle with zoom level 13 gmap = gmplot.GoogleMapPlotter(47.6061, -122.3328, 13, apikey=const.API_KEY) areatobeplotted = [ (47.64395531736767,-122.43597221319135), (47.51369277846956,-122.43597221319135), (47.51369277846956,-122.24156866779164), (47.64395531736767,-122.24156866779164), (47.64395531736767,-122.43597221319135) ] region_rect = S2LatLngRect( S2LatLng.FromDegrees(47.51369277846956,-122.43597221319135), S2LatLng.FromDegrees(47.64395531736767, -122.24156866779164)) coverer = S2RegionCoverer() coverer.set_min_level(8) coverer.set_max_level(15) covering = coverer.GetCovering(region_rect) geoms = 0 for cellid in covering: new_cell = S2Cell(cellid) vertices = [] for i in range(0, 4): vertex = new_cell.GetVertex(i) latlng = S2LatLng(vertex) vertices.append((latlng.lat().degrees(), latlng.lng().degrees())) gmap.polygon(*zip(*vertices), face_color='pink', edge_color='cornflowerblue', edge_width=5) geoms+=1 gmap.polygon(*zip(*areatobeplotted), face_color='red', edge_color='green', edge_width=5) print(f"Total Geometries: {geoms}") gmap.draw('/tmp/map.html')

 Output: Total Geometries: 273

В приведенном выше коде мы сначала центрируем плоттер Google Map вокруг района Сиэтла. В S2RegionCoverer мы инициализируем средство покрытия региона, чтобы оно имело динамические уровни от минимального уровня 8 до максимального уровня 15. Это позволяет S2 динамически подгонять все ячейки под определенные размеры ячеек для наилучшего соответствия. Метод GetCovering возвращает покрытие прямоугольника вокруг области Сиэтла.

Затем мы перебираем каждую ячейку, вычисляя вершины ячеек и нанося их на карту. Мы сохраняем количество сгенерированных ячеек примерно 273. Наконец, мы рисуем входной прямоугольник красным. Этот код будет отображать ячейки S2 на карте Сиэтла в /tmp/map.html , как показано ниже:

Сохраните несколько мест кофеен в базе данных

Давайте сгенерируем базу данных кофеен вместе с их идентификаторами ячеек S2. Вы можете хранить эти ячейки в любой базе данных по вашему выбору. В этом уроке мы будем использовать базу данных данных SQLite. В приведенном ниже примере кода мы подключаемся к базе данных SQLite, чтобы создать таблицу CoffeeShops с тремя полями Id , name и cell_id .

Как и в предыдущем примере, мы используем S2RegionCoverer для расчета ячеек, но на этот раз мы используем фиксированный уровень для построения точек. Наконец, вычисленный идентификатор преобразуется в строку и сохраняется в базе данных.

 import sqlite3 from s2 import S2CellId,S2LatLng,S2RegionCoverer # Connect to SQLite database conn = sqlite3.connect('/tmp/sqlite_cells.db') cursor = conn.cursor() # Create a table to store cell IDs cursor.execute('''CREATE TABLE IF NOT EXISTS CoffeeShops ( id INTEGER PRIMARY KEY AUTOINCREMENT, name TEXT, cell_id TEXT )''') coverer = S2RegionCoverer() # Function to generate S2 cell ID for a given latitude and longitude def generate_cell_id(latitude, longitude, level=16): cell=S2CellId(S2LatLng.FromDegrees(latitude, longitude)) return str(cell.parent(level)) # Function to insert cell IDs into the database def insert_cell_ids(name,lat,lng): cell_id = generate_cell_id(lat, lng) cursor.execute("INSERT INTO CoffeeShops (name, cell_id) VALUES (?, ?)", (name, cell_id)) conn.commit() # Insert cell IDs into the database insert_cell_ids("Overcast Coffee", 47.616656277302155, -122.31156460382837) insert_cell_ids("Seattle Sunshine", 47.67366852914391, -122.29051997415843) insert_cell_ids("Sip House", 47.6682364706238, -122.31328618043693) insert_cell_ids("Victoria Coffee",47.624408595334536, -122.3117362652041) # Close connection conn.close()

На данный момент у нас есть база данных, в которой хранятся кофейни вместе с их идентификаторами ячеек, определяемыми выбранным разрешением для уровня ячейки.

Запрос ближайших кофеен

Наконец, давайте запросим кафе в районе Университетского округа.

 import sqlite3 from s2 import S2RegionCoverer,S2LatLngRect, S2LatLng # Connect to SQLite database conn = sqlite3.connect('/tmp/sqlite_cells.db') cursor = conn.cursor() # Function to query database for cells intersecting with the given polygon def query_intersecting_cells(start_x,start_y,end_x,end_y): # Create S2RegionCoverer region_rect = S2LatLngRect( S2LatLng.FromDegrees(start_x,start_y), S2LatLng.FromDegrees(end_x,end_y)) coverer = S2RegionCoverer() coverer.set_min_level(8) coverer.set_max_level(15) covering = coverer.GetCovering(region_rect) # Query for intersecting cells intersecting_cells = set() for cell_id in covering: cursor.execute("SELECT name FROM CoffeeShops WHERE cell_id >= ? and cell_id<=?", (str(cell_id.range_min()),str(cell_id.range_max()),)) intersecting_cells.update(cursor.fetchall()) return intersecting_cells # Query for intersecting cells intersecting_cells = query_intersecting_cells(47.6527847,-122.3286438,47.6782181, -122.2797203) # Print intersecting cells print("Intersecting cells:") for cell_id in intersecting_cells: print(cell_id[0]) # Close connection conn.close()

 Output: Intersecting cells: Sip House Seattle Sunshine

Ниже приведено визуальное представление ячеек. Для краткости приведенный ниже код визуализации не добавляется.

Поскольку все дочерние и родительские ячейки имеют общий префикс, мы можем запросить диапазоны ячеек между минимальным и максимальным значениями, чтобы получить все ячейки между этими двумя значениями. В нашем примере мы используем тот же принцип для запроса кафе.

Заключение:

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

Обычно для запросов к базе данных с геоиндексацией вам придется выполнить некоторую дополнительную постобработку данных. Чтобы гарантировать, что мы не перегружаем узел, необходимо тщательно продумать логику запросов и постобработки.

Использованная литература:

• Плоттер Google — https://github.com/gmplot/gmplot/wiki/GoogleMapPlotter
• Руководство разработчика S2 — http://s2geometry.io/devguide/
• Sqlite – https://docs.python.org/3/library/sqlite3.html.
• Геопространственные данные Azure Cosmos DB — https://learn.microsoft.com/en-us/azure/cosmos-db/nosql/query/geospatial?tabs=javascript
• ГДАЛ - https://gdal.org/index.html
• Шепли- https://shapely.readthedocs.io/en/stable/manual.html
• Геопанды - https://geopandas.org
• Посгис - https://postgis.net/
• Гаечный ключ Google — https://cloud.google.com/spanner?hl=en .
• Редис — https://redis.io/