92 lines
No EOL
3.3 KiB
Python
92 lines
No EOL
3.3 KiB
Python
import duckdb
|
||
from pathlib import Path
|
||
import geohash
|
||
import osmium
|
||
import subprocess
|
||
from pprint import pprint
|
||
from typing import List, Tuple, Dict
|
||
|
||
def remove_points(input_file: str, output_file: str):
|
||
try:
|
||
input_path, output_path = Path(input_file), Path(output_file)
|
||
if not input_path.exists():
|
||
raise FileNotFoundError(f"File '{input_file}' not found.")
|
||
|
||
con = duckdb.connect()
|
||
con.execute("INSTALL spatial; LOAD spatial;")
|
||
|
||
con.execute(f"""
|
||
CREATE VIEW geoparquet_filtered AS
|
||
SELECT * FROM read_parquet('{input_path}')
|
||
WHERE ST_GeometryType(geometry) IS NOT NULL
|
||
AND ST_GeometryType(geometry) != 'POINT'
|
||
""")
|
||
|
||
if con.execute("SELECT COUNT(*) FROM geoparquet_filtered").fetchone()[0] == 0:
|
||
raise ValueError("No non-POINT geometries found.")
|
||
|
||
output_path.parent.mkdir(parents=True, exist_ok=True)
|
||
con.execute(f"COPY geoparquet_filtered TO '{output_path}' (FORMAT PARQUET)")
|
||
print(f"✅ Saved to {output_path}")
|
||
|
||
except Exception as e:
|
||
print(f"❌ Error: {e}")
|
||
|
||
finally:
|
||
if 'con' in locals():
|
||
con.close()
|
||
|
||
def bbox(file_path: str) -> Dict[str, float]:
|
||
try:
|
||
result = subprocess.run(
|
||
["osmium", "fileinfo", "-e", file_path],
|
||
capture_output=True, text=True, check=True
|
||
)
|
||
for line in result.stdout.split("\n"):
|
||
if "Bounding box:" in line:
|
||
bbox_str = line.strip().replace("Bounding box:", "").strip()
|
||
# Convert to a tuple of floats
|
||
bbox_tuple = tuple(map(float, bbox_str.strip("()").split(",")))
|
||
return bbox_tuple
|
||
except subprocess.CalledProcessError as e:
|
||
print(f"❌ Error running osmium: {e}")
|
||
|
||
def check_data(file: str):
|
||
con = duckdb.connect()
|
||
con.execute("INSTALL spatial; LOAD spatial;")
|
||
con.execute(f"CREATE VIEW geoparquet AS SELECT * FROM read_parquet('{file}')")
|
||
columns = con.execute("DESCRIBE geoparquet").fetchdf()
|
||
print("Columns in the file:\n", columns)
|
||
dc = con.execute("SELECT count(*) FROM geoparquet LIMIT 10").fetchdf()
|
||
print("\nTotal rows:", dc)
|
||
df = con.execute("SELECT * FROM geoparquet LIMIT 10").fetchdf()
|
||
print("\nFirst 10 rows:\n", df)
|
||
try:
|
||
geo_types = con.execute("SELECT DISTINCT ST_GeometryType(geometry) AS geom_type FROM geoparquet").fetchdf()
|
||
print("\nDetected Geometry Types:\n", geo_types)
|
||
except Exception as e:
|
||
print("\nCould not determine geometry type. Possible reason: No valid geometry column found.\n", e)
|
||
con.close()
|
||
|
||
def bbox_geohashes(bbox: Tuple[float, float, float, float]) -> List[str]:
|
||
min_x, min_y, max_x, max_y = bbox
|
||
geohashes = set()
|
||
|
||
# Define step size based on 4-character Geohash precision (~39km × 19.5km)
|
||
step_size = 0.25 # Approximate step size in degrees (varies slightly by latitude)
|
||
|
||
lat = min_y
|
||
while lat <= max_y:
|
||
lon = min_x
|
||
while lon <= max_x:
|
||
gh = geohash.encode(lat, lon, precision=4)
|
||
geohashes.add(gh)
|
||
lon += step_size # Move eastward
|
||
lat += step_size # Move northward
|
||
|
||
return sorted(geohashes)
|
||
|
||
remove_points('compact.parquet', 'processed.parquet')
|
||
bbox = bbox('israel.pbf')
|
||
hashes = bbox_geohashes(bbox)
|
||
print(hashes) |