crs-l2-svalbard-polar-areas

Rerun

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

def solve():
    # Load the glacier data
    gdf = gpd.read_file('svalbard_glaciers_wgs84.gpkg')
    
    # Check CRS
    print(f"Original CRS: {gdf.crs}")
    
    # Svalbard is at high latitude. 
    # A good choice for area measurement in Svalbard is a projection centered on Svalbard.
    # UTM Zone 33N (EPSG:32633) or similar? 
    # Actually, for Svalbard, Polar Stereographic is often preferred.
    # EPSG:3413 is NSIDC Sea Ice Polar Stereographic North.
    # Let's try EPSG:3413 or a custom North Polar Stereographic.
    # Actually, EPSG:3413 is a common choice for Arctic regions.
    # Alternatively, EPSG:3031 (ETRS89-extended / LAEA Europe) is good for area, 
    # but for high latitude, Polar Stereographic (EPSG:3413) is more standard for the Arctic.
    
    # Let's check if EPSG:3413 is available and suitable.
    # For Svalbard area, we want something that minimizes distortion.
    # Let's use EPSG:3413.
    target_crs = 'EPSG:3413'
    
    # Project the data
    gdf_projected = gdf.to_crs(target_crs)
    
    # Calculate area in square meters and convert to km2
    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()
    
    # Calculate bounding boxes in projected coordinates
    # bbox returns (minx, miny, maxx, maxy)
    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
    
    # Get EPSG code
    epsg_code = top_20.crs.to_epsg()
    top_20['crs_epsg'] = epsg_code
    
    # Prepare final dataframe with required columns
    # columns: name, area_km2, bbox_minx_polar, bbox_miny_polar, bbox_maxx_polar, bbox_maxy_polar, crs_epsg
    # We need to find what the 'name' column is. Let's check the columns first.
    print(f"Columns in gdf: {gdf.columns.tolist()}")
    
    # Assuming 'name' column exists or might be 'Name' or similar.
    # Let's identify it.
    name_col = None
    for col in ['name', 'Name', 'GLACIER_NAME', 'id']:
        if col in gdf.columns:
            name_col = col
            break
    
    if name_col is None:
        # Fallback if name column is different
        print("Could not find a 'name' column. Using index.")
        top_20['name'] = top_20.index.astype(str)
    else:
        top_20['name'] = top_20[name_col]

    output_df = top_20[['name', 'area_km2', 'bbox_minx_polar', 'bbox_miny_polar', 'bbox_maxx_polar', 'bbox_maxy_polar', 'crs_epsg']]
    
    # Export to CSV
    output_df.to_csv('svalbard_glaciers_top20.csv', index=False)
    print("Successfully created svalbard_glaciers_top20.csv")

if __name__ == "__main__":
    solve()

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

def solve():
    # Load the glacier data
    gdf = gpd.read_file('svalbard_glaciers_wgs84.gpkg')
    
    # Check CRS
    print(f"Original CRS: {gdf.crs}")
    
    # Svalbard is at high latitude. 
    # A good choice for area measurement in Svalbard is a projection centered on Svalbard.
    # EPSG:3413 is NSIDC Sea Ice Polar Stereographic North.
    target_crs = 'EPSG:3413'
    
    # Project the data
    gdf_projected = gdf.to_crs(target_crs)
    
    # Calculate area in square meters and convert to km2
    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()
    
    # Calculate bounding boxes in projected coordinates
    bounds = top_20.geometry.bounds
    top_20['bbox_minx_polar'] = bounds.minx
    top_20['bbox_miny_polar'] = bounds.miny
    top_20['bbox_maxx_polar'] = bounds.maxx
    top_20['bbox_maxy_polar'] = bounds.maxy
    
    # Get EPSG code
    epsg_code = top_20.crs.to_epsg()
    top_20['crs_epsg'] = epsg_code
    
    # Identify 'name' column
    print(f"Columns in gdf: {gdf.columns.tolist()}")
    name_col = None
    for col in ['name', 'Name', 'GLACIER_NAME', 'id', 'NAME']:
        if col in gdf.columns:
            name_col = col
            break
    
    if name_col is None:
        top_20['name'] = top_20.index.astype(str)
    else:
        top_20['name'] = top_20[name_col]

    output_df = top_20[['name', 'area_km2', 'bbox_minx_polar', 'bbox_miny_polar', 'bbox_maxx_polar', 'bbox_maxy_polar', 'crs_epsg']]
    
    # Export to CSV
    output_df.to_csv('svalbard_glaciers_top20.csv', index=False)
    print("Successfully created svalbard_glaciers_top20.csv")

if __name__ == "__main__":
    solve()