From 513fbc4fc300ad9a25d039a88852a1d91d3fdbf6 Mon Sep 17 00:00:00 2001 From: randogoth Date: Tue, 4 Mar 2025 08:07:21 +0000 Subject: [PATCH] all water coasline --- water.py | 237 ++++++++++++++++++++++--------------------------------- 1 file changed, 94 insertions(+), 143 deletions(-) diff --git a/water.py b/water.py index 9a0dd24..053ec9f 100644 --- a/water.py +++ b/water.py @@ -1,168 +1,119 @@ +import os import argparse -import duckdb -import geohash import geopandas as gpd -import pandas as pd -import json +import pygeohash as pgh import cbor2 -from pathlib import Path -from shapely import wkb -from typing import List, Dict +import pandas as pd +import logging +from shapely.geometry import box, Polygon +from tqdm import tqdm -def get_geohashes_from_bbox(min_x, min_y, max_x, max_y) -> List[str]: - """Generates all 4-character Geohashes that intersect a bounding box.""" - geohashes = set() - step_size = 0.25 # 4-character Geohash grid size +# Configure logging +logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.INFO) +logger = logging.getLogger(__name__) - 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 - lat += step_size +def geohash_bbox(geohash_code): + """Returns a Shapely Polygon representing the bounding box of a geohash.""" + lat, lon, lat_err, lon_err = pgh.decode_exactly(geohash_code) + lat_min, lon_min = lat - lat_err, lon - lon_err + lat_max, lon_max = lat + lat_err, lon + lon_err + return box(lon_min, lat_min, lon_max, lat_max) - return sorted(geohashes) +def generate_4digit_geohashes(): + """Generates all possible 4-digit geohashes covering the world.""" + base_geohashes = "0123456789bcdefghjkmnpqrstuvwxyz" + return [a + b + c + d for a in base_geohashes for b in base_geohashes for c in base_geohashes for d in base_geohashes] -def append_to_parquet(file_path: Path, new_data: gpd.GeoDataFrame): - """Appends new data to an existing Parquet file while ensuring correct merging and removing duplicates.""" - if file_path.exists(): - try: - existing_data = gpd.read_parquet(file_path) - combined_data = pd.concat([existing_data, new_data], ignore_index=True) - combined_data = combined_data.drop_duplicates(subset=["geometry"], keep="last").reset_index(drop=True) - combined_data.to_parquet(file_path, index=False) - print(f"✅ Appended new data to {file_path}, duplicates removed") - except Exception as e: - print(f"❌ Error while merging {file_path}: {e}") - else: - new_data.to_parquet(file_path, index=False) - print(f"✅ Created new file: {file_path}") +def intersect_geohash_tiles(water_polygons, coastline, geohashes): + """Finds and slices water polygons into 4-digit geohash tiles along the coastline.""" + results = {} -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" - gdf.to_file(geojson_file, driver="GeoJSON") - 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, coastline_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, keeping only coast-adjacent areas.""" - input_path, coastline_path, output_path = Path(input_file), Path(coastline_file), Path(output_dir) - if not input_path.exists() or not coastline_path.exists(): - raise FileNotFoundError(f"File '{input_file}' or '{coastline_file}' not found.") - output_path.mkdir(parents=True, exist_ok=True) - - print(f"📂 Processing '{input_file}', output will be saved in '{output_dir}'...") - - con = duckdb.connect() - con.execute("INSTALL spatial; LOAD spatial;") - - # Load coastline data and create a 30km buffer - con.execute(f""" - CREATE TEMP TABLE coastline AS - SELECT ST_Buffer(geometry, 30000) AS buffer_geom - FROM read_parquet('{coastline_path}'); - """) - - 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'") - - bbox_results = con.execute(""" - SELECT rowid, ST_XMin(ST_Envelope(geometry)), ST_YMin(ST_Envelope(geometry)), - ST_XMax(ST_Envelope(geometry)), ST_YMax(ST_Envelope(geometry)) - FROM geoparquet; - """).fetchall() - - geohash_mapping: Dict[str, List[int]] = {} - - for rowid, min_x, min_y, max_x, max_y in bbox_results: - intersecting_geohashes = get_geohashes_from_bbox(min_x, min_y, max_x, max_y) + logger.info("Processing geohashes for intersection with water polygons...") + for ghash in tqdm(geohashes, desc="Processing geohashes"): + ghash_poly = geohash_bbox(ghash) - # Check if geohash is near coastline buffer - valid_geohashes = [] - for gh in intersecting_geohashes: - gh_bbox = geohash.bbox(gh) - intersects = con.execute(f""" - SELECT COUNT(*) FROM coastline - WHERE ST_Intersects(buffer_geom, - ST_MakeEnvelope({gh_bbox['w']}, {gh_bbox['s']}, {gh_bbox['e']}, {gh_bbox['n']})); - """).fetchone()[0] > 0 - - if intersects: - valid_geohashes.append(gh) + # Select only water polygons that intersect with the geohash tile + water_subset = water_polygons[water_polygons.intersects(ghash_poly)] + if water_subset.empty: + continue - for geohash_code in valid_geohashes: - if geohash_code not in geohash_mapping: - geohash_mapping[geohash_code] = [] - geohash_mapping[geohash_code].append(rowid) - - for geohash_code, rowids in geohash_mapping.items(): - geohash_bbox = geohash.bbox(geohash_code) + # Clip with geohash bounding box + clipped_water = gpd.clip(water_subset, ghash_poly) + + # Remove areas already covered by coastline + coastline_subset = coastline[coastline.intersects(ghash_poly)] + if not coastline_subset.empty: + clipped_water = clipped_water.overlay(coastline_subset, how="difference") - filtered_data = con.execute(f""" - SELECT ST_AsWKB(ST_Intersection(geometry, - ST_MakeEnvelope({geohash_bbox['w']}, {geohash_bbox['s']}, - {geohash_bbox['e']}, {geohash_bbox['n']}))) - AS clipped_geom - FROM geoparquet - WHERE rowid IN ({','.join(map(str, rowids))}) - AND ST_Intersects(geometry, - ST_MakeEnvelope({geohash_bbox['w']}, {geohash_bbox['s']}, - {geohash_bbox['e']}, {geohash_bbox['n']})); - """).fetchdf() + if not clipped_water.empty: + results[ghash] = clipped_water + logger.debug(f"Geohash {ghash} contains {len(clipped_water)} water polygons.") - if not filtered_data.empty: - filtered_data = filtered_data[filtered_data["clipped_geom"].apply(lambda x: isinstance(x, (bytes, bytearray)))] - - if not filtered_data.empty: - 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) + logger.info(f"Completed geohash processing. {len(results)} geohash tiles contain water polygons.") + return results - filtered_data = filtered_data.dropna(subset=["geometry"]) +def save_data(geohash_tiles, output_dir, export_cbor, export_geojson): + """Saves geohash-sliced water polygons in requested formats.""" + os.makedirs(output_dir, exist_ok=True) + logger.info(f"Saving data to {output_dir}...") - if not filtered_data.empty: - filtered_data["natural"] = "water" - filtered_data["water"] = "sea" - - gdf = gpd.GeoDataFrame(filtered_data, geometry="geometry", crs="EPSG:4326") + for ghash, gdf in tqdm(geohash_tiles.items(), desc="Saving geohashes"): + parquet_path = os.path.join(output_dir, f"{ghash}.parquet") + gdf.to_parquet(parquet_path) + logger.debug(f"Saved {ghash}.parquet") - if export_parquet: - append_to_parquet(output_path / f"{geohash_code}.parquet", gdf) - if export_geojson: - save_as_geojson(geohash_code, gdf, output_path) - if export_cbor: - save_as_cbor(geohash_code, gdf, output_path) + if export_cbor: + cbor_path = os.path.join(output_dir, f"{ghash}.cbor") + with open(cbor_path, "wb") as f: + cbor2.dump(gdf.to_dict(), f) + logger.debug(f"Saved {ghash}.cbor") - con.close() - print("🎉 Processing complete!") + if export_geojson: + geojson_path = os.path.join(output_dir, f"{ghash}.geojson") + gdf.to_file(geojson_path, driver="GeoJSON") + logger.debug(f"Saved {ghash}.geojson") + + logger.info("Data saving process completed.") def main(): - parser = argparse.ArgumentParser(description="Split a GeoParquet file into coastal Geohash tiles with export format options.") - parser.add_argument("-i", "--input", required=True, help="Path to the input GeoParquet file.") - parser.add_argument("-c", "--coastline", required=True, help="Path to the coastline Parquet file.") - 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") + parser = argparse.ArgumentParser(description="Slice water polygons along the coastline into 4-digit geohash tiles.") + parser.add_argument("coastline", help="Path to the coastline Parquet file.") + parser.add_argument("water", help="Path to the water polygons Parquet file.") + parser.add_argument("output_dir", help="Directory to save the output files.") + parser.add_argument("--cbor", action="store_true", help="Export as CBOR.") + parser.add_argument("--geojson", action="store_true", help="Export as GeoJSON.") + parser.add_argument("--verbose", action="store_true", help="Enable detailed logging.") 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 + # Set logging level + if args.verbose: + logger.setLevel(logging.DEBUG) + logger.debug("Verbose logging enabled.") - try: - slice_and_split_geoparquet(args.input, args.coastline, args.output, args.parquet, args.geojson, args.cbor) - except Exception as e: - print(f"❌ Error: {e}") + logger.info("Loading input datasets...") + coastline = gpd.read_parquet(args.coastline) + water_polygons = gpd.read_parquet(args.water) + + # Ensure the geometries are valid + coastline = coastline[coastline.geometry.is_valid] + water_polygons = water_polygons[water_polygons.geometry.is_valid] + + logger.info(f"Loaded coastline dataset with {len(coastline)} features.") + logger.info(f"Loaded water polygons dataset with {len(water_polygons)} features.") + + # Generate 4-digit geohashes + logger.info("Generating all possible 4-digit geohashes...") + geohashes = generate_4digit_geohashes() + logger.info(f"Generated {len(geohashes)} geohash tiles.") + + # Intersect water polygons with geohash tiles + geohash_tiles = intersect_geohash_tiles(water_polygons, coastline, geohashes) + + # Save results + save_data(geohash_tiles, args.output_dir, args.cbor, args.geojson) + + logger.info("Processing complete. All files saved.") if __name__ == "__main__": main()