proper coastline
This commit is contained in:
randogoth
parent
8009cc6846
commit
dd68ff0131
1 changed files with 70 additions and 144 deletions
206
water.py
206
water.py
|
|
@ -1,180 +1,106 @@
|
||||||
import os
|
import os
|
||||||
import argparse
|
import argparse
|
||||||
import geopandas as gpd
|
|
||||||
import pygeohash as pgh
|
|
||||||
import cbor2
|
|
||||||
import pandas as pd
|
import pandas as pd
|
||||||
|
import geopandas as gpd
|
||||||
|
import geohash
|
||||||
import logging
|
import logging
|
||||||
from shapely.geometry import box
|
from shapely.geometry import box
|
||||||
from shapely.ops import unary_union
|
from shapely.ops import unary_union
|
||||||
from tqdm import tqdm
|
|
||||||
|
|
||||||
# Configure logging
|
# Configure logging
|
||||||
logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.INFO)
|
logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.INFO)
|
||||||
logger = logging.getLogger(__name__)
|
logger = logging.getLogger(__name__)
|
||||||
|
|
||||||
def geohash_bbox(geohash_code):
|
def load_geohash_files(folder):
|
||||||
"""Returns a Shapely Polygon representing the bounding box of a geohash."""
|
"""Finds all existing geohash parquet files in the folder."""
|
||||||
lat, lon, lat_err, lon_err = pgh.decode_exactly(geohash_code)
|
geohash_files = {}
|
||||||
lat_min, lon_min = lat - lat_err, lon - lon_err
|
for filename in os.listdir(folder):
|
||||||
lat_max, lon_max = lat + lat_err, lon + lon_err
|
|
||||||
return box(lon_min, lat_min, lon_max, lat_max)
|
|
||||||
|
|
||||||
def get_existing_geohashes(output_dir):
|
|
||||||
"""Returns a set of existing geohash parquet files in the output directory."""
|
|
||||||
existing_geohashes = set()
|
|
||||||
for filename in os.listdir(output_dir):
|
|
||||||
if filename.endswith(".parquet"):
|
if filename.endswith(".parquet"):
|
||||||
geohash_code = filename[:-8] # Remove '.parquet' extension
|
geohash_code = filename[:-8] # Remove '.parquet' extension
|
||||||
existing_geohashes.add(geohash_code)
|
geohash_files[geohash_code] = os.path.join(folder, filename)
|
||||||
logger.info(f"Found {len(existing_geohashes)} existing geohash files in {output_dir}.")
|
logger.info(f"Found {len(geohash_files)} geohash parquet files in {folder}.")
|
||||||
return existing_geohashes
|
return geohash_files
|
||||||
|
|
||||||
def subtract_existing_islands(water_gdf, existing_gdf):
|
def geohash_bbox(geohash_code):
|
||||||
"""Subtracts intersecting land polygons (islands) from water polygons."""
|
"""Returns a bounding box polygon for a geohash."""
|
||||||
if existing_gdf.empty:
|
bbox = geohash.bbox(geohash_code)
|
||||||
return water_gdf # No islands to subtract
|
return box(bbox["w"], bbox["s"], bbox["e"], bbox["n"])
|
||||||
|
|
||||||
# Collect all existing non-water geometries (islands or coastlines)
|
def ensure_crs_consistency(gdf):
|
||||||
land_geometries = existing_gdf[existing_gdf["tags"].apply(lambda tags: tags and tags.get("natural") != "water")]
|
"""Ensures the GeoDataFrame is in EPSG:4326 to prevent CRS mismatches."""
|
||||||
|
if gdf.crs is None or gdf.crs.to_string() != "EPSG:4326":
|
||||||
|
logger.info(f"Converting CRS of {gdf} to EPSG:4326")
|
||||||
|
gdf = gdf.to_crs("EPSG:4326")
|
||||||
|
return gdf
|
||||||
|
|
||||||
if land_geometries.empty:
|
def process_geohash_files(coastline_gdf, geohash_files):
|
||||||
return water_gdf # No islands to subtract
|
"""Processes each geohash file: removes old coastlines and inserts updated ones."""
|
||||||
|
for geohash_code, file_path in geohash_files.items():
|
||||||
|
logger.info(f"Processing {file_path}...")
|
||||||
|
|
||||||
# Merge islands into a single geometry
|
# Load existing geohash file
|
||||||
islands_union = unary_union(land_geometries.geometry)
|
try:
|
||||||
|
existing_gdf = gpd.read_parquet(file_path)
|
||||||
# Subtract islands from the water polygons
|
except Exception as e:
|
||||||
water_gdf["geometry"] = water_gdf.geometry.difference(islands_union)
|
logger.warning(f"Failed to read {file_path}. Skipping. Error: {e}")
|
||||||
|
|
||||||
# Remove empty geometries resulting from subtraction
|
|
||||||
water_gdf = water_gdf[~water_gdf.geometry.is_empty]
|
|
||||||
|
|
||||||
return water_gdf
|
|
||||||
|
|
||||||
def intersect_geohash_tiles(water_polygons, coastline, geohashes):
|
|
||||||
"""Finds and slices water polygons into geohash tiles along the coastline."""
|
|
||||||
results = {}
|
|
||||||
|
|
||||||
logger.info("Processing geohashes for intersection with water polygons...")
|
|
||||||
for ghash in tqdm(geohashes, desc="Processing geohashes"):
|
|
||||||
ghash_poly = geohash_bbox(ghash)
|
|
||||||
|
|
||||||
# Select only water polygons that intersect with the geohash tile
|
|
||||||
water_subset = water_polygons[water_polygons.intersects(ghash_poly)]
|
|
||||||
if water_subset.empty:
|
|
||||||
continue
|
continue
|
||||||
|
|
||||||
# Clip with geohash bounding box
|
existing_gdf = ensure_crs_consistency(existing_gdf)
|
||||||
clipped_water = gpd.clip(water_subset, ghash_poly)
|
|
||||||
|
|
||||||
# Remove areas already covered by coastline
|
# Remove existing coastline features
|
||||||
coastline_subset = coastline[coastline.intersects(ghash_poly)]
|
filtered_gdf = existing_gdf[~((existing_gdf.geometry.type.isin(["LineString", "MultiLineString"])) &
|
||||||
if not coastline_subset.empty:
|
(existing_gdf["tags"].apply(lambda tags: isinstance(tags, dict) and tags.get("natural") == "coastline")))]
|
||||||
clipped_water = clipped_water.overlay(coastline_subset, how="difference")
|
|
||||||
|
|
||||||
# Add water-specific tags
|
# Get coastline geometries for this geohash
|
||||||
if not clipped_water.empty:
|
bbox = geohash_bbox(geohash_code)
|
||||||
clipped_water["tags"] = [{"natural": "water", "water": "sea"}] * len(clipped_water)
|
new_coastline = coastline_gdf[coastline_gdf.intersects(bbox)].copy()
|
||||||
clipped_water = clipped_water[["tags", "geometry"]] # Drop unnecessary columns
|
|
||||||
results[ghash] = clipped_water
|
|
||||||
|
|
||||||
logger.debug(f"Geohash {ghash} contains {len(clipped_water)} water polygons.")
|
# Clip coastline to the geohash boundary
|
||||||
|
new_coastline["geometry"] = new_coastline["geometry"].apply(lambda geom: geom.intersection(bbox))
|
||||||
|
|
||||||
logger.info(f"Completed geohash processing. {len(results)} geohash tiles contain water polygons.")
|
if new_coastline.empty:
|
||||||
return results
|
logger.info(f"No coastline found for {geohash_code}. Skipping update.")
|
||||||
|
continue
|
||||||
|
|
||||||
def append_to_existing_parquet(geohash_tiles, output_dir):
|
# Assign "natural: coastline" tag
|
||||||
"""Appends sliced water polygons to existing parquet files after subtracting islands."""
|
new_coastline["tags"] = [{"natural": "coastline"}] * len(new_coastline)
|
||||||
for ghash, gdf in tqdm(geohash_tiles.items(), desc="Appending data"):
|
|
||||||
parquet_path = os.path.join(output_dir, f"{ghash}.parquet")
|
|
||||||
|
|
||||||
# Ensure input water polygons are in EPSG:4326
|
# Ensure columns match before merging
|
||||||
gdf = gdf.to_crs("EPSG:4326")
|
for col in ["feature_id", "tags"]:
|
||||||
|
if col not in filtered_gdf.columns:
|
||||||
|
filtered_gdf[col] = None
|
||||||
|
if col not in new_coastline.columns:
|
||||||
|
new_coastline[col] = None
|
||||||
|
|
||||||
if os.path.exists(parquet_path):
|
# Merge updated coastline into geohash file
|
||||||
try:
|
updated_gdf = gpd.GeoDataFrame(pd.concat([filtered_gdf, new_coastline], ignore_index=True), crs="EPSG:4326")
|
||||||
existing_gdf = gpd.read_parquet(parquet_path)
|
|
||||||
|
|
||||||
# Ensure the existing file is in EPSG:4326 before merging
|
# Save back to parquet
|
||||||
if existing_gdf.crs != "EPSG:4326":
|
updated_gdf.to_parquet(file_path, index=False)
|
||||||
existing_gdf = existing_gdf.to_crs("EPSG:4326")
|
logger.info(f"Updated coastline in {file_path}")
|
||||||
|
|
||||||
# Subtract islands (land features) from the new water polygons
|
|
||||||
gdf = subtract_existing_islands(gdf, existing_gdf)
|
|
||||||
|
|
||||||
# Merge the existing and new data
|
|
||||||
gdf = pd.concat([existing_gdf, gdf], ignore_index=True)
|
|
||||||
logger.debug(f"Appended new water polygons to {ghash}.parquet")
|
|
||||||
except Exception as e:
|
|
||||||
logger.warning(f"Failed to read {parquet_path}. Overwriting instead. Error: {e}")
|
|
||||||
|
|
||||||
# Save the updated file
|
|
||||||
gdf.to_parquet(parquet_path)
|
|
||||||
logger.debug(f"Updated {ghash}.parquet with new data.")
|
|
||||||
|
|
||||||
def save_optional_formats(geohash_tiles, output_dir, export_cbor, export_geojson):
|
|
||||||
"""Saves geohash-sliced water polygons in CBOR and GeoJSON formats."""
|
|
||||||
for ghash, gdf in tqdm(geohash_tiles.items(), desc="Saving optional formats"):
|
|
||||||
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")
|
|
||||||
|
|
||||||
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")
|
|
||||||
|
|
||||||
def main():
|
def main():
|
||||||
parser = argparse.ArgumentParser(description="Slice water polygons along the coastline into geohash tiles.")
|
parser = argparse.ArgumentParser(description="Update coastline geometries in existing geohash parquet files.")
|
||||||
parser.add_argument("coastline", help="Path to the coastline Parquet file.")
|
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("geohash_folder", help="Path to the folder containing <geohash>.parquet files.")
|
||||||
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()
|
args = parser.parse_args()
|
||||||
|
|
||||||
# Set logging level
|
# Load coastline dataset
|
||||||
if args.verbose:
|
logger.info(f"Loading coastline dataset from {args.coastline}...")
|
||||||
logger.setLevel(logging.DEBUG)
|
coastline_gdf = ensure_crs_consistency(gpd.read_parquet(args.coastline))
|
||||||
logger.debug("Verbose logging enabled.")
|
|
||||||
|
|
||||||
logger.info("Loading input datasets...")
|
# Get existing geohash parquet files
|
||||||
coastline = gpd.read_parquet(args.coastline)
|
geohash_files = load_geohash_files(args.geohash_folder)
|
||||||
water_polygons = gpd.read_parquet(args.water)
|
|
||||||
|
|
||||||
# Ensure the geometries are valid
|
if not geohash_files:
|
||||||
coastline = coastline[coastline.geometry.is_valid]
|
logger.warning("No geohash parquet files found. Exiting.")
|
||||||
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.")
|
|
||||||
|
|
||||||
# Get existing geohashes from the output directory
|
|
||||||
existing_geohashes = get_existing_geohashes(args.output_dir)
|
|
||||||
|
|
||||||
if not existing_geohashes:
|
|
||||||
logger.error("No existing geohash parquet files found in the output directory. Exiting.")
|
|
||||||
return
|
return
|
||||||
|
|
||||||
# Intersect water polygons with only the existing geohash tiles
|
# Process geohash files
|
||||||
geohash_tiles = intersect_geohash_tiles(water_polygons, coastline, existing_geohashes)
|
process_geohash_files(coastline_gdf, geohash_files)
|
||||||
|
|
||||||
if not geohash_tiles:
|
logger.info("Processing complete. All coastline updates applied.")
|
||||||
logger.info("No new data to append. Exiting.")
|
|
||||||
return
|
|
||||||
|
|
||||||
# Append results to existing parquet files
|
|
||||||
append_to_existing_parquet(geohash_tiles, args.output_dir)
|
|
||||||
|
|
||||||
# Save optional formats
|
|
||||||
save_optional_formats(geohash_tiles, args.output_dir, args.cbor, args.geojson)
|
|
||||||
|
|
||||||
logger.info("Processing complete. All files updated.")
|
|
||||||
|
|
||||||
if __name__ == "__main__":
|
if __name__ == "__main__":
|
||||||
main()
|
main()
|
||||||
|
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue