crs-l2-svalbard-polar-areas

Rerun

# /// script
# requires-python = ">=3.12"
# dependencies = ["geopandas", "shapely", "pyproj", "pandas"]
# ///
import geopandas as gpd
import pandas as pd

# Load the data
gdf = gpd.read_file("svalbard_glaciers_wgs84.gpkg")

# Check the CRS
print(f"Original CRS: {gdf.crs}")

# Svalbard is roughly at 78 deg N.
# For high latitude area measurements, we want a projection that minimizes distortion.
# A common choice for Svalbard is NSIDC Sea Ice Polar Stereographic (EPSG:3413)
# or UTM zones. However, Polar Stereographic is generally better for high latitude area.
# Let's use EPSG:3413 (NSIDC Sea Ice Polar Stereographic North).
# Another option is EPSG:3031 (WGS 84 / NSIDC Sea Ice Polar Stereographic North).
# Actually, EPSG:3413 is very common for polar regions.
# Let's try EPSG:3413.

target_crs = "EPSG:3413"

# Project the GeoDataFrame
gdf_projected = gdf.to_crs(target_crs)

# Calculate area in square meters and convert to square kilometers
gdf_projected['area_km2'] = gdf_projected.geometry.area / 10**6

# Get top 20 by area
top_20 = gdf_projected.sort_values(by='area_km2', ascending=False).head(20).copy()

# Get bounding box in the projected CRS
# bbox: minx, miny, maxx, maxy
# We need: bbox_minx_polar, bbox_miny_polar, bbox_maxx_polar, bbox_maxy_polar
top_20['bbox_minx_polar'] = top_20.geometry.bounds.minx
top_20['bbox_miny_polar'] = top_20.geometry.bounds.miny
top_20['bbox_maxx_polar'] = top_20.geometry.bounds.maxx
top_20['bbox_maxy_polar'] = top_20.geometry.bounds.maxy

# Extract EPSG code as integer
epsg_code = int(gdf_projected.crs.to_epsg())
top_20['crs_epsg'] = epsg_code

# Prepare the final DataFrame for CSV
# Columns: name, area_km2, bbox_minx_polar, bbox_miny_polar, bbox_maxx_polar, bbox_maxy_polar, crs_epsg
# Assuming 'name' is the column name for glacier name. Let's check columns first.
print(f"Columns in gdf: {gdf.columns.tolist()}")

# Let's re-examine the columns. I will check the column names after printing.
# If 'name' is not the column name, I'll adjust.

output_cols = ['name', 'area_km2', 'bbox_minx_polar', 'bbox_miny_polar', 'bbox_maxx_polar', 'bbox_maxy_polar', 'crs_epsg']
# Wait, I should check if 'name' exists.
if 'name' not in gdf.columns:
    # Check common names like 'NAME', 'glacier_name', etc.
    possible_names = ['NAME', 'name', 'Glacier', 'glacier', 'id']
    for p in possible_names:
        if p in gdf.columns:
            top_20 = top_20.rename(columns={p: 'name'})
            break

# Final selection and export
final_df = top_20[output_cols]
final_df.to_csv("svalbard_glaciers_top20.csv", index=False)

print("Successfully created svalbard_glaciers_top20.csv")