Cómo realizar cálculos SIG en bases de datos que no son SIG

Introducción:

Ocasionalmente, puede encontrarse con la necesidad de realizar funciones geoespaciales dentro de su aplicación, como mapear ubicaciones de usuarios o analizar datos geográficos. Existen numerosas bibliotecas específicas de idiomas disponibles para estas tareas, como GDAL, Shapely y Geopandas para Python.

Alternativamente, la funcionalidad geoespacial se puede implementar a través de bases de datos; por ejemplo, puede usar la extensión PostGIS con una base de datos relacional como PostgreSQL o aprovechar la compatibilidad nativa con tipos de datos espaciales en una base de datos distribuida como Azure CosmosDB.

Sin embargo, si su almacenamiento backend, como Redis o Google Spanner, no admite consultas espaciales de forma nativa y necesita manejar consultas geoespaciales a gran escala, este artículo está diseñado para usted.

¿Cuáles son mis opciones?

Siempre puedes crear otro microservicio para manejar datos espaciales, pero esta opción a menudo implica la sobrecarga de mantener una aplicación adicional. Otro enfoque es utilizar bibliotecas de indexación geográfica como S2 y H3. S2, desarrollado por Google, se basa en la curva de Hilbert, mientras que H3, desarrollado por Uber, se basa en un sistema de red global discreto geodésico. S2 y H3 comparten muchas similitudes: ambos dividen una región determinada en celdas y utilizan números enteros de 64 bits para indexar estas celdas.

Sin embargo, la principal diferencia radica en la forma de las células; S2 usa celdas de forma cuadrada, mientras que H3 usa celdas de forma hexagonal. Para algunas aplicaciones, H3 podría ofrecer un mejor rendimiento. Sin embargo, en general, cualquiera de las bibliotecas debería ser suficiente. En este artículo usaremos S2, pero puedes realizar funciones similares usando H3.

Conceptos básicos de la biblioteca Google S2

• Celdas: S2 divide la esfera en celdas, cada una con un identificador único de 64 bits.

• Niveles de celda: la jerarquía permite diferentes niveles de detalle, desde grandes regiones hasta pequeñas áreas precisas. Cada nivel representa una resolución diferente:

  
  

  • Nivel 0: Las células más grandes, que cubren una porción importante de la superficie terrestre.

    
    

  • Niveles superiores: las células se subdividen progresivamente en cuadrantes más pequeños. Por ejemplo, las celdas del Nivel 1 se dividen cada una en cuatro celdas del Nivel 2, y así sucesivamente.

    
    

  • Resolución y área: los niveles más altos corresponden a resoluciones más finas y áreas de celda más pequeñas. Esta jerarquía permite una indexación y consultas precisas en distintos niveles de detalle.

    
    

La siguiente tabla muestra varios niveles de celda junto con sus áreas correspondientes.

De la tabla proporcionada, es evidente que puede lograr una precisión de mapeo de hasta 0,44 cm^2 usando S2. Dentro de cada cuadrado de una celda S2, existe una celda secundaria que comparte el mismo padre, lo que indica una estructura jerárquica. El nivel de la celda puede ser un valor estático (el mismo nivel se aplica a todas las celdas) o puede ser dinámico donde S2 decide qué resolución funciona mejor.

Calcular los vecinos más cercanos

Comencemos con un ejemplo. Considere que estamos escribiendo una aplicación que proporciona funciones similares a un servicio de proximidad para el área de Seattle. Queremos devolver una lista de cafeterías dentro de la vecindad indicada. Para realizar estas operaciones dividiremos esta tarea en 4 subtareas:

• Cargando mapa de Seattle
• Visualice las celdas S2 en el mapa de Seattle
• Almacene algunas ubicaciones de cafeterías en la base de datos.
• Consulta de cafeterías más cercanas

Cargando mapa de Seattle

Para cargar un mapa de Google, usaríamos la biblioteca gmplot. Esta biblioteca requiere una clave API de Google Maps para cargarse. Para generar la clave API, siga las instrucciones aquí .

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

El código anterior genera un archivo map.html como se muestra a continuación:

Visualice células S2 en el mapa de Seattle

Ahora que tenemos el mapa, dibujemos algunas celdas S2 para mapas:

 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

En el código anterior, primero centramos el trazador de Google Map alrededor del área de Seattle. En S2RegionCoverer , inicializamos el cubridor de región para que tenga niveles dinámicos entre un nivel mínimo de 8 y un nivel máximo de 15. Esto permite que S2 ajuste dinámicamente todas las celdas en tamaños de celda específicos para lograr el mejor ajuste. El método GetCovering devuelve la cobertura de un rectángulo alrededor del área de Seattle.

Luego, iteramos sobre cada celda, calculamos los vértices de las celdas y los trazamos en el mapa. Mantenemos el recuento de celdas generadas en alrededor de 273. Finalmente, trazamos el rectángulo de entrada en rojo. Este código trazará las celdas S2 en el mapa de Seattle en /tmp/map.html , como se muestra a continuación:

Almacene algunas ubicaciones de cafeterías en la base de datos

Generemos una base de datos de cafeterías junto con sus identificadores de celda S2. Puede almacenar estas celdas en cualquier base de datos de su elección. Para este tutorial, usaremos una base de datos de datos SQLite. En el siguiente código de ejemplo, nos conectamos a la base de datos SQLite para crear una tabla CoffeeShops con 3 campos Id , name y cell_id .

De manera similar al ejemplo anterior, usamos S2RegionCoverer para calcular las celdas, pero esta vez usamos un nivel fijo para trazar puntos. Finalmente, la ID calculada se convierte en una cadena y se almacena en la base de datos.

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

En este punto, tenemos una base de datos que almacena cafeterías junto con sus ID de celda, determinados por la resolución seleccionada para el nivel de celda.

Consulta de cafeterías más cercanas

Finalmente, busquemos cafeterías en el área del Distrito Universitario.

 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

A continuación se muestra una representación visual de las celdas. Para mantener la brevedad, no se agrega el código de visualización que se muestra a continuación.

Dado que todas las celdas secundarias y principales comparten un prefijo, podemos consultar rangos de celdas entre el mínimo y el máximo para obtener todas las celdas entre esos dos valores. En nuestro ejemplo, utilizamos el mismo principio para consultar la cafetería.

Conclusión:

En este artículo, hemos demostrado cómo utilizar la indexación geográfica para almacenar y consultar datos geoespaciales en bases de datos que no admiten consultas geoespaciales. Esto se puede ampliar aún más a múltiples casos de uso, como calcular la ruta entre dos puntos o encontrar los vecinos más cercanos.

Normalmente, para consultas de bases de datos geoindexadas, tendría que realizar un posprocesamiento adicional de los datos. Es necesario considerar cuidadosamente la lógica de consulta y posprocesamiento para garantizar que no abrumemos el nodo.

Referencias:

• Trazador de mapas de Google: https://github.com/gmplot/gmplot/wiki/GoogleMapPlotter
• Guía para desarrolladores de S2: http://s2geometry.io/devguide/
• SQLite: https://docs.python.org/3/library/sqlite3.html
• Geoespacial de 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
• Posgis- https://postgis.net/
• Llave inglesa de Google: https://cloud.google.com/spanner?hl=en
• Redis- https://redis.io/