API, GeoJSON, CBOR Support

This commit is contained in:
randogoth 2025-02-27 10:33:19 +00:00
parent 08eca3e042
commit a9f98c80df
8 changed files with 113 additions and 40 deletions

1
.gitignore vendored
View file

@ -1,3 +1,4 @@
*.parquet *.parquet
*.pbf *.pbf
*.png *.png
*.pyc

View file

@ -4,8 +4,8 @@ extract input output:
download region output: download region output:
quackosm --geom-filter-geocode {{region}} --osm-extract-source Geofabrik --osm-tags-filter-file tags.json --compact --output {{output}} quackosm --geom-filter-geocode {{region}} --osm-extract-source Geofabrik --osm-tags-filter-file tags.json --compact --output {{output}}
convert input output: convert input output export_flags='--parquet':
python3 slice.py --input {{input}} --output {{output}} python3 slice.py --input {{input}} --output {{output}} {{export_flags}}
plot file: plot file:
#!/usr/bin/env bash #!/usr/bin/env bash
@ -26,3 +26,6 @@ lookup lat lon:
#!/usr/bin/env bash #!/usr/bin/env bash
lat=$(echo "{{lat}}" | tr ',' ' ') lat=$(echo "{{lat}}" | tr ',' ' ')
python3 lookup.py $lat {{lon}} python3 lookup.py $lat {{lon}}
serve:
uvicorn api:app --host 0.0.0.0 --port 8000 --reload

47
api.py Normal file
View file

@ -0,0 +1,47 @@
from fastapi import FastAPI, HTTPException
from fastapi.responses import Response, FileResponse
from fastapi.middleware.gzip import GZipMiddleware
import geopandas as gpd
import cbor2
from pathlib import Path
app = FastAPI()
# ✅ Add Gzip compression middleware
app.add_middleware(GZipMiddleware, minimum_size=500)
# Directory containing Parquet files
GEOHASH_FOLDER = Path("geohash")
@app.get("/get_geodata/")
async def get_geodata(geohash: str, format: str = "cbor"):
"""Fetch geospatial data by Geohash and return in the requested format with HTTP compression."""
if len(geohash) < 4:
raise HTTPException(status_code=400, detail="Geohash must be at least 4 characters long.")
base_geohash = geohash[:4]
parquet_file = GEOHASH_FOLDER / f"{base_geohash}.parquet"
if not parquet_file.exists():
raise HTTPException(status_code=404, detail=f"No data found for geohash '{base_geohash}'.")
if format == "parquet":
return FileResponse(parquet_file, media_type="application/octet-stream",
filename=f"{base_geohash}.parquet")
gdf = gpd.read_parquet(parquet_file)
if format == "geojson":
geojson_data = gdf.to_json()
return Response(content=geojson_data, media_type="application/geo+json",
headers={"Content-Disposition": f"attachment; filename={base_geohash}.geojson"})
elif format == "cbor":
geojson_dict = gdf.to_json()
cbor_data = cbor2.dumps(geojson_dict)
return Response(content=cbor_data, media_type="application/cbor",
headers={"Content-Disposition": f"attachment; filename={base_geohash}.cbor"})
raise HTTPException(status_code=400, detail="Invalid format. Use 'parquet', 'geojson', or 'cbor'.")

View file

@ -1,2 +1,2 @@
sudo dnf install osmium-tool jq python3-devel pip sudo dnf install osmium-tool jq python3-devel pip
pip install quackosm[cli] psycopg2-binary pyosmium geoalchemy2 geojson shapely pygeohash pyarrow pandas click pip install quackosm[cli] psycopg2-binary pyosmium geoalchemy2 geojson shapely pygeohash pyarrow pandas click fastapi uvicorn

View file

@ -90,6 +90,22 @@ Sample Output: `just lookup 52.496846793890256 13.435128880554208`
] ]
``` ```
### 4. Web API
API on port `8000` with endpoint `get_geodata` that takes the parameter `geojson` (min 4 digits) and `format` (`parquet|geojson|cbor`) and returns a file in the requested format. Supports optional Gzip compression during transfer.
Start server:
```bash
just serve
```
Request example:
```bash
curl -H "Accept-Encoding: gzip" -o geodata.cbor.gz -v http://localhost:8000/get_geodata/?geohash=svc5&format=geojson
```
## Advanced Use ## Advanced Use
### Extract and Filter Data ### Extract and Filter Data
@ -131,6 +147,6 @@ Sample Output:
- [X] GeoHash parquet file generation - [X] GeoHash parquet file generation
- [X] Append data to existing parquet files - [X] Append data to existing parquet files
- [ ] Create web API server - [X] Create web API server
- [ ] Dockerize it - [ ] Dockerize it
- [ ] Scheduled update routine - [ ] Scheduled update routine

View file

