Running SQL queries against your BDC Exports to check coverage

[Seth_Poche]

From time to time we receive requests from users who would like to programmatically check coverage for potential customers. This is possible by loading your BDC Exports into a Postgres (PostGIS) database and running SELECT queries with the lat/lon location(s) in question. The steps below define how this can be done.

(Note: This installation considers Windows OS.)

1. Run BDC exports using cnHeat and extract the results.
   

   

2. Install Postgres 17
   a. EDB: Open-Source, Enterprise Postgres Database Management
   b. Note: Remember to record the password during installation.

3. At the end of installation when prompted, use Stack Builder to install PostGIS.

4. Login to your database with the psql command prompt and enable PostGIS by running the following commands.
   
   

Enable Post GIS by copy/pasting these commands into the psql command prompt window. Enabling PostGIS | PostGIS

5. Install OSGeo4W which will be needed to import the .gpkg export file directly into your database.
   a. OSGeo4W
   
   

6. Pull up the OSGeo4W Shell and import the BDC exports into your database.
   

Example:

ogr2ogr -f PostgreSQL “PG:dbname=postgres user=postgres password=postgres” “c:$tmp\broadband_coverage.gpkg”

7. From there you can query against the database with the following statement.
   In this example we use pgAdmin 4 which ships with Postgres.
   Query: select * from merged where ST_Covers(geom, ST_Point(-98.586602, 33.852654, 4326)::geography);
   
   

From the results above we can see that we can provide 25/3 Mbps service to this location using our “Licensed by Rule” technology (CBRS-GAA).

For a more polished approach consider using dBeaver for your database client which has built-in geospatial functionality.

8. The last step (which isn’t fully covered here) is to query the database programmatically. This process is well documented on the Internet and will differ for each operator.

   Example: PostgreSQL Python: Querying Data

[Chris_Josey]

Cool! For those of you that prefer a way to do all of this with just Python, here you go.

Edit: Made it so that it returns the polygon info that the point falls in.

### Need to install the shapely and geopandas modules
### pip install shapely
### pip install geopandas

from shapely.geometry import Point
import geopandas as gpd


def check_coverage(latitude, longitude, coverage_file):
    """
    Checks if a given latitude and longitude coordinate falls within a polygon defined in a GeoPackage (.geopkg) file.
    If the point falls within a polygon, the function returns the latitude, longitude, result, and the information 
    of the polygon it falls within.

    Parameters:
    ----------
    latitude : float
        Latitude of the point to be checked.

    longitude : float
        Longitude of the point to be checked.

    coverage_file : str
        The file path to the GeoPackage (.geopkg) file containing the polygon geometry.

    Returns:
    -------
    dict
        A dictionary with the latitude, longitude, result string, and the attributes of the polygon the point falls 
        within. If no polygon is found, returns a dictionary indicating the point is "NOT In Coverage".

    Example:
    -------
    result = check_coverage(29.651634, -82.324829, "coverage_map.geopkg")
    print(result)
    # Output: {'latitude': 29.651634, 'longitude': -82.324829, 'result': 'In Coverage', 'polygon_info': {...}}

    Raises:
    ------
    FileNotFoundError
        If the .geopkg file specified by the coverage_file does not exist or cannot be opened.
    
    TypeError
        If the input latitude or longitude is not a valid float.
    """
    
    try:
        point = Point(float(longitude), float(latitude))
    except ValueError:
        raise TypeError("Latitude and longitude must be valid float values.")
    
    try:
        geopackage = gpd.read_file(coverage_file)
    except FileNotFoundError:
        raise FileNotFoundError(f"File {coverage_file} not found.")
    
    for _, row in geopackage.iterrows():
        polygon = row['geometry']
        if point.within(polygon):
            return {
                "latitude": latitude,
                "longitude": longitude,
                "result": "In Coverage",
                "polygon_info": row.to_dict()
            }
    
    return {"latitude": latitude, "longitude": longitude, "result": "NOT In Coverage"}