diff --git a/.gitignore b/.gitignore index 3ca0331..3687264 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ *.parquet *.pbf -*.png \ No newline at end of file +*.png +*.pyc \ No newline at end of file diff --git a/.justfile b/.justfile index f553331..446a273 100644 --- a/.justfile +++ b/.justfile @@ -4,8 +4,8 @@ extract input output: download region output: quackosm --geom-filter-geocode {{region}} --osm-extract-source Geofabrik --osm-tags-filter-file tags.json --compact --output {{output}} -convert input output: - python3 slice.py --input {{input}} --output {{output}} +convert input output export_flags='--parquet': + python3 slice.py --input {{input}} --output {{output}} {{export_flags}} plot file: #!/usr/bin/env bash @@ -26,3 +26,6 @@ lookup lat lon: #!/usr/bin/env bash lat=$(echo "{{lat}}" | tr ',' ' ') python3 lookup.py $lat {{lon}} + +serve: + uvicorn api:app --host 0.0.0.0 --port 8000 --reload \ No newline at end of file diff --git a/api.py b/api.py new file mode 100644 index 0000000..c6d6c8f --- /dev/null +++ b/api.py @@ -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'.") \ No newline at end of file diff --git a/install.txt b/install.txt index 1f37aec..b806f52 100644 --- a/install.txt +++ b/install.txt @@ -1,2 +1,2 @@ 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 diff --git a/readme.md b/readme.md index f75358f..e70974d 100644 --- a/readme.md +++ b/readme.md @@ -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 ### Extract and Filter Data @@ -131,6 +147,6 @@ Sample Output: - [X] GeoHash parquet file generation - [X] Append data to existing parquet files -- [ ] Create web API server +- [X] Create web API server - [ ] Dockerize it - [ ] Scheduled update routine \ No newline at end of file diff --git a/slice.py b/slice.py index dce9115..eb5db5d 100644 --- a/slice.py +++ b/slice.py @@ -5,6 +5,8 @@ import geopandas as gpd import pandas as pd import pyarrow.parquet as pq import pyarrow as pa +import json +import cbor2 from pathlib import Path from shapely import wkb 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.""" if file_path.exists(): try: - # ✅ Load existing data existing_data = gpd.read_parquet(file_path) - - # ✅ Merge with new data 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: combined_data = combined_data.drop_duplicates(subset="feature_id", keep="last") - # ✅ Reset index before saving combined_data = combined_data.reset_index(drop=True) - - # ✅ Save back to the same file combined_data.to_parquet(file_path, index=False) print(f"✅ Appended new data to {file_path}") - except Exception as e: print(f"❌ Error while merging {file_path}: {e}") - else: - # If the file doesn't exist, create it new_data.to_parquet(file_path, index=False) print(f"✅ Created new file: {file_path}") -def slice_and_split_geoparquet(input_file: str, output_dir: str): - """ - Splits the GeoParquet file into multiple files based on 4-character Geohash boxes. - Appends new data to existing Geohash Parquet files while preventing duplicate geometries. - """ +def save_as_geojson(geohash_code: str, gdf: gpd.GeoDataFrame, output_path: Path): + """Saves a GeoDataFrame as a GeoJSON file.""" + geojson_file = output_path / f"{geohash_code}.geojson" + 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) if not input_path.exists(): 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.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'") - # Compute bounding box for each geometry bbox_results = con.execute(""" SELECT feature_id, rowid, ST_XMin(ST_Envelope(geometry)), ST_YMin(ST_Envelope(geometry)), ST_XMax(ST_Envelope(geometry)), ST_YMax(ST_Envelope(geometry)) FROM geoparquet; """).fetchall() - # Compute all intersecting geohashes geohash_mapping: Dict[str, List[int]] = {} 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].append(rowid) - # Dictionary to store in-memory results before writing geohash_data: Dict[str, gpd.GeoDataFrame] = {} for geohash_code, rowids in geohash_mapping.items(): geohash_bbox = geohash.bbox(geohash_code) - # Clip geometries inside this geohash box filtered_data = con.execute(f""" SELECT feature_id, tags, ST_AsWKB(ST_Intersection(geometry, @@ -110,46 +112,51 @@ def slice_and_split_geoparquet(input_file: str, output_dir: str): """).fetchdf() if not filtered_data.empty: - # ✅ Remove non-WKB values filtered_data = filtered_data[filtered_data["clipped_geom"].apply(lambda x: isinstance(x, (bytes, bytearray)))] if not filtered_data.empty: - # ✅ Convert WKB to Shapely geometries safely filtered_data["geometry"] = filtered_data["clipped_geom"].apply( lambda x: wkb.loads(bytes(x)) if isinstance(x, (bytes, bytearray)) else None ) filtered_data.drop(columns=["clipped_geom"], inplace=True) - # Remove rows where geometry conversion failed filtered_data = filtered_data.dropna(subset=["geometry"]) if not filtered_data.empty: - # Convert to GeoPandas GeoDataFrame gdf = gpd.GeoDataFrame(filtered_data, geometry="geometry", crs="EPSG:4326") - - # Store in memory 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(): - geohash_file = output_path / f"{geohash_code}.parquet" - append_to_parquet(geohash_file, gdf) - print(f"✅ Updated: {geohash_file}") + if export_parquet: + geohash_file = output_path / f"{geohash_code}.parquet" + 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() print("🎉 Processing complete!") 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("-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() + if not (args.parquet or args.geojson or args.cbor): + args.parquet = True # Default to GeoParquet if no options are given + 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: print(f"❌ Error: {e}") if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/tags.json b/tags.json index 49ca31e..ffbc300 100644 --- a/tags.json +++ b/tags.json @@ -1,5 +1,4 @@ { - "highway": ["motorway", "trunk", "primary"], "aeroway": true, "boundary": ["aboriginal_lands", "border_zone", "forest", "hazard", "national_park", "protected_area", "disputed"], "geological": true, diff --git a/u33d.png b/u33d.png index be96a0f..c4dae55 100644 Binary files a/u33d.png and b/u33d.png differ