@ -5,6 +5,8 @@ import geopandas as gpd
import pandas as pd import pandas as pd
import pyarrow.parquet as pq import pyarrow.parquet as pq
import pyarrow as pa import pyarrow as pa
import json
import cbor2
from pathlib import Path from pathlib import Path
from shapely import wkb from shapely import wkb
from typing import List, Dict from typing import List, Dict
@ -29,36 +31,41 @@ def append_to_parquet(file_path: Path, new_data: gpd.GeoDataFrame):
"""Appends new data to an existing Parquet file while ensuring correct merging.""" """Appends new data to an existing Parquet file while ensuring correct merging."""
if file_path.exists(): if file_path.exists():
try: try:
# ✅ Load existing data
existing_data = gpd.read_parquet(file_path) existing_data = gpd.read_parquet(file_path)
# ✅ Merge with new data
combined_data = pd.concat([existing_data, new_data], ignore_index=True) combined_data = pd.concat([existing_data, new_data], ignore_index=True)
# ✅ Drop duplicates *only if feature_id exists*
if "feature_id" in combined_data.columns: if "feature_id" in combined_data.columns:
combined_data = combined_data.drop_duplicates(subset="feature_id", keep="last") combined_data = combined_data.drop_duplicates(subset="feature_id", keep="last")
# ✅ Reset index before saving
combined_data = combined_data.reset_index(drop=True) combined_data = combined_data.reset_index(drop=True)
# ✅ Save back to the same file
combined_data.to_parquet(file_path, index=False) combined_data.to_parquet(file_path, index=False)
print(f"✅ Appended new data to {file_path}") print(f"✅ Appended new data to {file_path}")
except Exception as e: except Exception as e:
print(f"❌ Error while merging {file_path}: {e}") print(f"❌ Error while merging {file_path}: {e}")
else: else:
# If the file doesn't exist, create it
new_data.to_parquet(file_path, index=False) new_data.to_parquet(file_path, index=False)
print(f"✅ Created new file: {file_path}") print(f"✅ Created new file: {file_path}")
def slice_and_split_geoparquet(input_file: str, output_dir: str): def save_as_geojson(geohash_code: str, gdf: gpd.GeoDataFrame, output_path: Path):
""" """Saves a GeoDataFrame as a GeoJSON file."""
Splits the GeoParquet file into multiple files based on 4-character Geohash boxes. geojson_file = output_path / f"{geohash_code}.geojson"
Appends new data to existing Geohash Parquet files while preventing duplicate geometries. geojson_dict = json.loads(gdf.to_json())
"""
with open(geojson_file, "w", encoding="utf-8") as f:
json.dump(geojson_dict, f)
print(f"✅ Saved GeoJSON: {geojson_file}")
def save_as_cbor(geohash_code: str, gdf: gpd.GeoDataFrame, output_path: Path):
"""Saves a GeoDataFrame as a CBOR file (compact binary format)."""
cbor_file = output_path / f"{geohash_code}.cbor"
geojson_dict = json.loads(gdf.to_json())
with open(cbor_file, "wb") as f:
cbor2.dump(geojson_dict, f)
print(f"✅ Saved CBOR: {cbor_file}")
def slice_and_split_geoparquet(input_file: str, output_dir: str, export_parquet: bool, export_geojson: bool, export_cbor: bool):
"""Splits the GeoParquet file into multiple files based on 4-character Geohash tiles."""
input_path, output_path = Path(input_file), Path(output_dir) input_path, output_path = Path(input_file), Path(output_dir)
if not input_path.exists(): if not input_path.exists():
raise FileNotFoundError(f"File '{input_file}' not found.") raise FileNotFoundError(f"File '{input_file}' not found.")
@ -69,17 +76,14 @@ def slice_and_split_geoparquet(input_file: str, output_dir: str):
con = duckdb.connect() con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial;") con.execute("INSTALL spatial; LOAD spatial;")
# Create a temporary table of all non POINT geometries
con.execute(f"CREATE TEMP TABLE geoparquet AS SELECT * FROM read_parquet('{input_path}') WHERE ST_GeometryType(geometry) IS NOT NULL AND ST_GeometryType(geometry) != 'POINT'") con.execute(f"CREATE TEMP TABLE geoparquet AS SELECT * FROM read_parquet('{input_path}') WHERE ST_GeometryType(geometry) IS NOT NULL AND ST_GeometryType(geometry) != 'POINT'")
# Compute bounding box for each geometry
bbox_results = con.execute(""" bbox_results = con.execute("""
SELECT feature_id, rowid, ST_XMin(ST_Envelope(geometry)), ST_YMin(ST_Envelope(geometry)), SELECT feature_id, rowid, ST_XMin(ST_Envelope(geometry)), ST_YMin(ST_Envelope(geometry)),
ST_XMax(ST_Envelope(geometry)), ST_YMax(ST_Envelope(geometry)) ST_XMax(ST_Envelope(geometry)), ST_YMax(ST_Envelope(geometry))
FROM geoparquet; FROM geoparquet;
""").fetchall() """).fetchall()
# Compute all intersecting geohashes
geohash_mapping: Dict[str, List[int]] = {} geohash_mapping: Dict[str, List[int]] = {}
for feature_id, rowid, min_x, min_y, max_x, max_y in bbox_results: for feature_id, rowid, min_x, min_y, max_x, max_y in bbox_results:
@ -89,13 +93,11 @@ def slice_and_split_geoparquet(input_file: str, output_dir: str):
geohash_mapping[geohash_code] = [] geohash_mapping[geohash_code] = []
geohash_mapping[geohash_code].append(rowid) geohash_mapping[geohash_code].append(rowid)
# Dictionary to store in-memory results before writing
geohash_data: Dict[str, gpd.GeoDataFrame] = {} geohash_data: Dict[str, gpd.GeoDataFrame] = {}
for geohash_code, rowids in geohash_mapping.items(): for geohash_code, rowids in geohash_mapping.items():
geohash_bbox = geohash.bbox(geohash_code) geohash_bbox = geohash.bbox(geohash_code)
# Clip geometries inside this geohash box
filtered_data = con.execute(f""" filtered_data = con.execute(f"""
SELECT feature_id, tags, SELECT feature_id, tags,
ST_AsWKB(ST_Intersection(geometry, ST_AsWKB(ST_Intersection(geometry,
@ -110,44 +112,49 @@ def slice_and_split_geoparquet(input_file: str, output_dir: str):
""").fetchdf() """).fetchdf()
if not filtered_data.empty: if not filtered_data.empty:
# ✅ Remove non-WKB values
filtered_data = filtered_data[filtered_data["clipped_geom"].apply(lambda x: isinstance(x, (bytes, bytearray)))] filtered_data = filtered_data[filtered_data["clipped_geom"].apply(lambda x: isinstance(x, (bytes, bytearray)))]
if not filtered_data.empty: if not filtered_data.empty:
# ✅ Convert WKB to Shapely geometries safely
filtered_data["geometry"] = filtered_data["clipped_geom"].apply( filtered_data["geometry"] = filtered_data["clipped_geom"].apply(
lambda x: wkb.loads(bytes(x)) if isinstance(x, (bytes, bytearray)) else None lambda x: wkb.loads(bytes(x)) if isinstance(x, (bytes, bytearray)) else None
) )
filtered_data.drop(columns=["clipped_geom"], inplace=True) filtered_data.drop(columns=["clipped_geom"], inplace=True)
# Remove rows where geometry conversion failed
filtered_data = filtered_data.dropna(subset=["geometry"]) filtered_data = filtered_data.dropna(subset=["geometry"])
if not filtered_data.empty: if not filtered_data.empty:
# Convert to GeoPandas GeoDataFrame
gdf = gpd.GeoDataFrame(filtered_data, geometry="geometry", crs="EPSG:4326") gdf = gpd.GeoDataFrame(filtered_data, geometry="geometry", crs="EPSG:4326")
# Store in memory
geohash_data[geohash_code] = gdf geohash_data[geohash_code] = gdf
# Step 4: Write each geohash's data **only once** and prevent duplicates
for geohash_code, gdf in geohash_data.items(): for geohash_code, gdf in geohash_data.items():
geohash_file = output_path / f"{geohash_code}.parquet" if export_parquet:
append_to_parquet(geohash_file, gdf) geohash_file = output_path / f"{geohash_code}.parquet"
print(f"✅ Updated: {geohash_file}") append_to_parquet(geohash_file, gdf)
if export_geojson:
save_as_geojson(geohash_code, gdf, output_path)
if export_cbor:
save_as_cbor(geohash_code, gdf, output_path)
con.close() con.close()
print("🎉 Processing complete!") print("🎉 Processing complete!")
def main(): def main():
parser = argparse.ArgumentParser(description="Split a GeoParquet file into Geohash tiles.") parser = argparse.ArgumentParser(description="Split a GeoParquet file into Geohash tiles with export format options.")
parser.add_argument("-i", "--input", required=True, help="Path to the input GeoParquet file.") parser.add_argument("-i", "--input", required=True, help="Path to the input GeoParquet file.")
parser.add_argument("-o", "--output", required=True, help="Directory to save output Geohash tiles.") parser.add_argument("-o", "--output", required=True, help="Directory to save output files.")
parser.add_argument("--parquet", action="store_true", help="Export as GeoParquet")
parser.add_argument("--geojson", action="store_true", help="Export as GeoJSON")
parser.add_argument("--cbor", action="store_true", help="Export as CBOR")
args = parser.parse_args() args = parser.parse_args()
if not (args.parquet or args.geojson or args.cbor):
args.parquet = True # Default to GeoParquet if no options are given
try: try:
slice_and_split_geoparquet(args.input, args.output) slice_and_split_geoparquet(args.input, args.output, args.parquet, args.geojson, args.cbor)
except Exception as e: except Exception as e:
print(f"❌ Error: {e}") print(f"❌ Error: {e}")

View file

@ -1,5 +1,4 @@
{ {
"highway": ["motorway", "trunk", "primary"],
"aeroway": true, "aeroway": true,
"boundary": ["aboriginal_lands", "border_zone", "forest", "hazard", "national_park", "protected_area", "disputed"], "boundary": ["aboriginal_lands", "border_zone", "forest", "hazard", "national_park", "protected_area", "disputed"],
"geological": true, "geological": true,

BIN
u33d.png

Binary file not shown.

Before

Width:  |  Height:  |  Size: 4 MiB

After

Width:  |  Height:  |  Size: 4 MiB

Before After
Before After