How to Perform GIS Computations on Non-GIS Databases

Introduction:

Occasionally, you may encounter the need to perform geospatial functions within your application, such as mapping user locations or analyzing geographic data. There are numerous language-specific libraries available for these tasks, such as GDAL, Shapely, and Geopandas for Python.

Alternatively, geospatial functionality can be implemented through databases; for example, you can use the PostGIS extension with a relational database like PostgreSQL, or leverage the native support for spatial data types in a distributed database like Azure CosmosDB.

However, if your backend storage, such as Redis or Google Spanner, doesn't natively support spatial queries and you need to handle large-scale geospatial queries, this article is tailored for you.

What Are My Options?

You can always build another microservice to handle spatial data, but this option often involves the overhead of maintaining an additional application. Another approach is to use geo-indexing libraries like S2 and H3. S2, developed by Google, is based on the Hilbert curve, while H3, developed by Uber, is based on a geodesic discrete global grid system. S2 and H3 share many similarities: both divide a given region into cells and use 64-bit integers to index these cells.

However, the main difference lies in the shape of the cells; S2 uses square-shaped cells, whereas H3 uses hexagon-shaped cells. For some applications, H3 might offer better performance. However, overall, either library should suffice. In this article, we will use S2, but you can perform similar functions using H3.

Basic Concepts of Google S2 Library

• Cells: S2 divides the sphere into cells, each with a unique 64-bit identifier.

• Cell Levels: The hierarchy allows different levels of detail, from large regions to small precise areas. Each level represents a different resolution:

  • Level 0: The largest cells, covering a significant portion of the Earth's surface.

  • Higher Levels: Cells are progressively subdivided into smaller quadrants. For instance, Level 1 cells are each divided into four Level 2 cells, and so on.

  • Resolution and Area: Higher levels correspond to finer resolutions and smaller cell areas. This hierarchy allows for precise indexing and querying at varying levels of detail.

The table below shows various cell levels along with their corresponding areas.

From the table provided, it's evident that you can achieve mapping precision down to 0.44 cm^2 using S2. Within each square of an S2 cell, there exists a child cell that shares the same parent, indicating a hierarchical structure. The level of the cell can be a static value (same level applied to all the cells) or can be dynamic where S2 decides what resolution works best.

Calculating Nearest Neighbors

Let's start with an example. Consider we are writing an application that provides proximity service-like features for the Seattle area. We want to return a list of coffee shops within the given vicinity. To perform these operations, we will divide this task into 4 sub-tasks:

• Loading Seattle map
• Visualize S2 cells on the Seattle map
• Store a few coffee shop locations in the database
• Query for nearest coffee shops

Loading Seattle Map

To load a Google map, we would use the gmplot library. This library requires a Google Maps API key to load. To generate the API key, follow the instructions here.

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

The above code generates a map.html file as shown below:

Visualize S2 Cells on the Seattle Map

Now that we have the map, let's draw some S2 cells for maps :

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

In the code above, we first center the Google Map plotter around the Seattle area. In S2RegionCoverer, we initialize the region coverer to have dynamic levels between a minimum level of 8 and a maximum level of 15. This allows S2 to dynamically fit all cells into specific cell sizes for the best fit. The GetCovering method returns the covering for a rectangle around the Seattle area.

Then, we iterate over each cell, calculating the vertices for the cells and plotting them on the map. We keep the count of cells generated to around 273. Finally, we plot the input rectangle in red. This code will plot the S2 cells on the Seattle map at /tmp/map.html, as shown below:

Store a Few Coffee Shop Locations In the Database

Let’s generate a database of coffee shops along with their S2 cell identifiers. You can store these cells in any database of your choice. For this tutorial, we will be using a SQLite data database. In the code sample below, we connect to the SQLite database to create a table CoffeeShops with 3 fields Id, name, and cell_id.

Similar to the previous example, we use S2RegionCoverer to calculate the cells but this time, we use a fixed level for plotting points. Finally, the calculated ID is converted to a string and stored in the database.

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

At this point, we have a database that stores coffee shops along with their cell IDs, determined by the selected resolution for the cell level.

Query for the Nearest Coffee Shops

Finally, let’s query for coffee shops in the University District area.

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

Below is a visual representation of the cells. To maintain brevity, the visualization code for below is not added.

Since all child and parent cells share a prefix, we can query for cell ranges between min and max to get all the cells between those two values. In our example, we use the same principle to query coffee shop

Conclusion:

In this article, we have demonstrated how to use geo-indexing to store and query geospatial data in databases that don't support geospatial querying. This can be further extended to multiple use cases such as calculating routing between 2 points or getting nearest neighbors.

Typically, for geo-indexed database querying, you would have to perform some additional post-processing on the data. Careful consideration for querying and post-processing logic is needed to ensure we don't overwhelm the node.

References:

• Google map plotter - https://github.com/gmplot/gmplot/wiki/GoogleMapPlotter
• S2 developer guide - http://s2geometry.io/devguide/
• Sqlite - https://docs.python.org/3/library/sqlite3.html
• Azure Cosmos DB geospatial - 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
• Posgis - https://postgis.net/
• Google spanner - https://cloud.google.com/spanner?hl=en
• Redis - https://redis.io/