如何在非 GIS 数据库上执行 GIS 计算

介绍：

有时，您可能会需要在应用程序中执行地理空间功能，例如映射用户位置或分析地理数据。有许多特定于语言的库可用于执行这些任务，例如 GDAL、Shapely 和 Python 的 Geopandas。

或者，可以通过数据库实现地理空间功能；例如，您可以将 PostGIS 扩展与 PostgreSQL 等关系数据库一起使用，或者利用 Azure CosmosDB 等分布式数据库对空间数据类型的本机支持。

但是，如果您的后端存储（例如 Redis 或 Google Spanner）本身不支持空间查询，并且您需要处理大规模地理空间查询，那么本文非常适合您。

我有哪些选择？

您始终可以构建另一个微服务来处理空间数据，但此选项通常涉及维护额外应用程序的开销。另一种方法是使用地理索引库，如 S2 和 H3。由 Google 开发的 S2 基于希尔伯特曲线，而由 Uber 开发的 H3 基于测地线离散全球网格系统。S2 和 H3 有许多相似之处：都将给定区域划分为单元并使用 64 位整数来索引这些单元。

但是，主要的区别在于单元格的形状；S2 使用方形单元格，而 H3 使用六边形单元格。对于某些应用程序，H3 可能会提供更好的性能。但是，总体而言，任何一个库都应该足够。在本文中，我们将使用 S2，但您可以使用 H3 执行类似的功能。

Google S2 库的基本概念

• 细胞：S2 将球体划分为细胞，每个细胞都有一个唯一的 64 位标识符。

• 单元级别：层次结构允许不同级别的细节，从大区域到小的精确区域。每个级别代表不同的分辨率：

  
  

  • 0 级：最大的细胞，覆盖地球表面的很大一部分。

    
    

  • 更高级别：单元格逐渐细分为更小的象限。例如，每个 1 级单元格被划分为四个 2 级单元格，依此类推。

    
    

  • 分辨率和面积：级别越高，分辨率越精细，单元格面积越小。此层次结构允许在不同细节级别上进行精确索引和查询。

    
    

下表显示了各个细胞级别及其对应的区域。

从提供的表格中可以看出，使用 S2 可以实现低至 0.44 cm^2 的映射精度。在 S2 单元格的每个方块内，都存在一个共享同一父级的子单元格，这表明存在层次结构。单元格的级别可以是静态值（同一级别应用于所有单元格），也可以是动态的，其中 S2 决定哪种分辨率效果最好。

计算最近邻居

让我们从一个例子开始。假设我们正在编写一个应用程序，为西雅图地区提供类似邻近服务的功能。我们想要返回给定区域内的咖啡店列表。为了执行这些操作，我们将此任务分为 4 个子任务：

• 正在加载西雅图地图
• 在西雅图地图上可视化 S2 单元格
• 在数据库中存储一些咖啡店的位置
• 查询最近的咖啡店

正在加载西雅图地图

要加载 Google 地图，我们将使用 gmplot 库。此库需要 Google 地图 API 密钥才能加载。要生成 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 地图绘图器置于西雅图地区的中心。在S2RegionCoverer中，我们初始化区域覆盖器，使其具有最低级别 8 和最高级别 15 之间的动态级别。这允许 S2 动态地将所有单元格调整GetCovering特定单元格大小，以实现最佳拟合。GetCovering 方法返回西雅图地区周围矩形的覆盖范围。

然后，我们迭代每个单元格，计算单元格的顶点并将它们绘制在地图上。我们将生成的单元格数量保持在 273 左右。最后，我们用红色绘制输入矩形。此代码将在西雅图地图/tmp/map.html上绘制 S2 单元格，如下所示：

在数据库中存储一些咖啡店位置

让我们生成一个咖啡店数据库及其 S2 单元标识符。您可以将这些单元存储在您选择的任何数据库中。在本教程中，我们将使用 SQLite 数据数据库。在下面的代码示例中，我们连接到 SQLite 数据库以创建一个表CoffeeShops其中包含 3 个字段Id 、 name和cell_id 。

与上一个示例类似，我们使用S2RegionCoverer来计算单元格，但这次我们使用固定级别来绘制点。最后，将计算出的 ID 转换为字符串并存储在数据库中。

 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()

此时，我们有一个数据库，其中存储了咖啡店及其单元 ID，由所选的单元级别分辨率决定。

查询最近的咖啡店

最后，我们来查询大学区内的咖啡店。

 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

以下是单元格的可视化表示。为了简洁起见，未添加下面的可视化代码。

由于所有子单元格和父单元格都共享一个前缀，因此我们可以查询最小值和最大值之间的单元格范围，以获取这两个值之间的所有单元格。在我们的示例中，我们使用相同的原理来查询咖啡店

结论：

在本文中，我们演示了如何使用地理索引在不支持地理空间查询的数据库中存储和查询地理空间数据。这可以进一步扩展到多种用例，例如计算两点之间的路线或获取最近邻居。

通常，对于地理索引数据库查询，您必须对数据执行一些额外的后处理。需要仔细考虑查询和后处理逻辑，以确保我们不会让节点不堪重负。

参考：

• 谷歌地图绘图仪 - 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
• GDAL- https://gdal.org/index.html
• Shapley- https://shapely.readthedocs.io/en/stable/manual.html
• Geopandas - https://geopandas.org
• Postgis - https://postgis.net/
• Google 扳手 - https://cloud.google.com/spanner?hl=en
• Redis—— https://redis.io/