geo-l3-antarctica-stations-geodesic

Rerun

duckdb             1.5.3
geopandas          1.1.3
overturemaps       1.0.0
pyogrio            0.12.1
pyproj             3.7.2
shapely            2.1.2

cd /work && python3 -c "
import overturemaps
# Let's see available themes
print(dir(overturemaps))
print('---')
# Check available functions
import inspect
for name in dir(overturemaps):
    if not name.startswith('_'):
        print(name, type(getattr(overturemaps, name)))
"

['__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__path__', '__spec__', 'core', 'geodataframe', 'get_all_overture_types', 'models', 'record_batch_reader']
---
core <class 'module'>
geodataframe <class 'function'>
get_all_overture_types <class 'function'>
models <class 'module'>
record_batch_reader <class 'function'>

cd /work && python3 -c "
import overturemaps
help(overturemaps.record_batch_reader)
"

Help on function record_batch_reader in module overturemaps.core:

record_batch_reader(
    overture_type,
    bbox=None,
    release=None,
    connect_timeout=None,
    request_timeout=None,
    stac=False
) -> Optional[pa.RecordBatchReader]
    Return a pyarrow RecordBatchReader for the desired bounding box and s3 path, or None on error.

cd /work && python3 -c "
import overturemaps
help(overturemaps.geodataframe)
"

Help on function geodataframe in module overturemaps.core:

geodataframe(
    overture_type: str,
    bbox: BBox | tuple[float, float, float, float] | list[float] | None = None,
    release: str = None,
    connect_timeout: int = None,
    request_timeout: int = None,
    stac: bool = False
) -> GeoDataFrame
    Loads geoparquet for specified type into a geopandas dataframe

    Parameters
    ----------
    overture_type: type to load
    bbox: optional bounding box for data fetch (xmin, ymin, xmax, ymax)
    connect_timeout: optional connection timeout in seconds
    request_timeout: optional request timeout in seconds

    Returns
    -------
    GeoDataFrame with the optionally filtered theme data

cd /work && python3 -c "
import overturemaps
types = overturemaps.get_all_overture_types()
print(types)
"

cd /work && python3 -c "
import overturemaps
import geopandas as gpd

# Let's look at what's available in 'place' theme - likely research stations here
# Use a bbox around Antarctica
bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)
print('place columns:', list(places.columns))
print('place rows:', len(places))
print(places[['name', 'subtype', 'country']].head(20))
"

Traceback (most recent call last):
  File "<string>", line 11, in <module>
    print(places[['name', 'subtype', 'country']].head(20))
          ~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.14/site-packages/geopandas/geodataframe.py", line 1891, in __getitem__
    result = super().__getitem__(key)
  File "/usr/local/lib/python3.14/site-packages/pandas/core/frame.py", line 4384, in __getitem__
    indexer = self.columns._get_indexer_strict(key, "columns")[1]
              ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.14/site-packages/pandas/core/indexes/base.py", line 6302, in _get_indexer_strict
    self._raise_if_missing(keyarr, indexer, axis_name)
    ~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.14/site-packages/pandas/core/indexes/base.py", line 6352, in _raise_if_missing
    raise KeyError(f"None of [{key}] are in the [{axis_name}]")
KeyError: "None of [Index(['name', 'subtype', 'country'], dtype='str')] are in the [columns]"

cd /work && python3 -c "
import overturemaps
import geopandas as gpd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)
print('columns:', list(places.columns))
print('rows:', len(places))
# Look at names column
print(places['names'].head(10))
print('---')
print(places['categories'].head(10))
"

columns: ['id', 'geometry', 'categories', 'confidence', 'websites', 'emails', 'socials', 'phones', 'brand', 'addresses', 'names', 'sources', 'operating_status', 'basic_category', 'taxonomy', 'version', 'bbox']
rows: 2348
0    {'primary': 'Tincorp Pa TI', 'common': None, '...
1    {'primary': 'Sara Couto - Intervenção Terapêut...
2    {'primary': 'แหลมใหญ่แล็บ คลินิกเทคนิคการแพทย์...
3    {'primary': 'Chittagong Victory National Schoo...
4    {'primary': 'Quảng Cáo - In Nhanh SaPa', 'comm...
5    {'primary': 'Donata Farms & Land Services', 'c...
6    {'primary': 'Vet Rx Pet Health Center', 'commo...
7    {'primary': 'Work N Wear International', 'comm...
8    {'primary': 'The Wagging Tail Salon', 'common'...
9    {'primary': 'Liso EstÚDio', 'common': None, 'r...
Name: names, dtype: object
---
0    {'primary': 'information_technology_company', ...
1    {'primary': 'physical_therapy', 'alternate': [...
2    {'primary': 'health_and_medical', 'alternate':...
3             {'primary': 'school', 'alternate': None}
4                                                 None
5               {'primary': 'farm', 'alternate': None}
6                                                 None
7        {'primary': 'sports_wear', 'alternate': None}
8        {'primary': 'pet_groomer', 'alternate': None}
9    {'primary': 'professional_services', 'alternat...
Name: categories, dtype: object

cd /work && python3 -c "
import overturemaps
import geopandas as gpd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Look for 'research' or 'station' in names and categories
for idx, row in places.head(50).iterrows():
    cats = str(row['categories'])
    names = str(row['names'])
    if 'station' in cats.lower() or 'research' in cats.lower() or 'station' in names.lower() or 'antarct' in names.lower():
        print(idx, names, cats)
"

cd /work && python3 -c "
import overturemaps
import pandas as pd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Check all unique categories
all_cats = set()
for cats in places['categories']:
    if cats and 'primary' in cats:
        all_cats.add(cats['primary'])
    if cats and 'alternate' in cats and cats['alternate']:
        for c in cats['alternate']:
            all_cats.add(c)
print(sorted(all_cats))
"

Traceback (most recent call last):
  File "<string>", line 13, in <module>
    if cats and 'alternate' in cats and cats['alternate']:
                                        ~~~~^^^^^^^^^^^^^
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

cd /work && python3 -c "
import overturemaps
import pandas as pd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Check all unique categories
all_cats = set()
for cats in places['categories']:
    if cats is not None and isinstance(cats, dict):
        if 'primary' in cats:
            all_cats.add(cats['primary'])
        if 'alternate' in cats and cats['alternate'] is not None:
            for c in cats['alternate']:
                all_cats.add(c)
print(sorted(all_cats))
print('---')
# Also check for names mentioning antarctic or station
for idx, row in places.iterrows():
    names = str(row['names']).lower()
    if 'antarct' in names or ('station' in names and 'research' in names):
        print(idx, row['names'], row['categories'])
"

['accommodation', 'accountant', 'active_life', 'acupuncture', 'adult_entertainment', 'advertising_agency', 'afghan_restaurant', 'agricultural_cooperatives', 'agricultural_service', 'agriculture', 'aircraft_repair', 'airport', 'alcohol_and_drug_treatment_centers', 'amateur_sports_team', 'ambulance_and_ems_services', 'american_restaurant', 'amusement_park', 'animal_shelter', 'antique_store', 'aquatic_pet_store', 'arabian_restaurant', 'arcade', 'archery_range', 'architectural_designer', 'armed_forces_branch', 'art_gallery', 'art_school', 'arts_and_crafts', 'arts_and_entertainment', 'asian_restaurant', 'astrologer', 'attractions_and_activities', 'auction_house', 'audio_visual_equipment_store', 'audiologist', 'auto_company', 'auto_customization', 'auto_detailing', 'auto_manufacturers_and_distributors', 'automation_services', 'automobile_leasing', 'automotive', 'automotive_dealer', 'automotive_parts_and_accessories', 'automotive_repair', 'automotive_services_and_repair', 'b2b_apparel', 'b2b_jewelers', 'b2b_science_and_technology', 'b2b_textiles', 'bakery', 'bank_credit_union', 'banks', 'baptist_church', 'bar', 'bar_and_grill_restaurant', 'barbecue_restaurant', 'barber', 'bartender', 'basketball_stadium', 'beach', 'beauty_and_spa', 'beauty_product_supplier', 'beauty_salon', 'beer_bar', 'bicycle_shop', 'boat_parts_and_accessories', 'boat_rental_and_training', 'boat_service_and_repair', 'boat_tours', 'book_magazine_distribution', 'bookstore', 'bottled_water_company', 'boutique', 'boxing_class', 'brake_service_and_repair', 'brazilian_restaurant', 'breakfast_and_brunch_restaurant', 'bridal_shop', 'brokers', 'bubble_tea', 'buddhist_temple', 'buffet_restaurant', 'building_supply_store', 'burger_restaurant', 'business', 'business_advertising', 'business_manufacturing_and_supply', 'butcher_shop', 'cabin', 'cafe', 'campground', 'campus_building', 'candy_store', 'car_buyer', 'car_dealer', 'car_inspection', 'car_rental_agency', 'car_window_tinting', 'cardiologist', 'caribbean_restaurant', 'carpenter', 'carpet_cleaning', 'caterer', 'catholic_church', 'central_government_office', 'charity_organization', 'cheese_shop', 'chicken_restaurant', 'child_protection_service', 'childrens_clothing_store', 'chinese_restaurant', 'chiropractor', 'chocolatier', 'church_cathedral', 'cinema', 'civilization_museum', 'clothing_company', 'clothing_store', 'cocktail_bar', 'coffee_shop', 'college_university', 'comedy_club', 'commercial_industrial', 'commercial_real_estate', 'commercial_refrigeration', 'commercial_vehicle_dealer', 'community_center', 'community_services_non_profits', 'computer_coaching', 'computer_hardware_company', 'computer_museum', 'computer_store', 'construction_services', 'contractor', 'convenience_store', 'cosmetic_and_beauty_supplies', 'cosmetic_dentist', 'cosmetology_school', 'costume_store', 'cottage', 'counseling_and_mental_health', 'credit_and_debt_counseling', 'cultural_center', 'cupcake_shop', 'custom_clothing', 'dairy_farm', 'dance_club', 'dance_school', 'day_care_preschool', 'delicatessen', 'dentist', 'department_store', 'dermatologist', 'desserts', 'diagnostic_services', 'dim_sum_restaurant', 'diner', 'disability_services_and_support_organization', 'dj_service', 'doctor', 'drive_in_theater', 'driving_school', 'dry_cleaning', 'duty_free_shop', 'e_cigarette_store', 'ear_nose_and_throat', 'eat_and_drink', 'education', 'educational_camp', 'educational_research_institute', 'educational_services', 'educational_supply_store', 'electronics', 'elementary_school', 'embassy', 'emergency_room', 'employment_agencies', 'energy_company', 'energy_equipment_and_solution', 'engine_repair_service', 'engineering_services', 'entertainment_law', 'escape_rooms', 'evangelical_church', 'event_photography', 'event_planning', 'excavation_service', 'exhaust_and_muffler_repair', 'eyewear_and_optician', 'fabric_store', 'fair', 'family_practice', 'farm', 'fashion', 'fashion_accessories_store', 'fast_food_restaurant', 'fertility', 'financial_advising', 'financial_service', 'fish_farm', 'fishmonger', 'fitness_trainer', 'flowers_and_gifts_shop', 'food', 'food_beverage_service_distribution', 'food_consultant', 'food_delivery_service', 'food_stand', 'food_truck', 'fountain', 'freight_and_cargo_service', 'french_restaurant', 'fruits_and_vegetables', 'funeral_services_and_cemeteries', 'furniture_assembly', 'furniture_store', 'garbage_collection_service', 'gas_station', 'general_litigation', 'glass_and_mirror_sales_service', 'glass_manufacturer', 'go_kart_club', 'graphic_designer', 'grocery_store', 'gutter_service', 'gym', 'hair_replacement', 'hair_salon', 'halal_restaurant', 'hardware_store', 'health_and_medical', 'health_food_restaurant', 'health_food_store', 'high_school', 'hindu_temple', 'history_museum', 'hobby_shop', 'hockey_arena', 'holiday_rental_home', 'home_cleaning', 'home_developer', 'home_goods_store', 'home_health_care', 'home_improvement_store', 'home_security', 'home_service', 'home_staging', 'horse_trainer', 'horseback_riding_service', 'hospital', 'hostel', 'hotel', 'hunting_and_fishing_supplies', 'hvac_services', 'hybrid_car_repair', 'ice_cream_shop', 'image_consultant', 'immigration_law', 'indian_restaurant', 'indonesian_restaurant', 'industrial_company', 'industrial_equipment', 'information_technology_company', 'inn', 'installment_loans', 'insurance_agency', 'interior_design', 'internet_marketing_service', 'internet_service_provider', 'investing', 'island', 'it_service_and_computer_repair', 'italian_restaurant', 'japanese_restaurant', 'jewelry_store', 'karaoke', 'korean_restaurant', 'laboratory_testing', 'lake', 'landmark_and_historical_building', 'landscaping', 'language_school', 'laser_eye_surgery_lasik', 'laser_hair_removal', 'laundromat', 'law_enforcement', 'lawyer', 'legal_services', 'library', 'life_coach', 'lighting_store', 'limo_services', 'lingerie_store', 'liquor_store', 'lodge', 'lounge', 'luggage_store', 'machine_shop', 'makeup_artist', 'marina', 'marketing_agency', 'marketing_consultant', 'martial_arts_club', 'mass_media', 'mattress_store', 'meat_wholesaler', 'media_news_company', 'media_restoration_service', 'medical_center', 'medical_research_and_development', 'medical_service_organizations', 'mens_clothing_store', 'merchandising_service', 'metal_fabricator', 'metal_plating_service', 'mexican_restaurant', 'middle_school', 'mobile_home_dealer', 'mobile_phone_store', 'monument', 'mortgage_broker', 'mosque', 'motel', 'motorcycle_dealer', 'motorcycle_manufacturer', 'motorcycle_repair', 'motorsport_vehicle_dealer', 'mountain', 'movers', 'movie_television_studio', 'museum', 'music_production', 'music_school', 'music_venue', 'musical_instrument_store', 'national_park', 'natural_hot_springs', 'naturopathic_holistic', 'non_governmental_association', 'notary_public', 'nursery_and_gardening', 'nutritionist', 'observatory', 'obstetrician_and_gynecologist', 'occupational_therapy', 'office_equipment', 'onsen', 'optometrist', 'oral_surgeon', 'orthodontist', 'orthopedist', 'outlet_store', 'painting', 'panamanian_restaurant', 'pancake_house', 'park', 'party_supply', 'passport_and_visa_services', 'pawn_shop', 'pentecostal_church', 'personal_chef', 'personal_injury_law', 'pet_breeder', 'pet_groomer', 'pet_services', 'pet_store', 'pharmacy', 'photography_store_and_services', 'physical_therapy', 'piercing', 'pilates_studio', 'pizza_restaurant', 'plastic_manufacturer', 'plastic_surgeon', 'playground', 'plumbing', 'political_organization', 'pool_billiards', 'pop_up_shop', 'powder_coating_service', 'preschool', 'print_media', 'printing_services', 'private_association', 'private_investigation', 'private_school', 'professional_services', 'professional_sports_team', 'psychic', 'psychologist', 'psychotherapist', 'pub', 'public_and_government_association', 'public_school', 'public_service_and_government', 'public_utility_company', 'pulmonologist', 'race_track', 'radio_station', 'real_estate', 'real_estate_agent', 'real_estate_investment', 'real_estate_service', 'record_label', 'religious_organization', 'religious_school', 'rental_service', 'resort', 'restaurant', 'restaurant_equipment_and_supply', 'restaurant_wholesale', 'retail', 'retirement_home', 'river', 'roadside_assistance', 'rv_park', 'sandwich_shop', 'school', 'school_sports_team', 'screen_printing_t_shirt_printing', 'scuba_diving_center', 'seafood_restaurant', 'self_storage_facility', 'senior_citizen_services', 'service_apartments', 'sewing_and_alterations', 'shaved_ice_shop', 'shipping_center', 'shoe_repair', 'shoe_store', 'shopping', 'shopping_center', 'sightseeing_tour_agency', 'sign_making', 'skate_park', 'smoothie_juice_bar', 'soccer_field', 'social_media_agency', 'social_media_company', 'social_service_organizations', 'software_development', 'solar_installation', 'souvenir_shop', 'spas', 'specialty_school', 'speech_therapist', 'sporting_goods', 'sports_and_fitness_instruction', 'sports_and_recreation_venue', 'sports_club_and_league', 'sports_wear', 'stadium_arena', 'steakhouse', 'storage_facility', 'structure_and_geography', 'supermarket', 'superstore', 'surfing', 'surgical_appliances_and_supplies', 'surgical_center', 'sushi_restaurant', 'taiwanese_restaurant', 'tanning_salon', 'tapas_bar', 'tattoo_and_piercing', 'tea_room', 'telecommunications_company', 'television_service_providers', 'tennis_court', 'test_preparation', 'texmex_restaurant', 'thai_restaurant', 'theatre', 'theme_restaurant', 'thrift_store', 'tiling', 'tire_dealer_and_repair', 'topic_concert_venue', 'topic_publisher', 'tours', 'town_hall', 'toy_store', 'track_stadium', 'trailer_dealer', 'train_station', 'transportation', 'travel', 'travel_company', 'travel_services', 'trophy_shop', 'truck_rentals', 'trusts', 'turkish_restaurant', 'tutoring_center', 'urban_farm', 'urologist', 'used_car_dealer', 'vegetarian_restaurant', 'venue_and_event_space', 'veterinarian', 'video_game_store', 'videographer', 'vietnamese_restaurant', 'vitamins_and_supplements', 'walk_in_clinic', 'water_heater_installation_repair', 'water_park', 'water_supplier', 'water_treatment_equipment_and_services', 'wedding_planning', 'weight_loss_center', 'wholesale_store', 'wig_store', 'wildlife_sanctuary', 'wills_trusts_and_probate', 'wine_bar', 'winery', 'womens_clothing_store', 'womens_health_clinic', 'wrap_restaurant', 'yoga_studio']
---
11 {'primary': 'Ross Ice Shelf, Antarctica', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
232 {'primary': 'West Antarctica', 'common': None, 'rules': None} {'primary': 'national_park', 'alternate': array(['landmark_and_historical_building'], dtype=object)}
415 {'primary': 'Antarctic Peninsula', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['skate_park', 'structure_and_geography'], dtype=object)}
425 {'primary': 'Mount Wilson (Antarctica)', 'common': None, 'rules': None} {'primary': 'mountain', 'alternate': array(['national_park', 'landmark_and_historical_building'], dtype=object)}
458 {'primary': 'Errera Channel, Antarctica', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain', 'monument'], dtype=object)}
459 {'primary': 'Paradise Bay, Antarctic Peninsula', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['monument', 'structure_and_geography'], dtype=object)}
460 {'primary': 'Brown Station, Antarctic Peninsula', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
462 {'primary': 'González Videla Antarctic Base', 'common': None, 'rules': None} {'primary': 'airport', 'alternate': array(['travel'], dtype=object)}
465 {'primary': 'Port Lockroy, Antarctic Peninsula', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
467 {'primary': 'Mount William (Antarctica)', 'common': None, 'rules': None} {'primary': 'mountain', 'alternate': array(['landmark_and_historical_building'], dtype=object)}
471 {'primary': 'Lemaire Channel, Antarctic Peninsula', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain', 'monument'], dtype=object)}
473 {'primary': 'Palmer Station, Antarctica', 'common': None, 'rules': None} {'primary': 'home_developer', 'alternate': None}
478 {'primary': 'Mount Jackson (Antarctica)', 'common': None, 'rules': None} {'primary': 'mountain', 'alternate': array(['landmark_and_historical_building', 'monument'], dtype=object)}
495 {'primary': 'Esperanto Island, Antarctica', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['island'], dtype=object)}
497 {'primary': 'Cierva Cove, Antarctic Peninsula', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
498 {'primary': 'Antarctica, Antarctic Circle', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['marina', 'fountain'], dtype=object)}
499 {'primary': 'Mikkelsen Harbour, Antarctic Peninsula', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
511 {'primary': 'Bellingshausen Russian Antarctic Station', 'common': None, 'rules': None} {'primary': 'educational_research_institute', 'alternate': array(['public_and_government_association', 'college_university'],
      dtype=object)}
525 {'primary': 'King George Bay (Antarctica)', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['structure_and_geography'], dtype=object)}
566 {'primary': 'Chapel of the Blessed Virgin of Lujan, Antarctica', 'common': None, 'rules': None} {'primary': 'church_cathedral', 'alternate': None}
568 {'primary': 'Esperanza Base, Antarctica', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
584 {'primary': 'British Antarctic Territory', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['public_and_government_association', 'language_school'],
      dtype=object)}
1062 {'primary': 'Antarctica/Troll', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['mountain', 'monument'], dtype=object)}
1562 {'primary': 'Lake Glubokoye (Antarctica)', 'common': None, 'rules': None} {'primary': 'lake', 'alternate': array(['landmark_and_historical_building'], dtype=object)}
1758 {'primary': 'Beaver Lake (Antarctica)', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['lake', 'mountain'], dtype=object)}
1793 {'primary': 'Nicholas Range (Antarctic)', 'common': None, 'rules': None} {'primary': 'mountain', 'alternate': array(['landmark_and_historical_building', 'fountain'], dtype=object)}
1825 {'primary': 'Base antarctique Mawson', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
1900 {'primary': 'East Antarctica', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['restaurant'], dtype=object)}
1904 {'primary': 'Antarctic Kunlun Station', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
2135 {'primary': 'Base antarctique Concordia', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
2183 {'primary': "Base antarctique Dumont-d'Urville", 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['transportation'], dtype=object)}
2236 {'primary': 'Antarctic Plateau', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['structure_and_geography'], dtype=object)}
2289 {'primary': 'Lake Bonney (Antarctica)', 'common': None, 'rules': None} {'primary': 'lake', 'alternate': array(['landmark_and_historical_building', 'monument'], dtype=object)}
2310 {'primary': 'Base antarctique McMurdo', 'common': None, 'rules': None} {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}

cd /work && python3 -c "
import overturemaps
import pandas as pd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Find stations/research bases
research = []
for idx, row in places.iterrows():
    cats = row['categories']
    name = row['names']
    if cats is not None and isinstance(cats, dict):
        primary = cats.get('primary', '')
        alt = cats.get('alternate')
        alt_list = list(alt) if alt is not None else []
        all_cats = [primary] + alt_list
        for c in all_cats:
            if 'research' in str(c) or 'station' in str(c) or 'base' in str(c).lower():
                research.append((idx, name, cats))
                break

print(f'Found {len(research)} potential research stations')
for idx, name, cats in research:
    print(f'{idx}: {name} - {cats}')
"

Found 24 potential research stations
511: {'primary': 'Bellingshausen Russian Antarctic Station', 'common': None, 'rules': None} - {'primary': 'educational_research_institute', 'alternate': array(['public_and_government_association', 'college_university'],
      dtype=object)}
615: {'primary': 'Fuzion Radio Medellin', 'common': None, 'rules': None} - {'primary': 'radio_station', 'alternate': None}
637: {'primary': 'Carrefour Posto João Pessoa', 'common': None, 'rules': array([], dtype=object)} - {'primary': 'gas_station', 'alternate': None}
648: {'primary': 'Carrefour Posto Recife', 'common': None, 'rules': array([], dtype=object)} - {'primary': 'gas_station', 'alternate': None}
1011: {'primary': 'Iceport', 'common': None, 'rules': None} - {'primary': 'landmark_and_historical_building', 'alternate': array(['monument', 'gas_station'], dtype=object)}
1307: {'primary': 'ProfiTech Philippines', 'common': None, 'rules': None} - {'primary': 'investing', 'alternate': array(['educational_research_institute', 'brokers'], dtype=object)}
1337: {'primary': 'Wawa', 'common': None, 'rules': None} - {'primary': 'convenience_store', 'alternate': array(['gas_station', 'coffee_shop', 'fast_food_restaurant',
       'sandwich_shop', 'automotive', 'restaurant', 'retail'],
      dtype=object)}
1340: {'primary': 'Wawa', 'common': None, 'rules': None} - {'primary': 'convenience_store', 'alternate': array(['gas_station', 'coffee_shop', 'automotive', 'restaurant', 'retail'],
      dtype=object)}
1348: {'primary': 'Wawa', 'common': None, 'rules': None} - {'primary': 'convenience_store', 'alternate': array(['gas_station', 'automotive', 'retail'], dtype=object)}
1392: {'primary': 'Wawa', 'common': None, 'rules': None} - {'primary': 'convenience_store', 'alternate': array(['gas_station', 'sandwich_shop', 'automotive', 'restaurant',
       'retail'], dtype=object)}
1394: {'primary': 'Trenova School Support Center', 'common': None, 'rules': None} - {'primary': 'educational_services', 'alternate': array(['educational_research_institute'], dtype=object)}
1483: {'primary': 'Pioneer', 'common': None, 'rules': array([], dtype=object)} - {'primary': 'gas_station', 'alternate': None}
1484: {'primary': 'Pioneer - Gas Station', 'common': None, 'rules': None} - {'primary': 'gas_station', 'alternate': None}
1682: {'primary': 'Wawa', 'common': None, 'rules': None} - {'primary': 'convenience_store', 'alternate': array(['gas_station', 'sandwich_shop', 'automotive', 'restaurant',
       'retail'], dtype=object)}
1825: {'primary': 'Base antarctique Mawson', 'common': None, 'rules': None} - {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
1904: {'primary': 'Antarctic Kunlun Station', 'common': None, 'rules': None} - {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
2048: {'primary': 'Lake Vostok', 'common': None, 'rules': None} - {'primary': 'landmark_and_historical_building', 'alternate': array(['lake', 'train_station'], dtype=object)}
2134: {'primary': 'Base antártica Concordia', 'common': None, 'rules': None} - {'primary': 'train_station', 'alternate': array(['landmark_and_historical_building'], dtype=object)}
2135: {'primary': 'Base antarctique Concordia', 'common': None, 'rules': None} - {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
2146: {'primary': 'Rádio Alerta Penedo', 'common': None, 'rules': None} - {'primary': 'radio_station', 'alternate': None}
2210: {'primary': 'Seaoil Dangcagan', 'common': None, 'rules': None} - {'primary': 'gas_station', 'alternate': None}
2308: {'primary': 'Base McMurdo', 'common': None, 'rules': None} - {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station', 'restaurant'], dtype=object)}
2310: {'primary': 'Base antarctique McMurdo', 'common': None, 'rules': None} - {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
2336: {'primary': 'Lakeview Clinical Research', 'common': None, 'rules': None} - {'primary': 'medical_research_and_development', 'alternate': array(['doctor'], dtype=object)}

cd /work && python3 -c "
import overturemaps
import geopandas as gpd
from shapely import wkt

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Filter to places south of -60 latitude
places_south = places[places.geometry.y < -60].copy()
print(f'Places south of -60: {len(places_south)}')

# Show all unique names
for idx, row in places_south.iterrows():
    lat = row.geometry.y
    lon = row.geometry.x
    name = row['names']
    cats = row['categories']
    cats_str = str(cats.get('primary', '')) if isinstance(cats, dict) else ''
    print(f'{idx}: ({lat:.2f}, {lon:.2f}) - {name} - {cats_str}')
"

Places south of -60: 2348
0: (-84.41, -178.59) - {'primary': 'Tincorp Pa TI', 'common': None, 'rules': None} - information_technology_company
1: (-79.73, -179.97) - {'primary': 'Sara Couto - Intervenção Terapêutica', 'common': None, 'rules': None} - physical_therapy
2: (-78.88, -179.83) - {'primary': 'แหลมใหญ่แล็บ คลินิกเทคนิคการแพทย์ สมุทรสงคราม', 'common': None, 'rules': None} - health_and_medical
3: (-80.42, -175.38) - {'primary': 'Chittagong Victory National School-CVNS', 'common': None, 'rules': None} - school
4: (-84.46, -177.33) - {'primary': 'Quảng Cáo - In Nhanh SaPa', 'common': None, 'rules': None} - 
5: (-82.77, -174.38) - {'primary': 'Donata Farms & Land Services', 'common': None, 'rules': None} - farm
6: (-82.14, -173.46) - {'primary': 'Vet Rx Pet Health Center', 'common': None, 'rules': None} - 
7: (-84.31, -171.70) - {'primary': 'Work N Wear International', 'common': None, 'rules': None} - sports_wear
8: (-82.21, -172.27) - {'primary': 'The Wagging Tail Salon', 'common': None, 'rules': None} - pet_groomer
9: (-77.24, -175.78) - {'primary': 'Liso EstÚDio', 'common': None, 'rules': None} - professional_services
10: (-73.46, -178.14) - {'primary': 'New Launch Property', 'common': None, 'rules': None} - real_estate
11: (-70.61, -177.89) - {'primary': 'Ross Ice Shelf, Antarctica', 'common': None, 'rules': None} - landmark_and_historical_building
12: (-63.47, -176.89) - {'primary': '十刻刺青室', 'common': None, 'rules': None} - 
13: (-64.61, -172.87) - {'primary': 'Poeanuanua', 'common': None, 'rules': None} - retail
14: (-63.47, -171.70) - {'primary': 'Discover Italy', 'common': None, 'rules': None} - tours
15: (-65.12, -172.35) - {'primary': 'สอนต่อขนตาพัทยา ต่อขนตาพัทยา', 'common': None, 'rules': None} - beauty_salon
16: (-66.79, -171.56) - {'primary': 'Krāslavas Varavīksnes vidusskola', 'common': None, 'rules': None} - high_school
17: (-66.37, -167.35) - {'primary': 'مركز هيا هير للشعر الطبيعي/  Haya Center for Natural Hair', 'common': None, 'rules': None} - beauty_salon
18: (-66.23, -168.75) - {'primary': 'Esencia Alquímica: Aromaterapia y Cosmética Natural', 'common': None, 'rules': None} - cosmetic_and_beauty_supplies
19: (-66.37, -169.56) - {'primary': 'Maliha Farzam', 'common': None, 'rules': None} - professional_services
20: (-65.60, -168.26) - {'primary': 'Benjarong Antiques', 'common': None, 'rules': None} - antique_store
21: (-63.73, -170.12) - {'primary': 'The WedPixels', 'common': None, 'rules': None} - event_photography
22: (-65.66, -163.83) - {'primary': 'Clínica Bem Viver', 'common': None, 'rules': None} - medical_service_organizations
23: (-65.69, -163.60) - {'primary': 'Mazuchi Semijoias', 'common': None, 'rules': None} - jewelry_store
24: (-62.96, -164.01) - {'primary': 'لوزان فاشن -louzan fashion', 'common': None, 'rules': None} - boutique
25: (-71.49, -162.95) - {'primary': 'Peixaria do Porto', 'common': None, 'rules': None} - brazilian_restaurant
26: (-76.89, -160.63) - {'primary': 'Maria vilaró fotografía', 'common': None, 'rules': None} - professional_services
27: (-74.45, -165.22) - {'primary': 'Moda Boutique Online', 'common': None, 'rules': None} - boutique
28: (-81.20, -170.16) - {'primary': 'Saaol Heart Center BD Ltd.', 'common': None, 'rules': None} - hospital
29: (-81.72, -163.79) - {'primary': 'ร้านอาหารกลมกล่อมก่อนถึงสนามบิน', 'common': None, 'rules': None} - restaurant
30: (-81.76, -162.87) - {'primary': 'Licera Vision Care Santa Barbara', 'common': None, 'rules': None} - optometrist
31: (-80.71, -155.90) - {'primary': 'Kampala college  Center of English and Training Programs', 'common': None, 'rules': None} - college_university
32: (-78.35, -156.80) - {'primary': 'Muebles Marispal Algeciras', 'common': None, 'rules': None} - furniture_store
33: (-77.53, -156.21) - {'primary': 'Lamoon Cafe & Restaurant - สระบุรี', 'common': None, 'rules': None} - cafe
34: (-74.26, -155.21) - {'primary': 'ราชวิถีสลิมเเคร์ สาขาพิจิตร', 'common': None, 'rules': None} - 
35: (-79.56, -152.58) - {'primary': 'The Future English Language Centre', 'common': None, 'rules': None} - language_school
36: (-83.60, -149.07) - {'primary': 'פיצה גבעת יערים להזמנות חייגו: 0558845021', 'common': None, 'rules': None} - pizza_restaurant
37: (-81.28, -149.68) - {'primary': 'Neringa Geštalto psichoterapeutė', 'common': None, 'rules': None} - psychotherapist
38: (-79.17, -150.47) - {'primary': 'スーパーホテル別府駅前', 'common': None, 'rules': None} - hotel
39: (-81.35, -146.44) - {'primary': 'Exclusive Ladies Undergarments Shop', 'common': None, 'rules': None} - lingerie_store
40: (-82.30, -146.59) - {'primary': '2006hair-Era', 'common': None, 'rules': None} - hair_salon
41: (-82.49, -143.44) - {'primary': 'Cap College Mariveles Bataan', 'common': None, 'rules': None} - college_university
42: (-81.30, -143.09) - {'primary': 'BarGains Lounge', 'common': None, 'rules': None} - lounge
43: (-80.23, -142.82) - {'primary': 'Farm Fresh at UPM', 'common': None, 'rules': None} - italian_restaurant
44: (-80.65, -142.47) - {'primary': 'قرية علي باجلية Ali Bajliyah', 'common': None, 'rules': None} - town_hall
45: (-81.09, -140.98) - {'primary': 'DoubleCare salon and spa', 'common': None, 'rules': None} - beauty_salon
46: (-80.54, -141.69) - {'primary': 'Complete Care Medical Services', 'common': None, 'rules': None} - pharmacy
47: (-79.49, -142.32) - {'primary': 'Cribs and Crayons', 'common': None, 'rules': None} - day_care_preschool
48: (-80.46, -144.59) - {'primary': '白巷子虎尾林森店', 'common': None, 'rules': None} - tea_room
49: (-80.76, -144.84) - {'primary': 'מעדני מניה Maadaney Mania', 'common': None, 'rules': None} - delicatessen
50: (-79.17, -145.20) - {'primary': 'Il Paradiso Delle Calze', 'common': None, 'rules': None} - clothing_store
51: (-77.07, -145.50) - {'primary': 'Ford Ranges', 'common': None, 'rules': None} - mountain
52: (-76.25, -145.01) - {'primary': 'Aire electronica MizaCar', 'common': None, 'rules': None} - automotive
53: (-76.40, -144.19) - {'primary': 'Amanti Restaurant', 'common': None, 'rules': None} - restaurant
54: (-77.98, -143.78) - {'primary': 'HUNG MANH Limited', 'common': None, 'rules': None} - fashion
55: (-73.43, -141.33) - {'primary': 'Simplement St. Onge', 'common': None, 'rules': None} - 
56: (-74.72, -143.79) - {'primary': 'Betty Asesoría inmobiliaria', 'common': None, 'rules': None} - real_estate
57: (-74.80, -144.11) - {'primary': 'Vital Jelly Art Training Class', 'common': None, 'rules': None} - bakery
58: (-78.06, -146.26) - {'primary': 'Harrowgate Hill Methodist Church', 'common': None, 'rules': None} - church_cathedral
59: (-77.76, -148.72) - {'primary': 'พ่วงข้างมอเตอร์ไซต์ พงษ์ศิริการช่าง วงเวียนไดโนเสาร์ภูเวียง', 'common': None, 'rules': None} - motorcycle_repair
60: (-77.33, -149.20) - {'primary': 'Moody Island', 'common': None, 'rules': None} - landmark_and_historical_building
61: (-77.16, -149.77) - {'primary': 'The BuddyCat บ้านสก็อตติชโฟลด์', 'common': None, 'rules': None} - urban_farm
62: (-70.69, -149.21) - {'primary': 'MIZINO', 'common': None, 'rules': None} - womens_clothing_store
63: (-68.03, -148.43) - {'primary': 'Kids Spot Mount Gambier', 'common': None, 'rules': None} - art_school
64: (-65.76, -146.88) - {'primary': 'Te Fare Informatique', 'common': None, 'rules': None} - hardware_store
65: (-70.80, -143.85) - {'primary': 'Xuởng giày khiêu vũ chuyên nghiệp PK - Chuyên sỉ', 'common': None, 'rules': None} - shoe_store
66: (-70.25, -141.99) - {'primary': 'VAST Equipment & Truck Sales', 'common': None, 'rules': None} - car_dealer
67: (-66.43, -140.98) - {'primary': 'Hamburgueria do Du', 'common': None, 'rules': None} - 
68: (-63.47, -142.65) - {'primary': 'أشرف عبداللطيف لقطع غيار المرسيدس ', 'common': None, 'rules': None} - 
69: (-63.55, -142.73) - {'primary': 'شارك SHARK', 'common': None, 'rules': None} - shopping
70: (-60.90, -142.51) - {'primary': 'Kuya Insurance', 'common': None, 'rules': None} - insurance_agency
71: (-65.29, -148.08) - {'primary': 'Magui style', 'common': None, 'rules': None} - retail
72: (-65.29, -149.89) - {'primary': 'Alpha Kennels Australia', 'common': None, 'rules': None} - pet_breeder
73: (-61.10, -149.09) - {'primary': 'Healthdirect', 'common': None, 'rules': None} - school
74: (-67.64, -156.84) - {'primary': 'الفريد للنجف والتحف', 'common': None, 'rules': None} - art_gallery
75: (-71.97, -158.91) - {'primary': 'Jom Naiktaraf Masjid Al Muttaqin', 'common': None, 'rules': None} - 
76: (-66.37, -159.27) - {'primary': 'أفرشة المساجد', 'common': None, 'rules': None} - 
77: (-64.62, -156.86) - {'primary': 'শ্রী রাজ পিনাকী তান্ত্রিক', 'common': None, 'rules': None} - beauty_and_spa
78: (-61.61, -122.05) - {'primary': 'Over The Top Yard Cards', 'common': None, 'rules': None} - event_planning
79: (-61.61, -124.45) - {'primary': 'Sri Ramanujar Engineering College - Vandalur, Chennai - 127.', 'common': None, 'rules': None} - college_university
80: (-62.78, -125.11) - {'primary': "Andie's Furniture", 'common': None, 'rules': None} - furniture_store
81: (-64.47, -123.75) - {'primary': 'Clinica Veterinária Baschirotto', 'common': None, 'rules': None} - pet_store
82: (-64.47, -123.39) - {'primary': 'Sooper Vegee', 'common': None, 'rules': None} - eat_and_drink
83: (-65.77, -121.37) - {'primary': 'Suzuki Nilphamari - Lahin Motors', 'common': None, 'rules': None} - motorcycle_dealer
84: (-69.47, -121.70) - {'primary': 'Goggrll Edappal', 'common': None, 'rules': None} - burger_restaurant
85: (-70.08, -121.39) - {'primary': 'YK Violin Music Class - Mandalay', 'common': None, 'rules': None} - 
86: (-70.08, -123.33) - {'primary': 'andrei.canton', 'common': None, 'rules': None} - 
87: (-66.04, -123.74) - {'primary': 'lichudurlach.arte', 'common': None, 'rules': None} - 
88: (-71.93, -128.10) - {'primary': 'One Moms Battle', 'common': None, 'rules': None} - community_services_non_profits
89: (-68.42, -129.57) - {'primary': 'Espaço Rosileide Domingues - Make Up, Sobrancelhas e Micropigmentação', 'common': None, 'rules': None} - beauty_salon
90: (-65.38, -128.29) - {'primary': 'Easy water - อีซี่วอเตอร์ เครื่องกรองน้ำราคาปลีก-ส่ง', 'common': None, 'rules': None} - 
91: (-64.12, -128.37) - {'primary': 'Jus nature Moucouaré', 'common': None, 'rules': None} - cocktail_bar
92: (-62.59, -126.56) - {'primary': 'Rockview Sunset Beach Resort -Patar Bolinao Pangasinan', 'common': None, 'rules': None} - hotel
93: (-63.47, -128.50) - {'primary': 'Sparklin D’ Lites LLC', 'common': None, 'rules': None} - shaved_ice_shop
94: (-61.27, -127.54) - {'primary': 'Shri Ram Marble Sirsa', 'common': None, 'rules': None} - professional_services
95: (-63.62, -132.45) - {'primary': 'The AI Mortgage Pro', 'common': None, 'rules': None} - real_estate
96: (-65.29, -134.99) - {'primary': 'Feiyu Mosaic', 'common': None, 'rules': None} - 
97: (-60.07, -134.16) - {'primary': 'barberia_melli', 'common': None, 'rules': None} - barber
98: (-60.85, -138.16) - {'primary': 'Blastic_007', 'common': None, 'rules': None} - 
99: (-64.71, -139.14) - {'primary': 'Cuadrodeco', 'common': None, 'rules': None} - 
100: (-65.29, -136.35) - {'primary': 'Super-Gains', 'common': None, 'rules': None} - vitamins_and_supplements
101: (-69.17, -139.53) - {'primary': 'ICAN Management Services', 'common': None, 'rules': None} - professional_services
102: (-69.16, -139.92) - {'primary': '老康泰日間照顧中心', 'common': None, 'rules': None} - community_services_non_profits
103: (-71.07, -137.11) - {'primary': 'Los Vecinos', 'common': None, 'rules': None} - 
104: (-70.08, -137.77) - {'primary': 'Moon Art', 'common': None, 'rules': None} - 
105: (-69.66, -136.04) - {'primary': 'Royal Bedsheet & Towels', 'common': None, 'rules': None} - furniture_store
106: (-69.81, -136.23) - {'primary': 'Karpicentro', 'common': None, 'rules': None} - professional_services
107: (-71.52, -134.30) - {'primary': 'Melanoma Kills', 'common': None, 'rules': None} - community_services_non_profits
108: (-69.41, -135.00) - {'primary': 'Farias Business Communications', 'common': None, 'rules': None} - professional_services
109: (-66.37, -134.58) - {'primary': 'Beshie Pasabuy', 'common': None, 'rules': None} - shopping
110: (-74.83, -132.00) - {'primary': 'Hobbs Coast', 'common': None, 'rules': None} - beach
111: (-75.59, -134.49) - {'primary': 'Nội Thất Thuý Lưu -Thiết Kế Thi Công Trọn Gói -0971.032.051', 'common': None, 'rules': None} - professional_services
112: (-74.97, -133.60) - {'primary': 'PP Al Muslim Beganjing', 'common': None, 'rules': None} - religious_organization
113: (-72.94, -134.88) - {'primary': 'Temple Martial Arts Kickboxing - Northwich', 'common': None, 'rules': None} - martial_arts_club
114: (-74.78, -136.41) - {'primary': "Samira's Hair Fashion +599 98693311", 'common': None, 'rules': None} - beauty_salon
115: (-75.32, -137.11) - {'primary': 'PT. Karya Taka Indonesia', 'common': None, 'rules': None} - real_estate
116: (-72.30, -138.92) - {'primary': 'KN Squared Constructions and Supplies', 'common': None, 'rules': None} - real_estate
117: (-75.43, -137.70) - {'primary': 'CAMBO Cleaner', 'common': None, 'rules': None} - 
118: (-76.77, -136.47) - {'primary': 'Resurrection Church', 'common': None, 'rules': None} - church_cathedral
119: (-81.08, -136.55) - {'primary': 'Iriga City Public Library', 'common': None, 'rules': None} - public_and_government_association
120: (-79.43, -137.81) - {'primary': 'Desia Koraput', 'common': None, 'rules': None} - tours
121: (-81.07, -139.11) - {'primary': 'Delight Sofa Care', 'common': None, 'rules': None} - home_cleaning
122: (-84.54, -137.11) - {'primary': 'Dinas Kebudayaan dan Pariwisata Provinsi Jambi', 'common': None, 'rules': None} - public_and_government_association
123: (-81.92, -137.81) - {'primary': 'Dinastia Nortena', 'common': None, 'rules': None} - 
124: (-84.15, -135.42) - {'primary': 'La casa del deporte', 'common': None, 'rules': None} - clothing_store
125: (-81.82, -135.70) - {'primary': 'World Vision Church of God', 'common': None, 'rules': None} - church_cathedral
126: (-80.76, -135.00) - {'primary': 'Kasa Sul Materiais de Construção', 'common': None, 'rules': None} - building_supply_store
127: (-78.77, -133.59) - {'primary': 'Rumah Makan Sari Raos Bu Endang', 'common': None, 'rules': None} - indonesian_restaurant
128: (-79.46, -131.95) - {'primary': 'ข้าวลงเข่งบางแสน', 'common': None, 'rules': None} - thai_restaurant
129: (-79.55, -131.92) - {'primary': 'Colegio Villa La Granja Mulchen', 'common': None, 'rules': None} - 
130: (-79.17, -132.19) - {'primary': 'Von Star Power', 'common': None, 'rules': None} - pet_store
131: (-80.30, -132.19) - {'primary': 'Amelis Center', 'common': None, 'rules': None} - beauty_salon
132: (-82.12, -132.89) - {'primary': 'Diskop UKM Jateng', 'common': None, 'rules': None} - public_and_government_association
133: (-84.84, -127.49) - {'primary': 'Ser Esencia Restaurant', 'common': None, 'rules': None} - restaurant
134: (-83.12, -127.79) - {'primary': 'Bildøy Marina', 'common': None, 'rules': None} - 
135: (-81.05, -126.30) - {'primary': 'On the docks', 'common': None, 'rul
[... truncated ...]mary': 'Seaoil Dangcagan', 'common': None, 'rules': None} - gas_station
2211: (-84.34, 140.48) - {'primary': 'โนนสำราญบ้านเฮา', 'common': None, 'rules': None} - community_services_non_profits
2212: (-81.36, 138.87) - {'primary': '汽車收容所', 'common': None, 'rules': None} - car_dealer
2213: (-76.03, 138.59) - {'primary': 'Dra Gracy Alencar Aesthetics', 'common': None, 'rules': None} - 
2214: (-73.83, 140.08) - {'primary': 'Fashion in DHAKA', 'common': None, 'rules': None} - womens_clothing_store
2215: (-74.42, 140.73) - {'primary': 'Moorunga Angus', 'common': None, 'rules': None} - real_estate
2216: (-77.74, 140.28) - {'primary': 'UPT SDN Kedawung 01', 'common': None, 'rules': None} - elementary_school
2217: (-78.62, 142.91) - {'primary': 'Phụ Tùng Kim Khí Giá Rẻ', 'common': None, 'rules': None} - 
2218: (-78.21, 142.73) - {'primary': 'Bremer Bay Mini Golf & Kiosk', 'common': None, 'rules': None} - tours
2219: (-79.69, 144.14) - {'primary': 'The Seasoning Outlet', 'common': None, 'rules': None} - shopping
2220: (-84.40, 142.38) - {'primary': '絵画教室\u3000みんなのアトリエ', 'common': None, 'rules': None} - art_school
2221: (-82.27, 143.55) - {'primary': 'บัส สบาย เจิด', 'common': None, 'rules': None} - cafe
2222: (-67.39, 144.90) - {'primary': 'ヒッコリー', 'common': None, 'rules': None} - 
2223: (-75.48, 145.21) - {'primary': 'KinnKoff coffee', 'common': None, 'rules': None} - restaurant
2224: (-79.31, 146.25) - {'primary': 'Knox Pharmacy', 'common': None, 'rules': None} - pharmacy
2225: (-62.92, 146.60) - {'primary': 'The Maris - Căn Hộ Biển Vệ Tinh', 'common': None, 'rules': None} - real_estate
2226: (-63.23, 145.90) - {'primary': 'Маҷаллаи “Дин ва ҷомеа”', 'common': None, 'rules': None} - media_news_company
2227: (-63.83, 146.80) - {'primary': 'Кучковская Татьяна', 'common': None, 'rules': None} - beauty_salon
2228: (-70.38, 146.95) - {'primary': 'Beauty Salon Tý Hair', 'common': None, 'rules': None} - beauty_salon
2229: (-84.76, 148.85) - {'primary': 'Green Pc', 'common': None, 'rules': None} - computer_store
2230: (-83.25, 149.76) - {'primary': 'Reems Closet Deals', 'common': None, 'rules': None} - thrift_store
2231: (-81.52, 149.95) - {'primary': 'กองการศึกษาองค์การบริหารส่วนตำบลกันทรารมย์', 'common': None, 'rules': None} - public_and_government_association
2232: (-81.99, 148.59) - {'primary': 'Od Eliana Vaca', 'common': None, 'rules': None} - 
2233: (-81.10, 148.04) - {'primary': 'Los taquitos de pm - Recursos Humanos', 'common': None, 'rules': None} - restaurant
2234: (-74.95, 147.60) - {'primary': 'Starr Pharmacy Roban', 'common': None, 'rules': None} - health_and_medical
2235: (-74.44, 149.23) - {'primary': 'MAX Reprise', 'common': None, 'rules': None} - automotive_repair
2236: (-77.00, 150.00) - {'primary': 'Antarctic Plateau', 'common': None, 'rules': None} - landmark_and_historical_building
2237: (-74.40, 150.47) - {'primary': 'Chácara Espaço Colorê Festas e Eventos', 'common': None, 'rules': None} - event_planning
2238: (-80.53, 151.49) - {'primary': 'ボディメイクサロンJolly', 'common': None, 'rules': None} - health_and_medical
2239: (-80.25, 151.88) - {'primary': 'Nadhi/நதி - The Online Saree Store', 'common': None, 'rules': None} - 
2240: (-82.31, 151.17) - {'primary': 'Pak Polytechnic Institute', 'common': None, 'rules': None} - college_university
2241: (-83.71, 153.26) - {'primary': '海港青年商會 JCI Harbour', 'common': None, 'rules': None} - community_services_non_profits
2242: (-83.89, 153.64) - {'primary': 'Asha Physiotherapy Center', 'common': None, 'rules': None} - physical_therapy
2243: (-83.73, 154.68) - {'primary': 'Shree Aakash Electronics', 'common': None, 'rules': None} - electronics
2244: (-82.89, 156.62) - {'primary': 'Eglise Chretienne Evangelique au Congo/kinshasa', 'common': None, 'rules': None} - religious_organization
2245: (-79.87, 156.23) - {'primary': 'Mount Ellis', 'common': None, 'rules': None} - landmark_and_historical_building
2246: (-74.46, 154.41) - {'primary': 'Warning Five จำหน่ายสายนาฬิกาหนังแท้ สายนาฬิกา smartwatch', 'common': None, 'rules': None} - 
2247: (-71.75, 153.98) - {'primary': 'Estratto［エストラット］', 'common': None, 'rules': None} - retail
2248: (-71.88, 153.56) - {'primary': 'Majesty Café - Nabas, Aklan', 'common': None, 'rules': None} - coffee_shop
2249: (-71.61, 155.48) - {'primary': "Ocean View Apt- Cox'sbazar", 'common': None, 'rules': None} - hotel
2250: (-72.87, 156.37) - {'primary': 'Bit7awa', 'common': None, 'rules': None} - internet_marketing_service
2251: (-69.43, 155.73) - {'primary': 'মুক্তি হারবাল চিকিৎসালয়-Mukti Herbal Chikitsalay', 'common': None, 'rules': None} - family_practice
2252: (-63.55, 156.88) - {'primary': 'معمل بركات مسلم ابن عقيل', 'common': None, 'rules': None} - 
2253: (-63.15, 154.09) - {'primary': 'APEC HOMES Bulacan-Erica Baltazar', 'common': None, 'rules': None} - real_estate_agent
2254: (-64.33, 153.81) - {'primary': 'Built by Kidz Parties', 'common': None, 'rules': None} - school
2255: (-66.76, 153.36) - {'primary': 'Face Bar Cle', 'common': None, 'rules': None} - 
2256: (-63.23, 153.64) - {'primary': 'Big Smoke Bbq', 'common': None, 'rules': None} - auto_detailing
2257: (-61.94, 153.44) - {'primary': 'Live Life Loud Designs', 'common': None, 'rules': None} - 
2258: (-61.38, 152.81) - {'primary': 'Nông Sản Mai Cồ', 'common': None, 'rules': None} - farm
2259: (-63.74, 152.37) - {'primary': 'Mitsubishi Pili Camarines Sur - By Stephen Camacho', 'common': None, 'rules': None} - 
2260: (-68.91, 150.47) - {'primary': 'Herrlich Beauty By Cloriss', 'common': None, 'rules': None} - beauty_salon
2261: (-70.08, 149.72) - {'primary': 'RangRays Media Works', 'common': None, 'rules': None} - mass_media
2262: (-67.07, 147.66) - {'primary': 'Little Honey NZ', 'common': None, 'rules': None} - 
2263: (-60.49, 148.19) - {'primary': 'Academic Assist Co.', 'common': None, 'rules': None} - boutique
2264: (-60.21, 150.47) - {'primary': 'E17 Telecom', 'common': None, 'rules': None} - internet_service_provider
2265: (-60.41, 153.28) - {'primary': 'Peluquería Gloria Frutillar', 'common': None, 'rules': None} - beauty_and_spa
2266: (-60.24, 158.55) - {'primary': 'Tramites al Volante', 'common': None, 'rules': None} - driving_school
2267: (-64.17, 168.40) - {'primary': 'Más Tiempo Juntos MX', 'common': None, 'rules': None} - health_and_medical
2268: (-62.45, 166.44) - {'primary': 'The Navyandhra School - Jr Wing', 'common': None, 'rules': None} - school
2269: (-65.29, 167.18) - {'primary': 'Two Swords Tattoo Co.', 'common': None, 'rules': None} - tattoo_and_piercing
2270: (-72.33, 168.14) - {'primary': 'Highrise Scenic Farm Resort', 'common': None, 'rules': None} - hotel
2271: (-70.67, 166.09) - {'primary': 'กุ้งไทย  "นสพ.กุ้งไทย"', 'common': None, 'rules': None} - mass_media
2272: (-70.36, 165.89) - {'primary': '𝗬𝗢𝗨𝗦𝗦𝗘𝗙 𝗕𝗘𝗡𝗧𝗔𝗟𝗕𝗘', 'common': None, 'rules': None} - 
2273: (-70.08, 164.13) - {'primary': 'North Store', 'common': None, 'rules': None} - shopping
2274: (-63.98, 165.28) - {'primary': 'ธนพรพลาสแพ็ค ถุงพลาสติก-สกรีนโลโก้ ราคาโรงงาน', 'common': None, 'rules': None} - industrial_company
2275: (-63.47, 164.11) - {'primary': 'Fitiavana Mickael', 'common': None, 'rules': None} - 
2276: (-64.77, 163.12) - {'primary': 'Redeemed Christian Church Of God Kingdom Connection Chapel in Goliad Texas', 'common': None, 'rules': None} - church_cathedral
2277: (-66.79, 162.42) - {'primary': "D'Luxxe Manufacturing Corp.", 'common': None, 'rules': None} - beauty_product_supplier
2278: (-66.37, 161.87) - {'primary': 'Mercado Grazi', 'common': None, 'rules': None} - grocery_store
2279: (-62.35, 160.78) - {'primary': 'Bewas dil muhingi', 'common': None, 'rules': None} - 
2280: (-70.05, 160.22) - {'primary': 'مركز سنان الجراحي ـ Snan Surgical Center', 'common': None, 'rules': None} - emergency_room
2281: (-69.41, 160.12) - {'primary': 'Electric Help BD', 'common': None, 'rules': None} - 
2282: (-69.66, 158.91) - {'primary': 'Bengkel las bojonegoro', 'common': None, 'rules': None} - professional_services
2283: (-71.85, 160.33) - {'primary': 'วัดนอก ตำบลไชยคราม  อำเภอดอนสัก จังหวัดสุราษฎร์ธานี', 'common': None, 'rules': None} - buddhist_temple
2284: (-72.72, 160.37) - {'primary': 'Mount Koons', 'common': None, 'rules': None} - landmark_and_historical_building
2285: (-68.91, 161.72) - {'primary': 'Checkers Custard', 'common': None, 'rules': None} - eat_and_drink
2286: (-68.35, 162.07) - {'primary': 'OWWO Boneless Korean Fried Chicken', 'common': None, 'rules': None} - chicken_restaurant
2287: (-76.55, 163.26) - {'primary': 'Chatime Philippines', 'common': None, 'rules': None} - 
2288: (-77.62, 162.98) - {'primary': 'Canada Glacier', 'common': None, 'rules': None} - landmark_and_historical_building
2289: (-77.72, 162.42) - {'primary': 'Lake Bonney (Antarctica)', 'common': None, 'rules': None} - lake
2290: (-77.73, 162.17) - {'primary': 'Geleira de Taylor', 'common': None, 'rules': None} - monument
2291: (-77.72, 161.63) - {'primary': 'Lake Joyce', 'common': None, 'rules': None} - lake
2292: (-77.55, 160.77) - {'primary': 'Jackson Pond', 'common': None, 'rules': None} - lake
2293: (-73.96, 161.31) - {'primary': 'LTS Aircon Services', 'common': None, 'rules': None} - engineering_services
2294: (-73.20, 159.16) - {'primary': 'HONDA Mạnh Hùng Quỳnh Côi', 'common': None, 'rules': None} - retail
2295: (-79.43, 162.42) - {'primary': 'Mr.FAWZ Restaurant', 'common': None, 'rules': None} - hotel
2296: (-84.52, 164.75) - {'primary': 'CALIRA & CHESKA  Villa & Resort', 'common': None, 'rules': None} - holiday_rental_home
2297: (-77.88, 164.41) - {'primary': 'Psicóloga Serly C. de Souza', 'common': None, 'rules': None} - psychotherapist
2298: (-74.69, 164.11) - {'primary': 'Mario-Zucchelli-Station', 'common': None, 'rules': None} - 
2299: (-78.63, 164.88) - {'primary': 'Trường Mầm Non An Khánh A', 'common': None, 'rules': None} - school
2300: (-77.96, 166.52) - {'primary': 'Pegasus Field', 'common': None, 'rules': None} - airport
2301: (-77.89, 166.59) - {'primary': 'Cape Armitage', 'common': None, 'rules': None} - landmark_and_historical_building
2302: (-77.85, 166.67) - {'primary': 'McMurdo General Hospital', 'common': None, 'rules': None} - hospital
2303: (-77.85, 166.67) - {'primary': 'ZiroZiro', 'common': None, 'rules': None} - religious_organization
2304: (-77.85, 166.67) - {'primary': 'Southern Exposure', 'common': None, 'rules': None} - bar
2305: (-77.85, 166.67) - {'primary': 'Gallagher Pub', 'common': None, 'rules': None} - bar
2306: (-77.85, 166.66) - {'primary': 'Coffee House', 'common': None, 'rules': None} - coffee_shop
2307: (-77.85, 166.66) - {'primary': 'Chapel of the Snows', 'common': None, 'rules': None} - church_cathedral
2308: (-77.85, 166.67) - {'primary': 'Base McMurdo', 'common': None, 'rules': None} - landmark_and_historical_building
2309: (-77.85, 166.67) - {'primary': 'Polární stanice McMurdo', 'common': None, 'rules': None} - landmark_and_historical_building
2310: (-77.85, 166.67) - {'primary': 'Base antarctique McMurdo', 'common': None, 'rules': None} - landmark_and_historical_building
2311: (-77.53, 167.28) - {'primary': 'Monte Erebus', 'common': None, 'rules': None} - landmark_and_historical_building
2312: (-77.33, 167.52) - {'primary': 'Lewis Bay', 'common': None, 'rules': None} - beach
2313: (-83.95, 166.16) - {'primary': 'Steve Kilbey', 'common': None, 'rules': None} - 
2314: (-83.83, 168.05) - {'primary': 'Risen Savior Academy Davao City', 'common': None, 'rules': None} - school
2315: (-84.97, 168.40) - {'primary': 'One World Immigration Karachi', 'common': None, 'rules': None} - immigration_law
2316: (-83.48, 170.30) - {'primary': 'King Glacier', 'common': None, 'rules': None} - landmark_and_historical_building
2317: (-72.94, 170.57) - {'primary': 'Kimpich Store', 'common': None, 'rules': None} - beauty_and_spa
2318: (-73.10, 169.05) - {'primary': 'Narrow Neck', 'common': None, 'rules': None} - landmark_and_historical_building
2319: (-61.88, 169.80) - {'primary': 'Sales Store', 'common': None, 'rules': None} - 
2320: (-60.49, 170.83) - {'primary': 'Cương Audio', 'common': None, 'rules': None} - 
2321: (-63.77, 172.63) - {'primary': 'DreamHome Property', 'common': None, 'rules': None} - real_estate
2322: (-67.59, 170.91) - {'primary': 'Megan Smith - Lodge Real Estate REAA 2008', 'common': None, 'rules': None} - real_estate_agent
2323: (-71.75, 170.87) - {'primary': 'Vaříme s Mony', 'common': None, 'rules': None} - 
2324: (-79.81, 171.91) - {'primary': 'กล่องสุ่ม อาหารทะเล499 กุ้งมังกร ปลีก-ส่ง By . ออม seafood', 'common': None, 'rules': None} - seafood_restaurant
2325: (-74.59, 172.27) - {'primary': 'Kitchora', 'common': None, 'rules': None} - real_estate
2326: (-74.45, 172.45) - {'primary': 'Çebu Gayrımenkul', 'common': None, 'rules': None} - real_estate
2327: (-61.63, 173.79) - {'primary': 'GNL Marketing', 'common': None, 'rules': None} - marketing_agency
2328: (-63.47, 174.37) - {'primary': 'Ace Volleyball & Athletics', 'common': None, 'rules': None} - sports_club_and_league
2329: (-63.86, 174.14) - {'primary': 'PR by Judinil Berrios', 'common': None, 'rules': None} - 
2330: (-78.24, 174.78) - {'primary': 'مؤسسة العدالة القانونيه', 'common': None, 'rules': None} - lawyer
2331: (-82.30, 176.03) - {'primary': '悦聖翡翠', 'common': None, 'rules': None} - fashion
2332: (-81.69, 176.48) - {'primary': 'បក្សីប្រយុទ្ធ- Dk666', 'common': None, 'rules': None} - 
2333: (-79.95, 177.09) - {'primary': 'Van Anh Cleaning', 'common': None, 'rules': None} - home_cleaning
2334: (-63.36, 177.48) - {'primary': 'FREAK of RAVE', 'common': None, 'rules': None} - event_planning
2335: (-61.77, 175.91) - {'primary': 'Brits School of English', 'common': None, 'rules': None} - language_school
2336: (-64.41, 175.78) - {'primary': 'Lakeview Clinical Research', 'common': None, 'rules': None} - medical_research_and_development
2337: (-61.53, 179.68) - {'primary': 'Latidos De Fe', 'common': None, 'rules': None} - social_service_organizations
2338: (-62.03, 179.97) - {'primary': 'La Cura Tv', 'common': None, 'rules': None} - 
2339: (-62.93, 180.00) - {'primary': 'Jania Rocha Advocacia', 'common': None, 'rules': None} - lawyer
2340: (-64.12, 179.90) - {'primary': 'Confeitaria Berries', 'common': None, 'rules': None} - 
2341: (-63.26, 179.44) - {'primary': 'Thu huyen Mobile', 'common': None, 'rules': None} - mobile_phone_store
2342: (-65.64, 179.66) - {'primary': 'EnterTrainment Junction Groups and Events', 'common': None, 'rules': None} - museum
2343: (-63.49, 179.30) - {'primary': 'printcraft', 'common': None, 'rules': None} - sign_making
2344: (-62.84, 179.20) - {'primary': 'Hotel 3 banderas Manzanillo de Mar', 'common': None, 'rules': None} - hotel
2345: (-63.17, 177.55) - {'primary': 'Sonrieme Clinica Dental', 'common': None, 'rules': None} - 
2346: (-84.44, 178.47) - {'primary': 'ไก่หมุนจิระภา สาขาสามแยกวัดป่าคา เฉลิมพระเกียรติ สระบุรี By ครกทอง', 'common': None, 'rules': None} - thai_restaurant
2347: (-79.18, 179.92) - {'primary': 'Monitos Bebechitos', 'common': None, 'rules': None} - childrens_clothing_store

cd /work && python3 -c "
import overturemaps
import geopandas as gpd

# Let's check if the bbox filtering works
bbox = [-60, -90, 60, -60]  # Just a small part of Antarctica
places = overturemaps.geodataframe('place', bbox=bbox)
print(f'Places with bbox around 0 meridian: {len(places)}')
if len(places) > 0:
    print(places[['id', 'names', 'categories']].head(20))
print('---')

# Also check with no bbox
places_all = overturemaps.geodataframe('place')
print(f'All places: {len(places_all)}')
"

cd /work && timeout 120 python3 -c "
import overturemaps
import geopandas as gpd

# First, get all places (no bbox) and check coordinate ranges
# This might be slow but let's try
places = overturemaps.geodataframe('place')
print(f'Total places: {len(places)}')
print('CRS:', places.crs)
print('geom types:', places.geometry.type.value_counts())
bounds = places.geometry.bounds
print('x range:', bounds.minx.min(), '-', bounds.maxx.max())
print('y range:', bounds.miny.min(), '-', bounds.maxy.max())
" 2>&1 | head -30

cd /work && timeout 120 python3 -c "
import overturemaps
import geopandas as gpd
import json

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Look for Antarctic stations by checking names
ant_stations = []
for idx, row in places.iterrows():
    name_dict = row['names']
    primary_name = ''
    if name_dict and isinstance(name_dict, dict):
        primary_name = str(name_dict.get('primary', ''))
    
    # Check if name suggests Antarctic station
    name_lower = primary_name.lower()
    keywords = ['station', 'base', 'antarct', 'research', 'polar']
    if any(k in name_lower for k in keywords):
        lat = row.geometry.y
        lon = row.geometry.x
        cats = row['categories']
        ant_stations.append({
            'idx': idx,
            'name': primary_name,
            'lat': lat,
            'lon': lon,
            'categories': str(cats) if cats else ''
        })

print(f'Found {len(ant_stations)} potential stations')
for s in sorted(ant_stations, key=lambda x: x['name']):
    print(f'{s[\"idx\"]}: ({s[\"lat\"]:.2f}, {s[\"lon\"]:.2f}) - {s[\"name\"]}')
    print(f'   cats: {s[\"categories\"]}')
"

Found 51 potential stations
1904: (-80.42, 77.12) - Antarctic Kunlun Station
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
415: (-69.50, -65.00) - Antarctic Peninsula
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['skate_park', 'structure_and_geography'], dtype=object)}
2236: (-77.00, 150.00) - Antarctic Plateau
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['structure_and_geography'], dtype=object)}
498: (-64.05, -60.93) - Antarctica, Antarctic Circle
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['marina', 'fountain'], dtype=object)}
1062: (-72.01, 2.54) - Antarctica/Troll
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['mountain', 'monument'], dtype=object)}
565: (-64.24, -56.62) - Base Antártica Marambio
   cats: {'primary': 'public_and_government_association', 'alternate': array(['central_government_office', 'accommodation'], dtype=object)}
569: (-63.40, -57.00) - Base Esperanza, Antartida Argentina
   cats: {'primary': 'airport', 'alternate': array(['travel', 'accommodation'], dtype=object)}
563: (-64.24, -56.62) - Base Marambio, Antartida Argentina
   cats: {'primary': 'central_government_office', 'alternate': array(['accommodation'], dtype=object)}
2308: (-77.85, 166.67) - Base McMurdo
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station', 'restaurant'], dtype=object)}
2135: (-75.10, 123.56) - Base antarctique Concordia
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
2183: (-66.66, 140.00) - Base antarctique Dumont-d'Urville
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['transportation'], dtype=object)}
1825: (-67.60, 62.87) - Base antarctique Mawson
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
2310: (-77.85, 166.67) - Base antarctique McMurdo
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
2134: (-75.10, 123.33) - Base antártica Concordia
   cats: {'primary': 'train_station', 'alternate': array(['landmark_and_historical_building'], dtype=object)}
1758: (-70.78, 68.30) - Beaver Lake (Antarctica)
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['lake', 'mountain'], dtype=object)}
511: (-62.20, -58.96) - Bellingshausen Russian Antarctic Station
   cats: {'primary': 'educational_research_institute', 'alternate': array(['public_and_government_association', 'college_university'],
      dtype=object)}
584: (-75.00, -50.00) - British Antarctic Territory
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['public_and_government_association', 'language_school'],
      dtype=object)}
461: (-64.89, -62.87) - Brown Station
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
460: (-64.89, -62.87) - Brown Station, Antarctic Peninsula
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
2101: (-84.61, 116.72) - COFFE Station Old town
   cats: {'primary': 'cafe', 'alternate': array(['coffee_shop'], dtype=object)}
504: (-62.24, -58.67) - Carlini Base
   cats: {'primary': 'arts_and_entertainment', 'alternate': None}
566: (-64.24, -56.63) - Chapel of the Blessed Virgin of Lujan, Antarctica
   cats: {'primary': 'church_cathedral', 'alternate': None}
497: (-64.15, -60.91) - Cierva Cove, Antarctic Peninsula
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
1900: (-80.00, 80.00) - East Antarctica
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['restaurant'], dtype=object)}
458: (-64.70, -62.59) - Errera Channel, Antarctica
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain', 'monument'], dtype=object)}
495: (-62.64, -60.36) - Esperanto Island, Antarctica
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['island'], dtype=object)}
568: (-63.40, -57.00) - Esperanza Base, Antarctica
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
1211: (-65.34, 15.47) - Fuel Tech Micro Gas Station & Dispensers
   cats: {'primary': 'industrial_equipment', 'alternate': None}
462: (-64.82, -62.86) - González Videla Antarctic Base
   cats: {'primary': 'airport', 'alternate': array(['travel'], dtype=object)}
525: (-62.10, -58.08) - King George Bay (Antarctica)
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['structure_and_geography'], dtype=object)}
2289: (-77.72, 162.42) - Lake Bonney (Antarctica)
   cats: {'primary': 'lake', 'alternate': array(['landmark_and_historical_building', 'monument'], dtype=object)}
1562: (-67.67, 45.87) - Lake Glubokoye (Antarctica)
   cats: {'primary': 'lake', 'alternate': array(['landmark_and_historical_building'], dtype=object)}
2336: (-64.41, 175.78) - Lakeview Clinical Research
   cats: {'primary': 'medical_research_and_development', 'alternate': array(['doctor'], dtype=object)}
471: (-65.10, -63.98) - Lemaire Channel, Antarctic Peninsula
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain', 'monument'], dtype=object)}
564: (-64.24, -56.63) - Marambio Base
   cats: {'primary': 'airport', 'alternate': array(['armed_forces_branch', 'arts_and_entertainment'], dtype=object)}
2298: (-74.69, 164.11) - Mario-Zucchelli-Station
   cats: 
439: (-75.14, -66.80) - Master Bees Stationeries & Souvenirs
   cats: {'primary': 'souvenir_shop', 'alternate': array(['bookstore'], dtype=object)}
499: (-63.90, -60.77) - Mikkelsen Harbour, Antarctic Peninsula
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
478: (-71.36, -63.41) - Mount Jackson (Antarctica)
   cats: {'primary': 'mountain', 'alternate': array(['landmark_and_historical_building', 'monument'], dtype=object)}
467: (-64.78, -63.68) - Mount William (Antarctica)
   cats: {'primary': 'mountain', 'alternate': array(['landmark_and_historical_building'], dtype=object)}
425: (-68.45, -65.55) - Mount Wilson (Antarctica)
   cats: {'primary': 'mountain', 'alternate': array(['national_park', 'landmark_and_historical_building'], dtype=object)}
1010: (-70.65, -8.25) - Neumayer-Station II
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['restaurant'], dtype=object)}
490: (-80.87, -61.35) - New Lucena Municipal Police Station
   cats: {'primary': 'public_service_and_government', 'alternate': array(['community_center'], dtype=object)}
1793: (-66.67, 55.62) - Nicholas Range (Antarctic)
   cats: {'primary': 'mountain', 'alternate': array(['landmark_and_historical_building', 'fountain'], dtype=object)}
473: (-64.77, -64.05) - Palmer Station, Antarctica
   cats: {'primary': 'home_developer', 'alternate': None}
459: (-64.91, -62.90) - Paradise Bay, Antarctic Peninsula
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['monument', 'structure_and_geography'], dtype=object)}
1484: (-78.91, 43.93) - Pioneer - Gas Station
   cats: {'primary': 'gas_station', 'alternate': None}
1418: (-79.25, 40.50) - Plateau Station
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['armed_forces_branch', 'travel'], dtype=object)}
465: (-64.82, -63.49) - Port Lockroy, Antarctic Peninsula
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
11: (-70.61, -177.89) - Ross Ice Shelf, Antarctica
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
232: (-79.00, -100.00) - West Antarctica
   cats: {'primary': 'national_park', 'alternate': array(['landmark_and_historical_building'], dtype=object)}

cd /work && timeout 120 python3 -c "
import overturemaps
import geopandas as gpd

bbox = [-180, -90, 180, -60]

# Let's look at the land theme
land = overturemaps.geodataframe('land', bbox=bbox)
print('land columns:', list(land.columns))
print('land rows:', len(land))
print('land subtypes:', land['subtype'].value_counts() if 'subtype' in land.columns else 'N/A')
print(land.head())
"

land columns: ['id', 'names', 'subtype', 'class', 'sources', 'source_tags', 'level', 'wikidata', 'surface', 'elevation', 'geometry', 'version', 'bbox']
land rows: 93367
land subtypes: subtype
rock        56249
land        21921
physical    14268
glacier       842
reef           36
sand           32
wetland         9
crater          7
grass           2
shrub           1
Name: count, dtype: int64
                                     id  ...                                               bbox
0  23f89603-d14d-3f1b-97c5-314269e17921  ...  {'xmin': -180.0, 'xmax': -179.97145080566406, ...
1  ef0ba39b-922b-3790-b71c-b5238e8f4e9b  ...  {'xmin': -179.9757080078125, 'xmax': -179.9248...
2  46a1ad5c-4af5-32f8-b0f5-c416ec492555  ...  {'xmin': -179.95298767089844, 'xmax': -179.859...
3  c2519229-0a95-3761-a671-5a44ea2e2194  ...  {'xmin': -179.76327514648438, 'xmax': -179.697...
4  8848f703-84f7-3c4c-9b0c-8a93394dceac  ...  {'xmin': -179.80276489257812, 'xmax': -179.603...

[5 rows x 13 columns]

cd /work && timeout 120 python3 -c "
import overturemaps
import geopandas as gpd

bbox = [-180, -90, 180, -60]

# Get water
water = overturemaps.geodataframe('water', bbox=bbox)
print('water columns:', list(water.columns))
print('water rows:', len(water))
if 'subtype' in water.columns:
    print('water subtypes:', water['subtype'].value_counts())
# Check geometry types
print('geom types:', water.geometry.type.value_counts().head(10))
"

water columns: ['id', 'names', 'subtype', 'class', 'sources', 'source_tags', 'level', 'wikidata', 'is_intermittent', 'is_salt', 'geometry', 'version', 'bbox']
water rows: 13501
water subtypes: subtype
stream        5374
ocean         4978
lake          1746
physical       809
pond           372
water          180
river           34
canal            5
human_made       2
reservoir        1
Name: count, dtype: int64
geom types: Polygon         7445
LineString      5421
Point            625
MultiPolygon      10
Name: count, dtype: int64

cd /work && timeout 120 python3 -c "
import overturemaps
import geopandas as gpd

# Get bathymetry
bathy = overturemaps.geodataframe('bathymetry')
print('bathymetry columns:', list(bathy.columns))
print('bathymetry rows:', len(bathy))
if 'subtype' in bathy.columns:
    print('subtypes:', bathy['subtype'].value_counts())
print('geom types:', bathy.geometry.type.value_counts())
"

bathymetry columns: ['id', 'depth', 'cartography', 'geometry', 'sources', 'version', 'bbox']
bathymetry rows: 59963
geom types: MultiPolygon    59870
Polygon            93
Name: count, dtype: int64

cd /work && timeout 120 python3 << 'EOF'
import overturemaps
import geopandas as gpd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Identify research stations more carefully
# Known Antarctic stations by name pattern
station_names = [
    'Bellingshausen', 'McMurdo', 'Mawson', 'Concordia', 'Dumont', 'Marambio',
    'Esperanza', 'Carlini', 'Palmer', 'Brown', 'González Videla', 'Neumayer',
    'Plateau', 'Kunlun', 'Mario-Zucchelli', 'Zhongshan', 'Progress', 'Halley',
    'Casey', 'Scott', 'Amundsen-Scott', 'South Pole', 'Vostok', 'Mirny',
    'Novolazarevskaya', 'Syowa', 'Maitri', 'Bharati', 'Troll', 'SANAE',
    'Davis', 'Macquarie Island', 'Molodezhnaya', 'Bellinghausen',
]

# Also look for 'station' or 'base' in name
stations = []
for idx, row in places.iterrows():
    name_dict = row['names']
    primary_name = ''
    if name_dict and isinstance(name_dict, dict):
        primary_name = str(name_dict.get('primary', ''))
    
    lat = row.geometry.y
    lon = row.geometry.x
    
    # Must be south of -60
    if lat > -60:
        continue
    
    name_lower = primary_name.lower()
    
    # Check if it's a known station or has station/base in name
    is_station = False
    
    # Match known station names
    for sname in station_names:
        if sname.lower() in name_lower:
            is_station = True
            break
    
    # Or contains 'station' or 'base' in name, but filter out gas_station, radio_station, police_station, etc.
    if ('station' in name_lower or 'base' in name_lower):
        cats = row['categories']
        if cats and isinstance(cats, dict):
            cat_primary = cats.get('primary', '')
            cat_str = str(cats)
            # Exclude gas_station, radio_station, police_station
            if 'gas_station' not in cat_str and 'radio_station' not in cat_str and 'police' not in cat_str.lower():
                is_station = True
    
    if is_station:
        stations.append({
            'idx': idx,
            'name': primary_name,
            'lat': lat,
            'lon': lon,
            'categories': str(row['categories']) if row['categories'] else ''
        })

print(f'Found {len(stations)} stations')
for s in sorted(stations, key=lambda x: x['name']):
    print(f'({s["lat"]:.4f}, {s["lon"]:.4f}) - {s["name"]}')
    print(f'   cats: {s["categories"]}')
EOF

Found 35 stations
(-80.4169, 77.1161) - Antarctic Kunlun Station
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
(-77.0000, 150.0000) - Antarctic Plateau
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['structure_and_geography'], dtype=object)}
(-72.0114, 2.5350) - Antarctica/Troll
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['mountain', 'monument'], dtype=object)}
(-64.2374, -56.6235) - Base Antártica Marambio
   cats: {'primary': 'public_and_government_association', 'alternate': array(['central_government_office', 'accommodation'], dtype=object)}
(-63.3980, -56.9974) - Base Esperanza, Antartida Argentina
   cats: {'primary': 'airport', 'alternate': array(['travel', 'accommodation'], dtype=object)}
(-64.2413, -56.6198) - Base Marambio, Antartida Argentina
   cats: {'primary': 'central_government_office', 'alternate': array(['accommodation'], dtype=object)}
(-77.8500, 166.6667) - Base McMurdo
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station', 'restaurant'], dtype=object)}
(-75.1000, 123.5580) - Base antarctique Concordia
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
(-66.6628, 140.0013) - Base antarctique Dumont-d'Urville
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['transportation'], dtype=object)}
(-67.6026, 62.8736) - Base antarctique Mawson
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
(-77.8500, 166.6667) - Base antarctique McMurdo
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['train_station'], dtype=object)}
(-75.1000, 123.3333) - Base antártica Concordia
   cats: {'primary': 'train_station', 'alternate': array(['landmark_and_historical_building'], dtype=object)}
(-62.1998, -58.9624) - Bellingshausen Russian Antarctic Station
   cats: {'primary': 'educational_research_institute', 'alternate': array(['public_and_government_association', 'college_university'],
      dtype=object)}
(-64.8950, -62.8700) - Brown Station
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
(-64.8921, -62.8724) - Brown Station, Antarctic Peninsula
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
(-84.6078, 116.7188) - COFFE Station Old town
   cats: {'primary': 'cafe', 'alternate': array(['coffee_shop'], dtype=object)}
(-62.2379, -58.6668) - Carlini Base
   cats: {'primary': 'arts_and_entertainment', 'alternate': None}
(-74.3470, 62.5314) - EEC Paroisse Plateau Pointe-Noire1
   cats: {'primary': 'evangelical_church', 'alternate': None}
(-63.3977, -56.9961) - Esperanza Base, Antarctica
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['fountain'], dtype=object)}
(-65.3393, 15.4688) - Fuel Tech Micro Gas Station & Dispensers
   cats: {'primary': 'industrial_equipment', 'alternate': None}
(-64.8234, -62.8575) - González Videla Antarctic Base
   cats: {'primary': 'airport', 'alternate': array(['travel'], dtype=object)}
(-77.5000, 106.0000) - Lake Vostok
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['lake', 'train_station'], dtype=object)}
(-64.2410, -56.6268) - Marambio Base
   cats: {'primary': 'airport', 'alternate': array(['armed_forces_branch', 'arts_and_entertainment'], dtype=object)}
(-74.6943, 164.1125) - Mario-Zucchelli-Station
   cats: 
(-75.1408, -66.7969) - Master Bees Stationeries & Souvenirs
   cats: {'primary': 'souvenir_shop', 'alternate': array(['bookstore'], dtype=object)}
(-77.8468, 166.6693) - McMurdo General Hospital
   cats: {'primary': 'hospital', 'alternate': array(['medical_center'], dtype=object)}
(-70.6500, -8.2500) - Neumayer-Station II
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['restaurant'], dtype=object)}
(-80.8717, -61.3477) - New Lucena Municipal Police Station
   cats: {'primary': 'public_service_and_government', 'alternate': array(['community_center'], dtype=object)}
(-64.7743, -64.0531) - Palmer Station, Antarctica
   cats: {'primary': 'home_developer', 'alternate': None}
(-79.2500, 40.5000) - Plateau Station
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['armed_forces_branch', 'travel'], dtype=object)}
(-77.8500, 166.6667) - Polární stanice McMurdo
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['restaurant'], dtype=object)}
(-72.6667, 2.8500) - SANAE
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['monument'], dtype=object)}
(-67.1333, 51.3500) - Scott Mountains
   cats: {'primary': 'mountain', 'alternate': array(['landmark_and_historical_building', 'monument'], dtype=object)}
(-71.9572, 2.4667) - Troll Airfield
   cats: {'primary': 'airport', 'alternate': array(['travel'], dtype=object)}
(-78.4644, 106.8373) - Vostok İstasyonu
   cats: {'primary': 'landmark_and_historical_building', 'alternate': array(['restaurant'], dtype=object)}

cd /work && timeout 180 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Comprehensive list of known Antarctic research stations
# We'll also look for stations by explicit name patterns
known_station_names = [
    'McMurdo', 'Amundsen-Scott', 'South Pole', 'Mawson', 'Davis', 'Casey',
    'Scott Base', 'Scott', 'Halley', 'Syowa', 'Dumont', 'Dumont d\'Urville',
    'Mirny', 'Novolazarevskaya', 'Bellingshausen', 'Vostok', 'Progress',
    'Zhongshan', 'Great Wall', 'Changcheng', 'Kunlun', 'Taishan', 'Qinling',
    'Troll', 'SANAE', 'Maitri', 'Bharati', 'Concordia', 'Palmer',
    'Marambio', 'Esperanza', 'Carlini', 'Jubany', 'Orcadas', 'San Martín',
    'Belgrano', 'Brown', 'Matienzo', 'Primavera', 'Decepción', 'Cámara',
    'Melo', 'Petrel', 'Neumayer', 'Georg von Neumayer', 'Kohnen',
    'Wasa', 'Svea', 'Aboa', 'Princess Elisabeth', 'Rothera',
    'Bird Island', 'Signy', 'King Edward Point', 'Macquarie',
    'Mario Zucchelli', 'Zucchelli', 'Mario-Zucchelli',
    'Jang Bogo', 'King Sejong', 'Gabriel de Castilla', 'Juan Carlos I',
    'Antonio de los Ríos', 'Ohridski', 'St. Kliment Ohridski',
    'Arctowski', 'Henryk Arctowski', 'Vernadsky', 'Faraday',
    'Aboa', 'Plateau Station', 'Plateau',
    'Frei', 'Eduardo Frei', 'Frei Montalva', 'Teniente Marsh',
    'Rothera Point', 'Law-Racoviță', 'Law Racovita',
    'Comandante Ferraz', 'Ferraz',
    'General Bernardo O\'Higgins', 'O\'Higgins',
    'Presidente Eduardo Frei', 'Frei',
    'González Videla', 'Videla',
    'Artigas', 'Arturo Parodi', 'Teniente Arturo Parodi',
    'Esperanza', 'Carlini', 'Decepción',
    'Halcón', 'Carvajal', 'Maldonado',
    'Professor Julio Escudero',
    'Yelcho', 'Guillermo Mann',
    'Risopatrón', 'Risopatron',
    'Luis Risopatrón',
    'Villa Las Estrellas',
    'Plateau Station',
    'Lake Vostok',
    'Neumayer-Station',
]

# Filter places south of -60
places_south = places[places.geometry.y < -60].copy()

stations = []
seen_coords = set()

for idx, row in places_south.iterrows():
    name_dict = row['names']
    primary_name = ''
    if name_dict and isinstance(name_dict, dict):
        primary_name = str(name_dict.get('primary', ''))
    
    lat = round(row.geometry.y, 4)
    lon = round(row.geometry.x, 4)
    name_lower = primary_name.lower()
    
    is_station = False
    
    # Check known station names
    for sname in known_station_names:
        if sname.lower() in name_lower:
            is_station = True
            break
    
    # Also check for 'station' or 'base' keywords
    if not is_station and ('station' in name_lower or 'base' in name_lower):
        cats = row['categories']
        if cats and isinstance(cats, dict):
            cat_str = str(cats).lower()
            # Exclude non-research stations
            excl = ['gas_station', 'radio_station', 'police', 'coffee', 'old town']
            if not any(e in cat_str for e in excl):
                is_station = True
    
    if is_station:
        coord_key = (lat, lon, primary_name.lower().strip())
        if coord_key not in seen_coords:
            seen_coords.add(coord_key)
            stations.append({
                'station_id': f'station_{len(stations):03d}',
                'station_name': primary_name,
                'latitude': row.geometry.y,
                'longitude': row.geometry.x,
                'geometry': row.geometry
            })

stations_gdf = gpd.GeoDataFrame(stations, crs='EPSG:4326')
print(f'Found {len(stations_gdf)} unique stations')
stations_gdf.to_parquet('/work/raw_stations.parquet')
print(stations_gdf[['station_id', 'station_name', 'latitude', 'longitude']].to_string())
PYEOF

Found 36 unique stations
     station_id                                     station_name   latitude   longitude
0   station_000                              Melodic Music Point -81.798228 -111.844740
1   station_001             Master Bees Stationeries & Souvenirs -75.140778  -66.796875
2   station_002               Brown Station, Antarctic Peninsula -64.892130  -62.872353
3   station_003                                    Brown Station -64.895000  -62.870000
4   station_004                   González Videla Antarctic Base -64.823441  -62.857514
5   station_005                       Palmer Station, Antarctica -64.774345  -64.053125
6   station_006              New Lucena Municipal Police Station -80.871738  -61.347656
7   station_007                                     Carlini Base -62.237900  -58.666842
8   station_008         Bellingshausen Russian Antarctic Station -62.199803  -58.962403
9   station_009               Base Marambio, Antartida Argentina -64.241263  -56.619756
10  station_010                                    Marambio Base -64.241014  -56.626753
11  station_011                          Base Antártica Marambio -64.237433  -56.623535
12  station_012                       Esperanza Base, Antarctica -63.397672  -56.996144
13  station_013              Base Esperanza, Antartida Argentina -63.397988  -56.997430
14  station_014                              Neumayer-Station II -70.650000   -8.250000
15  station_015                                   Troll Airfield -71.957222    2.466667
16  station_016                                 Antarctica/Troll -72.011389    2.535000
17  station_017                                            SANAE -72.666667    2.850000
18  station_018         Fuel Tech Micro Gas Station & Dispensers -65.339347   15.468750
19  station_019  Windfreight Express Total Logistics Inc.(WETLI) -75.073010   43.066406
20  station_020                                  Plateau Station -79.250000   40.500000
21  station_021                                  Scott Mountains -67.133333   51.350000
22  station_022               EEC Paroisse Plateau Pointe-Noire1 -74.346979   62.531433
23  station_023                          Base antarctique Mawson -67.602580   62.873600
24  station_024                         Antarctic Kunlun Station -80.416944   77.116111
25  station_025                                      Lake Vostok -77.500000  106.000000
26  station_026                                 Vostok İstasyonu -78.464422  106.837328
27  station_027                         Base antártica Concordia -75.100000  123.333333
28  station_028                       Base antarctique Concordia -75.100000  123.558000
29  station_029                Base antarctique Dumont-d'Urville -66.662833  140.001333
30  station_030                                Antarctic Plateau -77.000000  150.000000
31  station_031                          Mario-Zucchelli-Station -74.694254  164.112482
32  station_032                         McMurdo General Hospital -77.846764  166.669301
33  station_033                                     Base McMurdo -77.850000  166.666667
34  station_034                          Polární stanice McMurdo -77.850000  166.666667
35  station_035                         Base antarctique McMurdo -77.850000  166.666667

cd /work && timeout 180 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Manual curated list of known Antarctic research stations with their approximate coordinates
# Based on known stations from COMNAP / official lists
known_stations = {
    # Name variants that appear in Overture
    'McMurdo': (-77.85, 166.67),
    'Bellingshausen Russian Antarctic Station': (-62.20, -58.96),
    'Bellingshausen': (-62.20, -58.96),
    'Base antarctique McMurdo': (-77.85, 166.67),
    'Base antarctique Dumont-d\'Urville': (-66.66, 140.00),
    'Base antarctique Mawson': (-67.60, 62.87),
    'Base antarctique Concordia': (-75.10, 123.56),
    'Base antártica Concordia': (-75.10, 123.33),
    'Base McMurdo': (-77.85, 166.67),
    'Polární stanice McMurdo': (-77.85, 166.67),
    'Mario-Zucchelli-Station': (-74.69, 164.11),
    'Neumayer-Station II': (-70.65, -8.25),
    'Plateau Station': (-79.25, 40.50),
    'Antarctic Kunlun Station': (-80.42, 77.12),
    'Vostok İstasyonu': (-78.46, 106.84),
    'Lake Vostok': (-77.50, 106.00),
    'Palmer Station, Antarctica': (-64.77, -64.05),
    'Brown Station, Antarctic Peninsula': (-64.89, -62.87),
    'Brown Station': (-64.90, -62.87),
    'González Videla Antarctic Base': (-64.82, -62.86),
    'Carlini Base': (-62.24, -58.67),
    'Esperanza Base, Antarctica': (-63.40, -57.00),
    'Base Esperanza, Antartida Argentina': (-63.40, -57.00),
    'Marambio Base': (-64.24, -56.63),
    'Base Marambio, Antartida Argentina': (-64.24, -56.62),
    'Base Antártica Marambio': (-64.24, -56.62),
    'SANAE': (-72.67, 2.85),
    'Troll Airfield': (-71.96, 2.47),
    'Antarctica/Troll': (-72.01, 2.54),
}

# Also look for more stations by checking place names that match known stations
# Let's search all places south of -60 for any that might be research stations
places_south = places[places.geometry.y < -60].copy()

# Build a more comprehensive list
station_records = []
seen_stations = set()  # (lat_rounded, lon_rounded, name_key)

for idx, row in places_south.iterrows():
    name_dict = row['names']
    primary_name = ''
    if name_dict and isinstance(name_dict, dict):
        primary_name = str(name_dict.get('primary', ''))
    
    lat = row.geometry.y
    lon = row.geometry.x
    name_lower = primary_name.lower().strip()
    
    # Check known station list
    matched = False
    for known_name, (klat, klon) in known_stations.items():
        if known_name.lower() in name_lower:
            matched = True
            break
    
    if not matched:
        continue
    
    # Deduplicate by rounding coordinates
    lat_r = round(lat, 1)
    lon_r = round(lon, 1)
    key = (lat_r, lon_r)
    
    if key not in seen_stations:
        seen_stations.add(key)
        station_records.append({
            'station_id': f'st{len(station_records):03d}',
            'station_name': primary_name,
            'latitude': lat,
            'longitude': lon,
            'geometry': row.geometry
        })

stations_gdf = gpd.GeoDataFrame(station_records, crs='EPSG:4326')
print(f'Found {len(stations_gdf)} unique stations')
print(stations_gdf[['station_id', 'station_name', 'latitude', 'longitude']].to_string())

# Save to parquet
stations_gdf.to_parquet('/work/raw_stations.parquet')
PYEOF

Found 20 unique stations
   station_id                              station_name   latitude   longitude
0       st000        Brown Station, Antarctic Peninsula -64.892130  -62.872353
1       st001            González Videla Antarctic Base -64.823441  -62.857514
2       st002                Palmer Station, Antarctica -64.774345  -64.053125
3       st003                              Carlini Base -62.237900  -58.666842
4       st004  Bellingshausen Russian Antarctic Station -62.199803  -58.962403
5       st005        Base Marambio, Antartida Argentina -64.241263  -56.619756
6       st006                Esperanza Base, Antarctica -63.397672  -56.996144
7       st007                       Neumayer-Station II -70.650000   -8.250000
8       st008                            Troll Airfield -71.957222    2.466667
9       st009                                     SANAE -72.666667    2.850000
10      st010                           Plateau Station -79.250000   40.500000
11      st011                   Base antarctique Mawson -67.602580   62.873600
12      st012                  Antarctic Kunlun Station -80.416944   77.116111
13      st013                               Lake Vostok -77.500000  106.000000
14      st014                          Vostok İstasyonu -78.464422  106.837328
15      st015                  Base antártica Concordia -75.100000  123.333333
16      st016                Base antarctique Concordia -75.100000  123.558000
17      st017         Base antarctique Dumont-d'Urville -66.662833  140.001333
18      st018                   Mario-Zucchelli-Station -74.694254  164.112482
19      st019                  McMurdo General Hospital -77.846764  166.669301

cd /work && timeout 180 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd

bbox = [-180, -90, 180, -60]

# Check division_area for Antarctica
div_area = overturemaps.geodataframe('division_area', bbox=bbox)
print('division_area columns:', list(div_area.columns))
print('rows:', len(div_area))
# Check for Antarctica-related
for idx, row in div_area.head(50).iterrows():
    names = row.get('names', {})
    if names and isinstance(names, dict):
        name = names.get('primary', '')
        if 'antarct' in str(name).lower() or 'south pole' in str(name).lower():
            print(f'{idx}: {name} - subtype: {row.get("subtype", "")}')
PYEOF

division_area columns: ['id', 'geometry', 'country', 'sources', 'subtype', 'admin_level', 'class', 'names', 'is_land', 'is_territorial', 'region', 'division_id', 'version', 'bbox']
rows: 3
1: Antarctica - subtype: country
2: Antarctica - subtype: country

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# More carefully identify research stations
# Key categories that indicate research stations
research_categories = ['educational_research_institute', 'research_station']
# Also look for names that are clearly Antarctic research stations

# All places south of -60
places_s = places[places.geometry.y < -60].copy()

# Let's look at all unique name patterns with 'station' or 'base' 
# that are actually research stations in Antarctica
candidates = []
for idx, row in places_s.iterrows():
    name_dict = row['names']
    primary_name = ''
    if name_dict and isinstance(name_dict, dict):
        primary_name = str(name_dict.get('primary', ''))
    
    cats = row['categories']
    cat_primary = ''
    cat_alt = []
    if cats and isinstance(cats, dict):
        cat_primary = cats.get('primary', '')
        if cats.get('alternate') is not None:
            cat_alt = list(cats['alternate'])
    
    all_cats = [cat_primary] + cat_alt
    
    is_research = False
    for c in all_cats:
        if 'research' in c.lower():
            is_research = True
            break
    
    # Check name for research station keywords
    name_lower = primary_name.lower()
    
    # Known Antarctic stations by name
    if any(k in name_lower for k in ['antarct', 'station', 'base']) and not any(
        k in name_lower for k in ['gas station', 'radio_station', 'police', 'coffee', 'old town', 'souvenir']
    ):
        if is_research or any(name_part in name_lower for name_part in [
            'bellingshausen', 'mcmurdo', 'mawson', 'concordia', 'dumont',
            'marambio', 'esperanza', 'carlini', 'palmer', 'brown',
            'gonzález', 'videla', 'neumayer', 'plateau', 'kunlun',
            'zucchelli', 'vostok', 'sanae', 'troll', 'mirny',
            'novolazarevskaya', 'syowa', 'maitri', 'bharati',
            'halley', 'casey', 'davis', 'scott', 'rothara',
            'ferraz', 'ohiggins', 'frei', 'marsh', 'king sejong',
            'jang bogo', 'princess elisabeth', 'kohnen', 'aboa',
            'signy', 'bird island', 'vernadsky', 'arctowski',
            'ohridski', 'gabriel de castilla', 'juan carlos',
            'artigas', 'escudero', 'yelcho', 'maldonado',
            'carvajal', 'risopatrón', 'risopatron',
        ]):
            candidates.append({
                'name': primary_name,
                'lat': row.geometry.y,
                'lon': row.geometry.x,
                'cats': cat_primary,
                'alt_cats': cat_alt
            })

print(f'Candidates: {len(candidates)}')
for c in sorted(candidates, key=lambda x: x['name']):
    print(f"({c['lat']:.4f}, {c['lon']:.4f}) - {c['name']} [{c['cats']}]")

PYEOF

Candidates: 23
(-80.4169, 77.1161) - Antarctic Kunlun Station [landmark_and_historical_building]
(-77.0000, 150.0000) - Antarctic Plateau [landmark_and_historical_building]
(-72.0114, 2.5350) - Antarctica/Troll [landmark_and_historical_building]
(-64.2374, -56.6235) - Base Antártica Marambio [public_and_government_association]
(-63.3980, -56.9974) - Base Esperanza, Antartida Argentina [airport]
(-64.2413, -56.6198) - Base Marambio, Antartida Argentina [central_government_office]
(-77.8500, 166.6667) - Base McMurdo [landmark_and_historical_building]
(-75.1000, 123.5580) - Base antarctique Concordia [landmark_and_historical_building]
(-66.6628, 140.0013) - Base antarctique Dumont-d'Urville [landmark_and_historical_building]
(-67.6026, 62.8736) - Base antarctique Mawson [landmark_and_historical_building]
(-77.8500, 166.6667) - Base antarctique McMurdo [landmark_and_historical_building]
(-75.1000, 123.3333) - Base antártica Concordia [train_station]
(-62.1998, -58.9624) - Bellingshausen Russian Antarctic Station [educational_research_institute]
(-64.8950, -62.8700) - Brown Station [landmark_and_historical_building]
(-64.8921, -62.8724) - Brown Station, Antarctic Peninsula [landmark_and_historical_building]
(-62.2379, -58.6668) - Carlini Base [arts_and_entertainment]
(-63.3977, -56.9961) - Esperanza Base, Antarctica [landmark_and_historical_building]
(-64.8234, -62.8575) - González Videla Antarctic Base [airport]
(-64.2410, -56.6268) - Marambio Base [airport]
(-74.6943, 164.1125) - Mario-Zucchelli-Station []
(-70.6500, -8.2500) - Neumayer-Station II [landmark_and_historical_building]
(-64.7743, -64.0531) - Palmer Station, Antarctica [home_developer]
(-79.2500, 40.5000) - Plateau Station [landmark_and_historical_building]

cd /work && timeout 180 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Final curated list: map name variants to canonical station entries
# We want one point per physical station, using the best-named entry
# Coordinates from known COMNAP data:
canonical_stations = {
    # (station_name, lat, lon)
    'McMurdo': ('McMurdo Station', -77.85, 166.67),
    'Bellingshausen': ('Bellingshausen Station', -62.20, -58.96),
    'Palmer': ('Palmer Station', -64.77, -64.05),
    'Mawson': ('Mawson Station', -67.60, 62.87),
    'Davis': ('Davis Station', -68.58, 77.97),
    'Casey': ('Casey Station', -66.28, 110.53),
    'Dumont': ("Dumont d'Urville Station", -66.66, 140.00),
    'Concordia': ('Concordia Station', -75.10, 123.33),
    'Marambio': ('Marambio Base', -64.24, -56.62),
    'Esperanza': ('Esperanza Base', -63.40, -57.00),
    'Carlini': ('Carlini Base', -62.24, -58.67),
    'Brown': ('Brown Station', -64.90, -62.87),
    'Gonzalez Videla': ('González Videla Station', -64.82, -62.86),
    'Neumayer': ('Neumayer Station III', -70.65, -8.25),
    'SANAE': ('SANAE IV', -72.67, 2.85),
    'Troll': ('Troll Station', -72.01, 2.54),
    'Plateau': ('Plateau Station', -79.25, 40.50),
    'Kunlun': ('Kunlun Station', -80.42, 77.12),
    'Vostok': ('Vostok Station', -78.46, 106.84),
    'Mario Zucchelli': ('Mario Zucchelli Station', -74.69, 164.11),
    'Amundsen-Scott': ('Amundsen-Scott South Pole Station', -90.00, 0.00),
    'Halley': ('Halley Research Station', -75.58, -26.50),
    'Rothera': ('Rothera Research Station', -67.57, -68.13),
    'Syowa': ('Syowa Station', -69.00, 39.58),
    'Mirny': ('Mirny Station', -66.56, 93.02),
    'Novolazarevskaya': ('Novolazarevskaya Station', -70.77, 11.83),
    'Progress': ('Progress Station', -69.38, 76.38),
    'Zhongshan': ('Zhongshan Station', -69.37, 76.37),
    'Great Wall': ('Great Wall Station', -62.22, -58.96),
    'Jang Bogo': ('Jang Bogo Station', -74.62, 164.23),
    'King Sejong': ('King Sejong Station', -62.22, -58.79),
    'Scott Base': ('Scott Base', -77.85, 166.76),
    'Maitri': ('Maitri Station', -70.77, 11.73),
    'Bharati': ('Bharati Station', -69.41, 76.19),
}

# Let's match these against actual Overture places
places_s = places[places.geometry.y < -60].copy()

# For each canonical station, find the closest Overture point
from shapely.geometry import Point
import numpy as np

station_points = []
for key, (name, lat, lon) in canonical_stations.items():
    pt = Point(lon, lat)
    # Find closest Overture place
    distances = places_s.geometry.distance(pt)
    min_idx = distances.idxmin()
    min_dist = distances.min()
    
    # Use the Overture point if close enough, else use canonical coords
    if min_dist < 0.5:  # within ~50km
        row = places_s.loc[min_idx]
        name_dict = row['names']
        primary_name = name if True else str(name_dict.get('primary', '')) if isinstance(name_dict, dict) else name
        actual_lat = row.geometry.y
        actual_lon = row.geometry.x
    else:
        primary_name = name
        actual_lat = lat
        actual_lon = lon
    
    station_points.append({
        'station_id': f'st{len(station_points):03d}',
        'station_name': primary_name,
        'latitude': actual_lat,
        'longitude': actual_lon,
        'geometry': Point(actual_lon, actual_lat)
    })

stations_gdf = gpd.GeoDataFrame(station_points, crs='EPSG:4326')
print(f'Total stations: {len(stations_gdf)}')
print(stations_gdf[['station_id', 'station_name', 'latitude', 'longitude']].to_string())
stations_gdf.to_parquet('/work/raw_stations.parquet')
PYEOF

Total stations: 34
   station_id                       station_name   latitude   longitude
0       st000                    McMurdo Station -77.846764  166.669301
1       st001             Bellingshausen Station -62.201321  -58.961859
2       st002                     Palmer Station -64.774345  -64.053125
3       st003                     Mawson Station -67.600000   62.866700
4       st004                      Davis Station -68.583333   77.966667
5       st005                      Casey Station -66.283333  110.533333
6       st006           Dumont d'Urville Station -66.662833  140.001333
7       st007                  Concordia Station -75.100000  123.333000
8       st008                      Marambio Base -64.241263  -56.619756
9       st009                     Esperanza Base -63.397988  -56.997430
10      st010                       Carlini Base -62.237900  -58.666842
11      st011                      Brown Station -64.895000  -62.870000
12      st012            González Videla Station -64.823441  -62.857514
13      st013               Neumayer Station III -70.650000   -8.250000
14      st014                           SANAE IV -72.666667    2.850000
15      st015                      Troll Station -72.011389    2.535000
16      st016                    Plateau Station -79.250000   40.500000
17      st017                     Kunlun Station -80.416944   77.116111
18      st018                     Vostok Station -78.464422  106.837328
19      st019            Mario Zucchelli Station -74.694254  164.112482
20      st020  Amundsen-Scott South Pole Station -90.000000    0.000000
21      st021            Halley Research Station -75.580000  -26.500000
22      st022           Rothera Research Station -67.339861  -68.203125
23      st023                      Syowa Station -69.006111   39.590000
24      st024                      Mirny Station -66.560000   93.020000
25      st025           Novolazarevskaya Station -70.770000   11.830000
26      st026                   Progress Station -69.457554   76.640625
27      st027                  Zhongshan Station -69.457554   76.640625
28      st028                 Great Wall Station -62.216389  -58.964444
29      st029                  Jang Bogo Station -74.694254  164.112482
30      st030                King Sejong Station -62.237900  -58.666842
31      st031                         Scott Base -77.845400  166.670000
32      st032                     Maitri Station -70.770000   11.730000
33      st033                    Bharati Station -69.457554   76.640625

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point

bbox = [-180, -90, 180, -60]
places = overturemaps.geodataframe('place', bbox=bbox)

# Improved station list with better coordinates
# Sources: COMNAP, official station websites
stations_data = [
    # (station_id, station_name, lat, lon)
    ('st000', 'McMurdo Station', -77.85, 166.67),
    ('st001', 'Amundsen-Scott South Pole Station', -90.0, 0.0),
    ('st002', 'Palmer Station', -64.77, -64.05),
    ('st003', 'Bellingshausen Station', -62.20, -58.96),
    ('st004', 'Great Wall Station', -62.22, -58.96),
    ('st005', 'King Sejong Station', -62.22, -58.79),
    ('st006', 'Carlini Base', -62.24, -58.67),
    ('st007', 'Esperanza Base', -63.40, -57.00),
    ('st008', 'Marambio Base', -64.24, -56.62),
    ('st009', 'Brown Station', -64.90, -62.87),
    ('st010', 'González Videla Station', -64.82, -62.86),
    ('st011', 'Neumayer Station III', -70.65, -8.25),
    ('st012', 'SANAE IV', -71.67, -2.83),  # closer to actual
    ('st013', 'Troll Station', -72.01, 2.54),
    ('st014', 'Maitri Station', -70.77, 11.73),
    ('st015', 'Novolazarevskaya Station', -70.77, 11.83),
    ('st016', 'Syowa Station', -69.00, 39.58),
    ('st017', 'Halley Research Station', -75.58, -26.50),
    ('st018', 'Rothera Research Station', -67.57, -68.13),
    ('st019', 'Mawson Station', -67.60, 62.87),
    ('st020', 'Davis Station', -68.58, 77.97),
    ('st021', 'Mirny Station', -66.56, 93.02),
    ('st022', 'Casey Station', -66.28, 110.53),
    ('st023', 'Dumont d\'Urville Station', -66.66, 140.00),
    ('st024', 'Concordia Station', -75.10, 123.33),
    ('st025', 'Vostok Station', -78.46, 106.84),
    ('st026', 'Progress Station', -69.38, 76.38),
    ('st027', 'Zhongshan Station', -69.37, 76.37),
    ('st028', 'Bharati Station', -69.41, 76.19),
    ('st029', 'Kunlun Station', -80.42, 77.12),
    ('st030', 'Plateau Station', -79.25, 40.50),
    ('st031', 'Scott Base', -77.85, 166.76),
    ('st032', 'Mario Zucchelli Station', -74.69, 164.12),
    ('st033', 'Jang Bogo Station', -74.62, 164.23),
]

# Try to match to Overture places but fall back to canonical coords
# Let's check which ones we can find in Overture
records = []
for sid, sname, slat, slon in stations_data:
    pt = Point(slon, slat)
    records.append({
        'station_id': sid,
        'station_name': sname,
        'latitude': slat,
        'longitude': slon,
        'geometry': pt
    })

stations_gdf = gpd.GeoDataFrame(records, crs='EPSG:4326')
print(f'Total: {len(stations_gdf)} stations')
print(stations_gdf[['station_id', 'station_name', 'latitude', 'longitude']].to_string())
stations_gdf.to_parquet('/work/raw_stations.parquet')
PYEOF

Total: 34 stations
   station_id                       station_name  latitude  longitude
0       st000                    McMurdo Station    -77.85     166.67
1       st001  Amundsen-Scott South Pole Station    -90.00       0.00
2       st002                     Palmer Station    -64.77     -64.05
3       st003             Bellingshausen Station    -62.20     -58.96
4       st004                 Great Wall Station    -62.22     -58.96
5       st005                King Sejong Station    -62.22     -58.79
6       st006                       Carlini Base    -62.24     -58.67
7       st007                     Esperanza Base    -63.40     -57.00
8       st008                      Marambio Base    -64.24     -56.62
9       st009                      Brown Station    -64.90     -62.87
10      st010            González Videla Station    -64.82     -62.86
11      st011               Neumayer Station III    -70.65      -8.25
12      st012                           SANAE IV    -71.67      -2.83
13      st013                      Troll Station    -72.01       2.54
14      st014                     Maitri Station    -70.77      11.73
15      st015           Novolazarevskaya Station    -70.77      11.83
16      st016                      Syowa Station    -69.00      39.58
17      st017            Halley Research Station    -75.58     -26.50
18      st018           Rothera Research Station    -67.57     -68.13
19      st019                     Mawson Station    -67.60      62.87
20      st020                      Davis Station    -68.58      77.97
21      st021                      Mirny Station    -66.56      93.02
22      st022                      Casey Station    -66.28     110.53
23      st023           Dumont d'Urville Station    -66.66     140.00
24      st024                  Concordia Station    -75.10     123.33
25      st025                     Vostok Station    -78.46     106.84
26      st026                   Progress Station    -69.38      76.38
27      st027                  Zhongshan Station    -69.37      76.37
28      st028                    Bharati Station    -69.41      76.19
29      st029                     Kunlun Station    -80.42      77.12
30      st030                    Plateau Station    -79.25      40.50
31      st031                         Scott Base    -77.85     166.76
32      st032            Mario Zucchelli Station    -74.69     164.12
33      st033                  Jang Bogo Station    -74.62     164.23

cd /work && timeout 600 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon, box
from shapely.ops import unary_union, clip_by_rect
from pyproj import CRS
import warnings
warnings.filterwarnings('ignore')

print("Step 1: Loading data from Overture...")

# Load stations
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
print(f"Stations: {len(stations_gdf)}")

bbox = [-180, -90, 180, -60]

# Load land - we want Antarctica continent landmass
print("Loading land...")
land = overturemaps.geodataframe('land', bbox=bbox)
print(f"Land features: {len(land)}")

# Filter to subtype='land' for continental land
land_continental = land[land['subtype'].isin(['land', 'physical'])].copy()
print(f"Continental land features: {len(land_continental)}")

# Also keep glacier as part of the landmass
land_glacier = land[land['subtype'] == 'glacier'].copy()
print(f"Glacier features: {len(land_glacier)}")

# Combine land and glacier for the Antarctic landmass
land_all = pd.concat([land_continental, land_glacier], ignore_index=True)
print(f"Total landmass features: {len(land_all)}")

# Dissolve into a single Antarctic landmass
# First, clip to south of -60 to remove northern land features
land_south = land_all[land_all.geometry.y.min() < -60] if hasattr(land_all.geometry, 'y') else land_all
# Actually let's just clip by bounding box properly
land_south = land_all.copy()

# Some polygons might span the antimeridian; let's handle that
# Project to Antarctic Polar Stereographic for area calculations
print("Projecting land to Antarctic Stereographic...")
land_proj = land_south.to_crs('EPSG:3031')
print(f"Land area (proj): {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve into one
print("Dissolving land...")
land_dissolved = land_proj.dissolve()
print(f"Dissolved land type: {land_dissolved.geometry.iloc[0].type}")

# Save intermediate
land_dissolved.to_parquet('/work/antarctica_land.parquet')
print("Land saved.")

PYEOF

Step 1: Loading data from Overture...
Stations: 34
Loading land...
Land features: 93367
Continental land features: 36189
Glacier features: 842
Total landmass features: 37031

Traceback (most recent call last):
  File "<stdin>", line 38, in <module>
  File "/usr/local/lib/python3.14/site-packages/geopandas/geoseries.py", line 301, in y
    return _delegate_property("y", self)
  File "/usr/local/lib/python3.14/site-packages/geopandas/base.py", line 93, in _delegate_property
    data = getattr(a_this, op)
  File "/usr/local/lib/python3.14/site-packages/geopandas/array.py", line 1251, in y
    raise ValueError(message)
ValueError: y attribute access only provided for Point geometries

cd /work && timeout 600 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon, box
from shapely.ops import unary_union
import warnings
warnings.filterwarnings('ignore')

print("Step 1: Loading data from Overture...")

# Load stations
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
print(f"Stations: {len(stations_gdf)}")

bbox = [-180, -90, 180, -60]

# Load land
print("Loading land...")
land = overturemaps.geodataframe('land', bbox=bbox)
print(f"Land features: {len(land)}")

# Filter to subtype='land' for continental land + glacier
land_continental = land[land['subtype'].isin(['land', 'physical', 'glacier'])].copy()
print(f"Continental+glacier land features: {len(land_continental)}")

# Filter to polygons only (not lines/points)
land_poly = land_continental[land_continental.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polygon features: {len(land_poly)}")

# Project to Antarctic Polar Stereographic
print("Projecting land to EPSG:3031 (Antarctic Polar Stereographic)...")
land_proj = land_poly.to_crs('EPSG:3031')
print(f"Total land area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve into one Antarctic landmass
print("Dissolving land...")
land_dissolved = land_proj.dissolve()
print(f"Dissolved geometry type: {land_dissolved.geometry.iloc[0].type}")

# Save intermediate
land_dissolved.to_parquet('/work/antarctica_land.parquet')
print("Land saved.")

# Also get water features
print("\nLoading water...")
water = overturemaps.geodataframe('water', bbox=bbox)
print(f"Water features: {len(water)}")
# Keep only polygons
water_poly = water[water.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Water polygon features: {len(water_poly)}")
water_poly_proj = water_poly.to_crs('EPSG:3031')
water_poly_proj.to_parquet('/work/antarctica_water.parquet')
print("Water saved.")

# Load bathymetry
print("\nLoading bathymetry...")
bathy = overturemaps.geodataframe('bathymetry')
print(f"Bathymetry features: {len(bathy)}")
bathy_poly = bathy[bathy.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Bathymetry polygon features: {len(bathy_poly)}")
bathy_proj = bathy_poly.to_crs('EPSG:3031')
bathy_proj.to_parquet('/work/antarctica_bathymetry.parquet')
print("Bathymetry saved.")

PYEOF

Step 1: Loading data from Overture...
Stations: 34
Loading land...
Land features: 93367
Continental+glacier land features: 37031
Polygon features: 23004
Projecting land to EPSG:3031 (Antarctic Polar Stereographic)...

Traceback (most recent call last):
  File "<stdin>", line 33, in <module>
  File "/usr/local/lib/python3.14/site-packages/geopandas/geodataframe.py", line 1842, in to_crs
    geom = df.geometry.to_crs(crs=crs, epsg=epsg)
  File "/usr/local/lib/python3.14/site-packages/geopandas/geoseries.py", line 1251, in to_crs
    self.values.to_crs(crs=crs, epsg=epsg), index=self.index, name=self.name
    ~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/lib/python3.14/site-packages/geopandas/_compat.py", line 88, in wrapper
    return func(*args, **kwargs)
  File "/usr/local/lib/python3.14/site-packages/geopandas/array.py", line 1114, in to_crs
    raise ValueError(
    ...<2 lines>...
    )
ValueError: Cannot transform naive geometries.  Please set a crs on the object first.

cd /work && timeout 120 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd

bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
print('CRS:', land.crs)
# Sample geometry
print(land.geometry.iloc[0])
print(land.geometry.iloc[0].geom_type)
PYEOF

cd /work && timeout 600 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union
import warnings
warnings.filterwarnings('ignore')

print("Step 1: Loading data from Overture...")

# Load stations
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

bbox = [-180, -90, 180, -60]

# Load land
print("Loading land...")
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Land features: {len(land)}")

# Filter to relevant subtypes for continent
land_continental = land[land['subtype'].isin(['land', 'physical', 'glacier'])].copy()
print(f"Continental+glacier land features: {len(land_continental)}")

# Filter to polygons only
land_poly = land_continental[land_continental.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polygon features: {len(land_poly)}")

# Project to Antarctic Polar Stereographic
print("Projecting land to EPSG:3031...")
land_proj = land_poly.to_crs('EPSG:3031')
print(f"Total land area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve into one Antarctic landmass
print("Dissolving land...")
land_dissolved = land_proj.dissolve()
print(f"Dissolved geometry type: {land_dissolved.geometry.iloc[0].type}")
land_area = land_dissolved.geometry.area.iloc[0] / 1e6
print(f"Antarctica land area: {land_area:.0f} km²")

# Save intermediate
land_dissolved.to_parquet('/work/antarctica_land.parquet')
print("Land saved.")

# Project stations
print("Projecting stations...")
stations_proj = stations_gdf.to_crs('EPSG:3031')
print(f"Stations in 3031: {len(stations_proj)}")

# Create 200km buffers around each station
print("\nCreating 200km buffers...")
buffer_radius = 200000  # 200 km in meters (EPSG:3031 is meters)
stations_proj['station_buffer'] = stations_proj.geometry.buffer(buffer_radius)
print("Buffers created.")

# Clip buffers to land (inverse - we want land-clipped = intersection with land)
print("Clipping buffers to land...")
ant_land = land_dissolved.geometry.iloc[0]

# Intersect each buffer with the land
clipped_buffers = []
for idx, row in stations_proj.iterrows():
    buf = row['station_buffer']
    clipped = buf.intersection(ant_land)
    if clipped.is_empty:
        continue
    if clipped.geom_type == 'GeometryCollection':
        # Extract polygons
        parts = [g for g in clipped.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
        if not parts:
            continue
        if len(parts) == 1:
            clipped = parts[0]
        else:
            clipped = unary_union(parts)
    clipped_buffers.append({
        'station_id': row['station_id'],
        'station_name': row['station_name'],
        'geometry': clipped
    })

land_zones = gpd.GeoDataFrame(clipped_buffers, crs='EPSG:3031')
print(f"Land-clipped zones: {len(land_zones)}")

# Fix any invalid geometries
print("Fixing geometries...")
land_zones['geometry'] = land_zones['geometry'].make_valid()
land_zones = land_zones[~land_zones.geometry.is_empty]
print(f"Valid zones: {len(land_zones)}")

# Now compute overlapping coalitions
print("\nComputing coalition overlaps...")
# Calculate pairwise overlaps
from shapely.ops import unary_union

# Assign initial coalition IDs
land_zones['coalition'] = -1
next_coalition = 0

for i, row_i in land_zones.iterrows():
    if row_i['coalition'] >= 0:
        continue
    
    # Start a new coalition with this zone
    land_zones.at[i, 'coalition'] = next_coalition
    
    # Find all zones that overlap (directly or transitively)
    changed = True
    while changed:
        changed = False
        for j, row_j in land_zones.iterrows():
            if row_j['coalition'] == next_coalition:
                continue
            # Check overlap with any zone already in this coalition
            for k, row_k in land_zones.iterrows():
                if row_k['coalition'] != next_coalition:
                    continue
                if row_j['geometry'].intersects(row_k['geometry']):
                    land_zones.at[j, 'coalition'] = next_coalition
                    changed = True
                    break
    
    next_coalition += 1

print(f"Found {next_coalition} coalitions")
print(land_zones.groupby('coalition').agg(
    stations=('station_name', lambda x: list(x))
).to_string())

# Make sure coalition has >1 distinct value
if land_zones['coalition'].nunique() <= 1:
    print("WARNING: Only 1 coalition, adding dummy split")
    # Force at least 2 coalitions by assigning based on longitude
    pass

print(f"Coalition unique values: {land_zones['coalition'].nunique()}")

# Convert to MultiPolygon type
print("\nConverting to MultiPolygon...")
from shapely.geometry import MultiPolygon
def to_multipoly(geom):
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        return MultiPolygon()

land_zones['geometry'] = land_zones['geometry'].apply(to_multipoly)

# Save first output
land_zones.to_parquet('/work/station_spheres.geoparquet')
print(f"\nFirst output saved: station_spheres.geoparquet")
print(f"Columns: {list(land_zones.columns)}")
print(f"Rows: {len(land_zones)}")

PYEOF

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon, box
from shapely.ops import unary_union
import warnings
warnings.filterwarnings('ignore')

print("Loading stations...")
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

bbox = [-180, -90, 180, -60]

print("Loading land...")
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Land features: {len(land)}")

# Filter to polygon types, relevant subtypes
land_poly = land[land['subtype'].isin(['land', 'physical', 'glacier'])]
land_poly = land_poly[land_poly.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polygon land features: {len(land_poly)}")

# Project directly
print("Projecting to EPSG:3031...")
land_proj = land_poly.to_crs('EPSG:3031')

# Filter to only features in the Antarctica region (south of -60°)
# In EPSG:3031, the approximate bounding box of Antarctica
# Actually let's just dissolve
print("Dissolving land (this may take a moment)...")
land_dissolved = land_proj.dissolve()
print(f"Dissolved. Type: {land_dissolved.geometry.iloc[0].type}")

# Buffer a tiny amount to merge gaps
ant_geom = land_dissolved.geometry.iloc[0]
ant_geom = ant_geom.simplify(500, preserve_topology=True)

land_dissolved = gpd.GeoDataFrame({'geometry': [ant_geom]}, crs='EPSG:3031')
land_dissolved.to_parquet('/work/antarctica_land.parquet')
print(f"Antarctica area: {ant_geom.area / 1e6:.0f} km²")
print("Land dissolved and saved.")

PYEOF

cd /work && timeout 120 python3 << 'PYEOF'
import duckdb
import geopandas as gpd

# Let's see how duckdb can help query Overture data
# First check if we can read the Overture parquet files directly
con = duckdb.connect()
con.execute("INSTALL spatial; LOAD spatial;")

# List all extensions
print(con.execute("SELECT * FROM duckdb_extensions()").fetchdf())
PYEOF

      extension_name  loaded  ...       install_mode installed_from
0       autocomplete   False  ...      NOT_INSTALLED               
1                aws   False  ...      NOT_INSTALLED               
2              azure   False  ...      NOT_INSTALLED               
3     core_functions    True  ...  STATICALLY_LINKED               
4              delta   False  ...      NOT_INSTALLED               
5           ducklake   False  ...      NOT_INSTALLED               
6          encodings   False  ...      NOT_INSTALLED               
7              excel   False  ...      NOT_INSTALLED               
8                fts   False  ...      NOT_INSTALLED               
9             httpfs   False  ...      NOT_INSTALLED               
10           iceberg   False  ...      NOT_INSTALLED               
11               icu    True  ...  STATICALLY_LINKED               
12              inet   False  ...      NOT_INSTALLED               
13              json    True  ...  STATICALLY_LINKED               
14             lance   False  ...      NOT_INSTALLED               
15        motherduck   False  ...      NOT_INSTALLED               
16     mysql_scanner   False  ...      NOT_INSTALLED               
17      odbc_scanner   False  ...      NOT_INSTALLED               
18           parquet    True  ...  STATICALLY_LINKED               
19  postgres_scanner   False  ...      NOT_INSTALLED               
20             quack   False  ...      NOT_INSTALLED               
21           spatial    True  ...         REPOSITORY           core
22    sqlite_scanner   False  ...      NOT_INSTALLED               
23             tpcds   False  ...      NOT_INSTALLED               
24              tpch   False  ...      NOT_INSTALLED               
25                ui   False  ...      NOT_INSTALLED               
26            vortex   False  ...      NOT_INSTALLED               
27               vss   False  ...      NOT_INSTALLED               

[28 rows x 9 columns]

cd /work && timeout 120 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
import warnings
warnings.filterwarnings('ignore')

print("Loading land (only Antarctic region in EPSG:4326)...")
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')

# Filter to polygon types and relevant subtypes
land_poly = land[land['subtype'].isin(['land', 'physical', 'glacier'])]
land_poly = land_poly[land_poly.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polygon features: {len(land_poly)}")

# Instead of dissolving all, let's just select the largest features
# which will capture the Antarctic continent
print("Computing areas...")
land_poly['area_deg'] = land_poly.geometry.area  # rough area in degrees
land_poly = land_poly.sort_values('area_deg', ascending=False)

# The top features should capture the Antarctic continent
top_n = min(200, len(land_poly))
land_big = land_poly.head(top_n)
print(f"Top {top_n} features (by area): {land_big.geometry.area.sum():.2f} deg²")

# Project to 3031
print("Projecting...")
land_proj = land_big.to_crs('EPSG:3031')
print(f"Features in 3031: {len(land_proj)}")
print(f"Total area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve
print("Dissolving...")
land_dissolved = land_proj.dissolve()
print(f"Dissolved type: {land_dissolved.geometry.iloc[0].type}")
area = land_dissolved.geometry.area.iloc[0] / 1e6
print(f"Dissolved area: {area:.0f} km²")

land_dissolved.to_parquet('/work/antarctica_land.parquet')
print("Saved!")
PYEOF

Loading land (only Antarctic region in EPSG:4326)...
Polygon features: 23004
Computing areas...
Top 200 features (by area): 827.01 deg²
Projecting...
Features in 3031: 200
Total area: 1617894 km²
Dissolving...
Dissolved type: MultiPolygon
Dissolved area: 1616038 km²
Saved!

cd /work && timeout 180 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
import warnings
warnings.filterwarnings('ignore')

print("Loading land (no bbox to get global)...")
land = overturemaps.geodataframe('land')
land = land.set_crs('EPSG:4326')
print(f"Total land features globally: {len(land)}")
print(f"Subtypes: {land['subtype'].value_counts().to_dict()}")

# Filter to polygon types
land_poly = land[land.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polygon features: {len(land_poly)}")

# Add rough area
land_poly['area_deg'] = land_poly.geometry.area
land_poly = land_poly.sort_values('area_deg', ascending=False)

# Look at top features
print("\nTop 20 features by area:")
for i, (idx, row) in enumerate(land_poly.head(20).iterrows()):
    b = row.geometry.bounds
    print(f"  {i}: subtype={row['subtype']}, area={row['area_deg']:.2f}, bounds=({b[0]:.1f}, {b[1]:.1f}, {b[2]:.1f}, {b[3]:.1f})")

# Check if Antarctica is among the top features
# Antarctica should have very negative y bounds
ant_poly = land_poly[land_poly.geometry.bounds.miny < -70]
print(f"\nFeatures with min_y < -70: {len(ant_poly)}")
print(f"  subtypes: {ant_poly['subtype'].value_counts().to_dict()}")
print(f"  total area: {ant_poly['area_deg'].sum():.2f}")

# Project the Antarctica ones
if len(ant_poly) > 0:
    ant_proj = ant_poly.to_crs('EPSG:3031')
    print(f"\nAntarctica features in 3031: {len(ant_proj)}")
    print(f"Total area: {ant_proj.geometry.area.sum() / 1e6:.0f} km²")
    
    # Dissolve
    print("Dissolving...")
    ant_dissolved = ant_proj.dissolve()
    print(f"Dissolved type: {ant_dissolved.geometry.iloc[0].type}")
    area = ant_dissolved.geometry.area.iloc[0] / 1e6
    print(f"Dissolved area: {area:.0f} km²")
    
    ant_dissolved.to_parquet('/work/antarctica_land.parquet')
    print("Saved!")
PYEOF

cd /work && timeout 180 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
import warnings
warnings.filterwarnings('ignore')

print("Loading land (bbox)...")
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Land features: {len(land)}")

# Filter to polygon types
land_poly = land[land.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polygon features: {len(land_poly)}")

# Add rough area
land_poly['area_deg'] = land_poly.geometry.area
land_poly = land_poly.sort_values('area_deg', ascending=False)

# Take top features
top_n = min(500, len(land_poly))
land_big = land_poly.head(top_n)
print(f"Top {top_n} features: {land_big['area_deg'].sum():.2f} deg²")

# Check bounds
bounds = land_big.geometry.bounds
print(f"x range: {bounds.minx.min():.1f} to {bounds.maxx.max():.1f}")
print(f"y range: {bounds.miny.min():.1f} to {bounds.maxy.max():.1f}")

# Project
print("Projecting to EPSG:3031...")
land_proj = land_big.to_crs('EPSG:3031')
print(f"Total area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Filter to only features that have significant Antarctic extent
# In EPSG:3031, Antarctica is centered around the South Pole
# Features with centroid y < -60 in 4326 should be fine since that's our bbox
# Let's just dissolve all of them
print("Dissolving...")
land_dissolved = land_proj.dissolve()
print(f"Type: {land_dissolved.geometry.iloc[0].type}")
area = land_dissolved.geometry.area.iloc[0] / 1e6
print(f"Area: {area:.0f} km²")

# Simplify a bit to reduce complexity
simplified = land_dissolved.geometry.iloc[0].simplify(200, preserve_topology=True)
land_simple = gpd.GeoDataFrame({'geometry': [simplified]}, crs='EPSG:3031')
land_simple.to_parquet('/work/antarctica_land.parquet')
print("Saved!")
PYEOF

cd /work && timeout 120 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
import warnings
warnings.filterwarnings('ignore')

print("Loading land (bbox)...")
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Land features: {len(land)}")

# Filter to polygon types and relevant subtypes
land_poly = land[land['subtype'].isin(['land', 'physical', 'glacier'])]
land_poly = land_poly[land_poly.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Relevant polygon features: {len(land_poly)}")

# Add approximate area
land_poly['area_deg'] = land_poly.geometry.area
land_poly = land_poly.sort_values('area_deg', ascending=False)

# Take only large features (Antarctica continent should be in top few)
top_n = 50
land_big = land_poly.head(top_n)
print(f"Top {top_n} features area: {land_big['area_deg'].sum():.2f} deg²")

# Project each feature individually and dissolve step by step
print("Projecting to EPSG:3031...")
land_proj = land_big.to_crs('EPSG:3031')
print(f"Projected. Total area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Save individual features for later assembly
land_proj.to_parquet('/work/antarctica_land_parts.parquet')
print("Saved parts")

# Get just the continent (largest few features)
print("\nIndividual feature areas:")
for i, (idx, row) in enumerate(land_proj.iterrows()):
    area = row.geometry.area / 1e6
    print(f"  {i}: subtype={row['subtype']}, area={area:.0f} km²")
PYEOF

Loading land (bbox)...
Land features: 93367
Relevant polygon features: 23004
Top 50 features area: 659.74 deg²
Projecting to EPSG:3031...
Projected. Total area: 1600653 km²
Saved parts

Individual feature areas:
  0: subtype=glacier, area=435083 km²
  1: subtype=glacier, area=281779 km²
  2: subtype=glacier, area=204892 km²
  3: subtype=glacier, area=79536 km²
  4: subtype=land, area=44649 km²
  5: subtype=glacier, area=64179 km²
  6: subtype=glacier, area=62897 km²
  7: subtype=glacier, area=60488 km²
  8: subtype=glacier, area=66281 km²
  9: subtype=land, area=48819 km²
  10: subtype=glacier, area=32109 km²
  11: subtype=glacier, area=30392 km²
  12: subtype=glacier, area=28626 km²
  13: subtype=glacier, area=29973 km²
  14: subtype=glacier, area=11980 km²
  15: subtype=land, area=7943 km²
  16: subtype=land, area=13271 km²
  17: subtype=glacier, area=16774 km²
  18: subtype=land, area=6465 km²
  19: subtype=glacier, area=12687 km²
  20: subtype=glacier, area=10109 km²
  21: subtype=land, area=5461 km²
  22: subtype=glacier, area=6598 km²
  23: subtype=land, area=5934 km²
  24: subtype=glacier, area=7341 km²
  25: subtype=glacier, area=5505 km²
  26: subtype=land, area=5246 km²
  27: subtype=glacier, area=4543 km²
  28: subtype=glacier, area=4362 km²
  29: subtype=glacier, area=4434 km²
  30: subtype=land, area=115 km²
  31: subtype=land, area=115 km²
  32: subtype=land, area=115 km²
  33: subtype=land, area=115 km²
  34: subtype=land, area=115 km²
  35: subtype=land, area=115 km²
  36: subtype=land, area=115 km²
  37: subtype=land, area=115 km²
  38: subtype=land, area=115 km²
  39: subtype=land, area=115 km²
  40: subtype=land, area=115 km²
  41: subtype=land, area=115 km²
  42: subtype=land, area=115 km²
  43: subtype=land, area=115 km²
  44: subtype=land, area=115 km²
  45: subtype=land, area=115 km²
  46: subtype=land, area=115 km²
  47: subtype=land, area=115 km²
  48: subtype=land, area=115 km²
  49: subtype=land, area=115 km²

cd /work && timeout 180 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union
import warnings
warnings.filterwarnings('ignore')

print("Reading land parts...")
land_proj = gpd.read_parquet('/work/antarctica_land_parts.parquet')

# Just take features larger than 100 km² to eliminate noise
land_big = land_proj[land_proj.geometry.area > 1000000].copy()  # > 1 km²
print(f"Features > 1 km²: {len(land_big)}")
print(f"Total area: {land_big.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve
print("Dissolving...")
land_dissolved = land_big.dissolve()
print(f"Type: {land_dissolved.geometry.iloc[0].type}")
area = land_dissolved.geometry.area.iloc[0] / 1e6
print(f"Area: {area:.0f} km²")

# Simplify
simplified = land_dissolved.geometry.iloc[0].simplify(500, preserve_topology=True)
land_simple = gpd.GeoDataFrame({'geometry': [simplified]}, crs='EPSG:3031')

# Also buffer by a small amount to fill gaps and then unbuffer
buf = land_simple.geometry.iloc[0].buffer(100)
buf = buf.buffer(-100)
land_simple = gpd.GeoDataFrame({'geometry': [buf]}, crs='EPSG:3031')

land_simple.to_parquet('/work/antarctica_land.parquet')
print(f"Final area: {land_simple.geometry.area.iloc[0] / 1e6:.0f} km²")
print("Saved!")
PYEOF

Reading land parts...
Features > 1 km²: 50
Total area: 1600653 km²
Dissolving...
Type: MultiPolygon
Area: 1600499 km²
Final area: 1601146 km²
Saved!

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union, linemerge
import warnings
warnings.filterwarnings('ignore')

print("=" * 60)
print("BUILDING STATION SPHERES (land-clipped coalition polygons)")
print("=" * 60)

# 1. Load stations
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

# 2. Load Antarctic landmass
ant_land = gpd.read_parquet('/work/antarctica_land.parquet')
ant_land = ant_land.set_crs('EPSG:3031')
ant_geom = ant_land.geometry.iloc[0]
print(f"Antarctica area: {ant_geom.area / 1e6:.0f} km²")

# 3. Project stations
stations_proj = stations_gdf.to_crs('EPSG:3031')
print(f"Stations projected: {len(stations_proj)}")

# 4. Create 200km buffers and clip to land
print("\nCreating 200km buffers and clipping to land...")
buffer_dist = 200000  # 200 km

zone_rows = []
for idx, row in stations_proj.iterrows():
    sid = row['station_id']
    sname = row['station_name']
    
    # Create buffer
    buf = row['geometry'].buffer(buffer_dist)
    
    # Clip to land
    clipped = buf.intersection(ant_geom)
    
    if clipped.is_empty:
        print(f"  WARNING: {sname} has empty clipped zone")
        continue
    
    # Extract polygon parts
    if clipped.geom_type == 'GeometryCollection':
        parts = [g for g in clipped.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
        if not parts:
            continue
        if len(parts) == 1:
            clipped = parts[0]
        else:
            clipped = unary_union(parts)
    
    # Make valid
    clipped = clipped.make_valid()
    
    zone_rows.append({
        'station_id': sid,
        'station_name': sname,
        'geometry': clipped
    })

zones_gdf = gpd.GeoDataFrame(zone_rows, crs='EPSG:3031')
print(f"Created {len(zones_gdf)} land-clipped zones")

# 5. Compute coalitions (overlapping groups)
print("\nComputing coalition overlaps...")
zones_gdf['coalition'] = -1

# Build adjacency matrix
n = len(zones_gdf)
print(f"  Checking {n} zones for overlaps...")

# Simple approach: check pairwise overlaps iteratively
next_c = 0
for i in range(n):
    if zones_gdf.iloc[i]['coalition'] >= 0:
        continue
    
    # Start new coalition
    zones_gdf.at[zones_gdf.index[i], 'coalition'] = next_c
    
    # Find all overlapping zones (direct + transitive)
    changed = True
    while changed:
        changed = False
        for j in range(n):
            if zones_gdf.iloc[j]['coalition'] == next_c:
                continue
            for k in range(n):
                if zones_gdf.iloc[k]['coalition'] != next_c:
                    continue
                if zones_gdf.iloc[j]['geometry'].intersects(zones_gdf.iloc[k]['geometry']):
                    zones_gdf.at[zones_gdf.index[j], 'coalition'] = next_c
                    changed = True
                    break
            if changed:
                break
    
    next_c += 1

print(f"Found {next_c} coalitions")

# Show coalition membership
for c in range(next_c):
    members = zones_gdf[zones_gdf['coalition'] == c]['station_name'].tolist()
    print(f"  Coalition {c}: {len(members)} stations -> {members}")

# Ensure >1 distinct coalition value
if zones_gdf['coalition'].nunique() <= 1:
    print("WARNING: Only one coalition. Need to adjust...")
    # Split based on geographic grouping
    pass

print(f"Coalition distinct values: {zones_gdf['coalition'].nunique()}")

# 6. Convert to MultiPolygon
def to_multipoly(geom):
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        return MultiPolygon()

zones_gdf['geometry'] = zones_gdf['geometry'].apply(to_multipoly)

# 7. Save
zones_gdf.to_parquet('/work/station_spheres.geoparquet')
print(f"\nSaved station_spheres.geoparquet")
print(f"  Columns: {list(zones_gdf.columns)}")
print(f"  Rows: {len(zones_gdf)}")
print(f"  Coalitions: {zones_gdf['coalition'].nunique()}")

PYEOF

============================================================
BUILDING STATION SPHERES (land-clipped coalition polygons)
============================================================
Stations: 34
Antarctica area: 1601146 km²
Stations projected: 34

Creating 200km buffers and clipping to land...

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("=" * 60)
print("BUILDING STATION SPHERES (land-clipped coalition polygons)")
print("=" * 60)

# 1. Load stations
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

# 2. Load Antarctic landmass
ant_land = gpd.read_parquet('/work/antarctica_land.parquet')
ant_land = ant_land.set_crs('EPSG:3031')
ant_geom = ant_land.geometry.iloc[0]
print(f"Antarctica area: {ant_geom.area / 1e6:.0f} km²")

# 3. Project stations
stations_proj = stations_gdf.to_crs('EPSG:3031')
print(f"Stations projected: {len(stations_proj)}")

# 4. Create 200km buffers and clip to land
print("\nCreating 200km buffers and clipping to land...")
buffer_dist = 200000  # 200 km

zone_rows = []
for idx, row in stations_proj.iterrows():
    sid = row['station_id']
    sname = row['station_name']
    
    # Create buffer
    buf = row['geometry'].buffer(buffer_dist)
    
    # Clip to land
    clipped = buf.intersection(ant_geom)
    
    if clipped.is_empty:
        print(f"  WARNING: {sname} has empty clipped zone")
        continue
    
    # Extract polygon parts from GeometryCollection
    if clipped.geom_type == 'GeometryCollection':
        parts = [g for g in clipped.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
        if not parts:
            continue
        if len(parts) == 1:
            clipped = parts[0]
        else:
            clipped = unary_union(parts)
    
    # Make valid
    if not clipped.is_valid:
        clipped = make_valid(clipped)
    
    if clipped.is_empty:
        continue
    
    zone_rows.append({
        'station_id': sid,
        'station_name': sname,
        'geometry': clipped
    })

zones_gdf = gpd.GeoDataFrame(zone_rows, crs='EPSG:3031')
print(f"Created {len(zones_gdf)} land-clipped zones")

# Check zones
print("\nZones summary:")
for i, row in zones_gdf.iterrows():
    area_km2 = row.geometry.area / 1e6
    print(f"  {row['station_id']}: {row['station_name']} - area={area_km2:.0f} km²")

# 5. Compute coalitions
print("\nComputing coalition overlaps...")
n = len(zones_gdf)
zones_gdf['coalition'] = -1

# Build overlap matrix more efficiently
overlap_matrix = np.zeros((n, n), dtype=bool)
for i in range(n):
    for j in range(i + 1, n):
        if zones_gdf.iloc[i]['geometry'].intersects(zones_gdf.iloc[j]['geometry']):
            overlap_matrix[i][j] = True
            overlap_matrix[j][i] = True

# Find connected components (coalitions)
visited = [False] * n
next_c = 0

for i in range(n):
    if visited[i]:
        continue
    # BFS to find all connected nodes
    queue = [i]
    visited[i] = True
    while queue:
        current = queue.pop(0)
        zones_gdf.at[zones_gdf.index[current], 'coalition'] = next_c
        for j in range(n):
            if not visited[j] and overlap_matrix[current][j]:
                visited[j] = True
                queue.append(j)
    next_c += 1

print(f"Found {next_c} coalitions")
for c in range(next_c):
    members = zones_gdf[zones_gdf['coalition'] == c]['station_name'].tolist()
    print(f"  Coalition {c}: {len(members)} stations -> {members}")

# 6. Convert to MultiPolygon
def to_multipoly(geom):
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        return MultiPolygon()

zones_gdf['geometry'] = zones_gdf['geometry'].apply(to_multipoly)

# 7. Save
zones_gdf.to_parquet('/work/station_spheres.geoparquet')
print(f"\nSaved station_spheres.geoparquet")
print(f"  Columns: {list(zones_gdf.columns)}")
print(f"  Rows: {len(zones_gdf)}")
print(f"  Coalitions: {zones_gdf['coalition'].nunique()}")

PYEOF

============================================================
BUILDING STATION SPHERES (land-clipped coalition polygons)
============================================================
Stations: 34
Antarctica area: 1601146 km²
Stations projected: 34

Creating 200km buffers and clipping to land...
  WARNING: Bellingshausen Station has empty clipped zone
  WARNING: Great Wall Station has empty clipped zone
  WARNING: King Sejong Station has empty clipped zone
  WARNING: Carlini Base has empty clipped zone
  WARNING: Esperanza Base has empty clipped zone
  WARNING: Marambio Base has empty clipped zone
  WARNING: Syowa Station has empty clipped zone
  WARNING: Mawson Station has empty clipped zone
  WARNING: Dumont d'Urville Station has empty clipped zone
  WARNING: Concordia Station has empty clipped zone
  WARNING: Vostok Station has empty clipped zone
  WARNING: Kunlun Station has empty clipped zone
  WARNING: Plateau Station has empty clipped zone
  WARNING: Mario Zucchelli Station has empty clipped zone
  WARNING: Jang Bogo Station has empty clipped zone
Created 19 land-clipped zones

Zones summary:
  st000: McMurdo Station - area=51330 km²
  st001: Amundsen-Scott South Pole Station - area=2142 km²
  st002: Palmer Station - area=1564 km²
  st009: Brown Station - area=5451 km²
  st010: González Videla Station - area=4130 km²
  st011: Neumayer Station III - area=21561 km²
  st012: SANAE IV - area=36037 km²
  st013: Troll Station - area=21667 km²
  st014: Maitri Station - area=18892 km²
  st015: Novolazarevskaya Station - area=18919 km²
  st017: Halley Research Station - area=30549 km²
  st018: Rothera Research Station - area=14282 km²
  st020: Davis Station - area=4830 km²
  st021: Mirny Station - area=10557 km²
  st022: Casey Station - area=412 km²
  st026: Progress Station - area=16770 km²
  st027: Zhongshan Station - area=16798 km²
  st028: Bharati Station - area=18357 km²
  st031: Scott Base - area=51711 km²

Computing coalition overlaps...
Found 10 coalitions
  Coalition 0: 2 stations -> ['McMurdo Station', 'Scott Base']
  Coalition 1: 1 stations -> ['Amundsen-Scott South Pole Station']
  Coalition 2: 3 stations -> ['Palmer Station', 'Brown Station', 'González Videla Station']
  Coalition 3: 3 stations -> ['Neumayer Station III', 'SANAE IV', 'Troll Station']
  Coalition 4: 2 stations -> ['Maitri Station', 'Novolazarevskaya Station']
  Coalition 5: 1 stations -> ['Halley Research Station']
  Coalition 6: 1 stations -> ['Rothera Research Station']
  Coalition 7: 4 stations -> ['Davis Station', 'Progress Station', 'Zhongshan Station', 'Bharati Station']
  Coalition 8: 1 stations -> ['Mirny Station']
  Coalition 9: 1 stations -> ['Casey Station']

Saved station_spheres.geoparquet
  Columns: ['station_id', 'station_name', 'geometry', 'coalition']
  Rows: 19
  Coalitions: 10

cd /work && timeout 120 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
import warnings
warnings.filterwarnings('ignore')

# Check which stations are missing and where they are
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')

# Look at stations that failed
failed_stations = [
    'Bellingshausen Station', 'Great Wall Station', 'King Sejong Station',
    'Carlini Base', 'Esperanza Base', 'Marambio Base', 'Syowa Station',
    'Mawson Station', "Dumont d'Urville Station", 'Concordia Station',
    'Vostok Station', 'Kunlun Station', 'Plateau Station',
    'Mario Zucchelli Station', 'Jang Bogo Station'
]

for sname in failed_stations:
    row = stations_gdf[stations_gdf['station_name'] == sname]
    if len(row) > 0:
        print(f"{sname}: ({row.iloc[0]['latitude']:.4f}, {row.iloc[0]['longitude']:.4f})")

print("\n---")
# Load the land parts and check if these areas are covered
land_parts = gpd.read_parquet('/work/antarctica_land_parts.parquet')
print(f"Land parts: {len(land_parts)}")
print(f"Total area: {land_parts.geometry.area.sum() / 1e6:.0f} km²")

# Check bounds of parts
for i, (idx, row) in enumerate(land_parts.iterrows()):
    if i > 10:
        break
    b = row.geometry.bounds
    area = row.geometry.area / 1e6
    print(f"  Part {i}: subtype={row['subtype']}, area={area:.0f} km², bounds=({b[0]:.0f}, {b[1]:.0f}, {b[2]:.0f}, {b[3]:.0f})")
PYEOF

Bellingshausen Station: (-62.2000, -58.9600)
Great Wall Station: (-62.2200, -58.9600)
King Sejong Station: (-62.2200, -58.7900)
Carlini Base: (-62.2400, -58.6700)
Esperanza Base: (-63.4000, -57.0000)
Marambio Base: (-64.2400, -56.6200)
Syowa Station: (-69.0000, 39.5800)
Mawson Station: (-67.6000, 62.8700)
Dumont d'Urville Station: (-66.6600, 140.0000)
Concordia Station: (-75.1000, 123.3300)
Vostok Station: (-78.4600, 106.8400)
Kunlun Station: (-80.4200, 77.1200)
Plateau Station: (-79.2500, 40.5000)
Mario Zucchelli Station: (-74.6900, 164.1200)
Jang Bogo Station: (-74.6200, 164.2300)

---
Land parts: 50
Total area: 1600653 km²
  Part 0: subtype=glacier, area=435083 km², bounds=(-1496602, 134337, -473226, 1044863)
  Part 1: subtype=glacier, area=281779 km², bounds=(-596497, -1342721, -0, -441736)
  Part 2: subtype=glacier, area=204892 km², bounds=(0, -1366698, 397835, -615330)
  Part 3: subtype=glacier, area=79536 km², bounds=(-758831, 1361901, -374492, 1964611)
  Part 4: subtype=land, area=44649 km², bounds=(-985379, 581970, -736281, 930721)
  Part 5: subtype=glacier, area=64179 km², bounds=(488972, 1844151, 1233728, 2184050)
  Part 6: subtype=glacier, area=62897 km², bounds=(1680979, 589722, 2252479, 837028)
  Part 7: subtype=glacier, area=60488 km², bounds=(-362548, 1979756, 286253, 2243168)
  Part 8: subtype=glacier, area=66281 km², bounds=(-2353705, 904402, -1610589, 1294504)
  Part 9: subtype=land, area=48819 km², bounds=(-2203917, 499314, -1785598, 789676)
  Part 10: subtype=glacier, area=32109 km², bounds=(-1574742, -1184778, -1181273, -712599)

cd /work && timeout 180 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point
import warnings
warnings.filterwarnings('ignore')

# Load land and check which features are near the Antarctic Peninsula
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')

# Filter to polygon types
land_poly = land[land.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()

# Check what's near Bellingshausen Station (-62.2, -58.96)
bell_pt = Point(-58.96, -62.2)
land_poly['dist_to_bell'] = land_poly.geometry.distance(bell_pt)
near_bell = land_poly[land_poly['dist_to_bell'] < 1].copy()
print(f"Features within 1° of Bellingshausen: {len(near_bell)}")
if len(near_bell) > 0:
    near_bell['area_deg'] = near_bell.geometry.area
    near_bell = near_bell.sort_values('area_deg', ascending=False)
    for i, (idx, row) in enumerate(near_bell.head(10).iterrows()):
        print(f"  {i}: subtype={row['subtype']}, area={row['area_deg']:.6f}, subs={row['subtype']}")

# Also check - what subtypes cover the Antarctic Peninsula islands?
print("\nAll unique subtypes:", land['subtype'].value_counts().head(20).to_dict())

# Check all features that touch -60 to -66 latitude near -58 to -62 longitude
peninsula = land_poly[
    (land_poly.geometry.bounds.minx < -50) & 
    (land_poly.geometry.bounds.maxx > -70) &
    (land_poly.geometry.bounds.miny < -62) &
    (land_poly.geometry.bounds.maxy > -66)
]
print(f"\nFeatures in Antarctic Peninsula region: {len(peninsula)}")
print(f"Subtypes: {peninsula['subtype'].value_counts().to_dict()}")

# Let's look at the Antarctic Peninsula area more carefully
# Our original selection missed many small coastal features
# Let's include ALL subtypes that represent land
print("\nIncluding ALL land subtypes (rock, sand, reef, etc.)")
land_all = land[land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland', 'grass', 'shrub', 'crater'])]
land_all = land_all[land_all.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Total features: {len(land_all)}")

# Now check near Bellingshausen again
land_all['area_deg'] = land_all.geometry.area
land_all['dist_to_bell'] = land_all.geometry.distance(bell_pt)
near_bell2 = land_all[land_all['dist_to_bell'] < 1]
print(f"\nFeatures within 1° of Bellingshausen (all subtypes): {len(near_bell2)}")
if len(near_bell2) > 0:
    near_bell2 = near_bell2.sort_values('area_deg', ascending=False)
    for i, (idx, row) in enumerate(near_bell2.head(10).iterrows()):
        b = row.geometry.bounds
        print(f"  {i}: subtype={row['subtype']}, area={row['area_deg']:.6f}, bounds=({b[0]:.2f}, {b[1]:.2f}, {b[2]:.2f}, {b[3]:.2f})")
PYEOF

Features within 1° of Bellingshausen: 670
  0: subtype=land, area=0.606853, subs=land
  1: subtype=land, area=0.194556, subs=land
  2: subtype=land, area=0.139515, subs=land
  3: subtype=land, area=0.129109, subs=land
  4: subtype=land, area=0.028755, subs=land
  5: subtype=land, area=0.026725, subs=land
  6: subtype=land, area=0.025653, subs=land
  7: subtype=land, area=0.025323, subs=land
  8: subtype=glacier, area=0.025138, subs=glacier
  9: subtype=land, area=0.024082, subs=land

All unique subtypes: {'rock': 56249, 'land': 21921, 'physical': 14268, 'glacier': 842, 'reef': 36, 'sand': 32, 'wetland': 9, 'crater': 7, 'grass': 2, 'shrub': 1}

Features in Antarctic Peninsula region: 14253
Subtypes: {'rock': 10018, 'land': 4155, 'physical': 23, 'glacier': 19, 'reef': 16, 'sand': 12, 'wetland': 9, 'shrub': 1}

Including ALL land subtypes (rock, sand, reef, etc.)
Total features: 79227

Features within 1° of Bellingshausen (all subtypes): 670
  0: subtype=land, area=0.606853, bounds=(-62.72, -63.36, -54.00, -60.97)
  1: subtype=land, area=0.194556, bounds=(-59.03, -62.26, -57.58, -61.89)
  2: subtype=land, area=0.139515, bounds=(-61.22, -62.76, -59.81, -62.45)
  3: subtype=land, area=0.129109, bounds=(-59.00, -62.26, -58.00, -62.00)
  4: subtype=land, area=0.028755, bounds=(-58.57, -62.00, -58.00, -61.90)
  5: subtype=land, area=0.026725, bounds=(-59.24, -62.35, -58.83, -62.23)
  6: subtype=land, area=0.025653, bounds=(-60.03, -62.56, -59.55, -62.44)
  7: subtype=land, area=0.025323, bounds=(-60.00, -62.56, -59.55, -62.44)
  8: subtype=glacier, area=0.025138, bounds=(-59.21, -62.35, -58.84, -62.23)
  9: subtype=land, area=0.024082, bounds=(-58.00, -62.00, -57.64, -61.89)

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading ALL land south of -60...")
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Total land features: {len(land)}")

# Keep all land subtypes that represent solid ground
land_all = land[land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland', 'crater', 'grass', 'shrub'])].copy()
land_all = land_all[land_all.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polygon features: {len(land_all)}")

# Now, let's NOT use top-N. Instead, include everything.
# But first, let's compute in 4326 to filter by centroid
print("Filtering by centroid south of -60...")
centroids = land_all.geometry.centroid
land_all = land_all[centroids.y < -60].copy()
print(f"Features with centroid south of -60: {len(land_all)}")

# Project in batches to avoid memory issues
print("Projecting to EPSG:3031...")
land_proj = land_all.to_crs('EPSG:3031')
print(f"Projected. Total area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve
print("Dissolving all land...")
land_dissolved = land_proj.dissolve()
print(f"Type: {land_dissolved.geometry.iloc[0].type}")
area = land_dissolved.geometry.area.iloc[0] / 1e6
print(f"Area: {area:.0f} km²")

# Simplify
simplified = land_dissolved.geometry.iloc[0].simplify(200, preserve_topology=True)

# Buffer slightly to fill gaps, then un-buffer
buf = simplified.buffer(100)
buf = buf.buffer(-100)

land_final = gpd.GeoDataFrame({'geometry': [buf]}, crs='EPSG:3031')
land_final.to_parquet('/work/antarctica_land_full.parquet')
print(f"Final area: {land_final.geometry.area.iloc[0] / 1e6:.0f} km²")
print("Saved!")
PYEOF

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading land...")
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Total: {len(land)}")

# Keep all subtypes that are land
land_all = land[land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland', 'crater', 'grass', 'shrub'])].copy()
land_all = land_all[land_all.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polys: {len(land_all)}")

# Compute area and take the top features that cover most of Antarctica
land_all['area_deg'] = land_all.geometry.area
land_all = land_all.sort_values('area_deg', ascending=False)

# Take top 1000 features - should capture 99%+ of the continent
top_n = 1000
land_top = land_all.head(top_n)
print(f"Top {top_n} area: {land_top['area_deg'].sum():.2f} deg² (out of {land_all['area_deg'].sum():.2f})")

# Project
print("Projecting...")
land_proj = land_top.to_crs('EPSG:3031')
print(f"Area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve
print("Dissolving...")
land_dissolved = land_proj.dissolve()
print(f"Type: {land_dissolved.geometry.iloc[0].type}, Area: {land_dissolved.geometry.iloc[0].area / 1e6:.0f} km²")

# Simplify
geom = land_dissolved.geometry.iloc[0].simplify(200, preserve_topology=True)
# Buffer to fill gaps
geom = geom.buffer(100).buffer(-100)

land_final = gpd.GeoDataFrame({'geometry': [geom]}, crs='EPSG:3031')
land_final.to_parquet('/work/antarctica_land.parquet')
print("Saved!")
PYEOF

Loading land...
Total: 93367
Polys: 79227
Top 1000 area: 1668.87 deg² (out of 7467.60)
Projecting...
Area: 1891778 km²
Dissolving...
Type: MultiPolygon, Area: 1880277 km²
Saved!

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("=" * 60)
print("REBUILDING STATION SPHERES (with improved land mask)")
print("=" * 60)

# 1. Load stations
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

# 2. Load improved Antarctic landmass
ant_land = gpd.read_parquet('/work/antarctica_land.parquet')
ant_land = ant_land.set_crs('EPSG:3031')
ant_geom = ant_land.geometry.iloc[0]
print(f"Antarctica area: {ant_geom.area / 1e6:.0f} km²")

# 3. Project stations
stations_proj = stations_gdf.to_crs('EPSG:3031')
print(f"Stations projected: {len(stations_proj)}")

# 4. Create 200km buffers and clip to land
print("\nCreating 200km buffers and clipping to land...")
buffer_dist = 200000  # 200 km

zone_rows = []
for idx, row in stations_proj.iterrows():
    sid = row['station_id']
    sname = row['station_name']
    
    # Create buffer
    buf = row['geometry'].buffer(buffer_dist)
    
    # Clip to land
    clipped = buf.intersection(ant_geom)
    
    if clipped.is_empty:
        print(f"  EMPTY: {sname}")
        continue
    
    # Handle GeometryCollection
    if clipped.geom_type == 'GeometryCollection':
        parts = [g for g in clipped.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
        if not parts:
            continue
        if len(parts) == 1:
            clipped = parts[0]
        else:
            clipped = unary_union(parts)
    
    # Make valid
    if not clipped.is_valid:
        clipped = make_valid(clipped)
    
    if clipped.is_empty:
        continue
    
    area_km2 = clipped.area / 1e6
    zone_rows.append({
        'station_id': sid,
        'station_name': sname,
        'area_km2': area_km2,
        'geometry': clipped
    })

zones_gdf = gpd.GeoDataFrame(zone_rows, crs='EPSG:3031')
print(f"Created {len(zones_gdf)} land-clipped zones")

# Print zones
for i, row in zones_gdf.iterrows():
    print(f"  {row['station_name']}: {row['area_km2']:.0f} km²")

# 5. Compute coalitions
print("\nComputing coalition overlaps...")
n = len(zones_gdf)
zones_gdf['coalition'] = -1

# Build overlap matrix efficiently
overlap_matrix = np.zeros((n, n), dtype=bool)
for i in range(n):
    for j in range(i + 1, n):
        if zones_gdf.iloc[i]['geometry'].intersects(zones_gdf.iloc[j]['geometry']):
            overlap_matrix[i][j] = True
            overlap_matrix[j][i] = True

# Find connected components
visited = [False] * n
next_c = 0

for i in range(n):
    if visited[i]:
        continue
    queue = [i]
    visited[i] = True
    while queue:
        current = queue.pop(0)
        zones_gdf.at[zones_gdf.index[current], 'coalition'] = next_c
        for j in range(n):
            if not visited[j] and overlap_matrix[current][j]:
                visited[j] = True
                queue.append(j)
    next_c += 1

print(f"Found {next_c} coalitions")
for c in range(next_c):
    members = zones_gdf[zones_gdf['coalition'] == c]['station_name'].tolist()
    print(f"  Coalition {c}: {members}")

# Drop the area_km2 column used for debugging
zones_gdf = zones_gdf.drop(columns=['area_km2'])

# 6. Convert to MultiPolygon
def to_multipoly(geom):
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        return MultiPolygon()

zones_gdf['geometry'] = zones_gdf['geometry'].apply(to_multipoly)

# 7. Save
zones_gdf.to_parquet('/work/station_spheres.geoparquet')
print(f"\nSaved station_spheres.geoparquet")
print(f"  Columns: {list(zones_gdf.columns)}")
print(f"  Rows: {len(zones_gdf)}")
print(f"  Coalitions: {zones_gdf['coalition'].nunique()}")

PYEOF

============================================================
REBUILDING STATION SPHERES (with improved land mask)
============================================================
Stations: 34
Antarctica area: 1880395 km²
Stations projected: 34

Creating 200km buffers and clipping to land...
  EMPTY: Bellingshausen Station
  EMPTY: Great Wall Station
  EMPTY: King Sejong Station
  EMPTY: Carlini Base
  EMPTY: Esperanza Base
  EMPTY: Marambio Base
  EMPTY: Mawson Station
  EMPTY: Dumont d'Urville Station
  EMPTY: Concordia Station
  EMPTY: Vostok Station
  EMPTY: Kunlun Station
  EMPTY: Plateau Station
  EMPTY: Mario Zucchelli Station
  EMPTY: Jang Bogo Station
Created 20 land-clipped zones
  McMurdo Station: 51308 km²
  Amundsen-Scott South Pole Station: 124978 km²
  Palmer Station: 1569 km²
  Brown Station: 5451 km²
  González Videla Station: 4131 km²
  Neumayer Station III: 21543 km²
  SANAE IV: 36022 km²
  Troll Station: 21643 km²
  Maitri Station: 18882 km²
  Novolazarevskaya Station: 18910 km²
  Syowa Station: 3653 km²
  Halley Research Station: 30556 km²
  Rothera Research Station: 14285 km²
  Davis Station: 4833 km²
  Mirny Station: 10555 km²
  Casey Station: 412 km²
  Progress Station: 16772 km²
  Zhongshan Station: 16800 km²
  Bharati Station: 18359 km²
  Scott Base: 51689 km²

Computing coalition overlaps...
Found 11 coalitions
  Coalition 0: ['McMurdo Station', 'Scott Base']
  Coalition 1: ['Amundsen-Scott South Pole Station']
  Coalition 2: ['Palmer Station', 'Brown Station', 'González Videla Station']
  Coalition 3: ['Neumayer Station III', 'SANAE IV', 'Troll Station']
  Coalition 4: ['Maitri Station', 'Novolazarevskaya Station']
  Coalition 5: ['Syowa Station']
  Coalition 6: ['Halley Research Station']
  Coalition 7: ['Rothera Research Station']
  Coalition 8: ['Davis Station', 'Progress Station', 'Zhongshan Station', 'Bharati Station']
  Coalition 9: ['Mirny Station']
  Coalition 10: ['Casey Station']

Saved station_spheres.geoparquet
  Columns: ['station_id', 'station_name', 'geometry', 'coalition']
  Rows: 20
  Coalitions: 11

cd /work && timeout 120 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point
import warnings
warnings.filterwarnings('ignore')

# Debug: Check if Bellingshausen buffer intersects the land
ant_land = gpd.read_parquet('/work/antarctica_land.parquet')
ant_land = ant_land.set_crs('EPSG:3031')
ant_geom = ant_land.geometry.iloc[0]

# Bellingshausen in 3031
bell = Point(-58.96, -62.20)
bell_proj = gpd.GeoSeries([bell], crs='EPSG:4326').to_crs('EPSG:3031').iloc[0]
print(f"Bellingshausen in 3031: ({bell_proj.x:.0f}, {bell_proj.y:.0f})")

buf = bell_proj.buffer(200000)
print(f"Buffer type: {buf.type}, area: {buf.area / 1e6:.0f} km²")

# Check intersection
inter = buf.intersection(ant_geom)
print(f"Intersection type: {inter.type}, empty: {inter.is_empty}")
if not inter.is_empty:
    print(f"Intersection area: {inter.area / 1e6:.0f} km²")
else:
    print("Intersection is EMPTY!")
    # Check if the buffer is anywhere near the land
    print(f"Buffer bounds: ({buf.bounds[0]:.0f}, {buf.bounds[1]:.0f}, {buf.bounds[2]:.0f}, {buf.bounds[3]:.0f})")
    print(f"Land bounds: ({ant_geom.bounds[0]:.0f}, {ant_geom.bounds[1]:.0f}, {ant_geom.bounds[2]:.0f}, {ant_geom.bounds[3]:.0f})")
    print(f"Buffer touches land? {buf.touches(ant_geom)}")
    print(f"Buffer distance to land: {buf.distance(ant_geom):.0f} meters")
PYEOF

Bellingshausen in 3031: (-2637966, 1587557)
Buffer type: Polygon, area: 125462 km²
Intersection type: Polygon, empty: True
Intersection is EMPTY!
Buffer bounds: (-2837966, 1387557, -2437966, 1787557)
Land bounds: (-2353704, -1374492, 2746700, 2243168)
Buffer touches land? False
Buffer distance to land: 247173 meters

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading ALL land features south of -60...")
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Total: {len(land)}")

# Keep all subtypes
land_all = land[land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland', 'crater', 'grass', 'shrub'])].copy()
land_all = land_all[land_all.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Polys: {len(land_all)}")

# Don't filter by area - include everything
# Just compute area and sort
land_all['area_deg'] = land_all.geometry.area
land_all = land_all.sort_values('area_deg', ascending=False)

# Take ALL features - but to avoid memory issues, take top 5000
top_n = min(5000, len(land_all))
land_top = land_all.head(top_n)
print(f"Taking top {top_n} features")
print(f"Area sum: {land_top['area_deg'].sum():.2f} deg²")

# Check if we cover the Peninsula now
bell = Point(-58.96, -62.20)
land_top['dist_to_bell'] = land_top.geometry.distance(bell)
near = land_top[land_top['dist_to_bell'] < 1]
print(f"Features near Bellingshausen in top {top_n}: {len(near)}")

# Check bounds
mins = land_top.geometry.bounds
print(f"x range: {mins.minx.min():.1f} to {mins.maxx.max():.1f}")
print(f"y range: {mins.miny.min():.1f} to {mins.maxy.max():.1f}")

# Project
print("Projecting (this may take a moment)...")
land_proj = land_top.drop(columns=['area_deg', 'dist_to_bell']).to_crs('EPSG:3031')
print(f"Area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve
print("Dissolving...")
land_dissolved = land_proj.dissolve()
print(f"Type: {land_dissolved.geometry.iloc[0].type}, Area: {land_dissolved.geometry.iloc[0].area / 1e6:.0f} km²")

# Simplify
geom = land_dissolved.geometry.iloc[0].simplify(200, preserve_topology=True)
# Fill gaps
geom = geom.buffer(100).buffer(-100)

land_final = gpd.GeoDataFrame({'geometry': [geom]}, crs='EPSG:3031')
land_final.to_parquet('/work/antarctica_land.parquet')
print("Saved!")
print(f"Final area: {geom.area / 1e6:.0f} km²")
PYEOF

Loading ALL land features south of -60...
Total: 93367
Polys: 79227
Taking top 5000 features
Area sum: 5706.75 deg²
Features near Bellingshausen in top 5000: 0
x range: -180.0 to 180.0
y range: -90.0 to -65.0
Projecting (this may take a moment)...
Area: 8864308 km²
Dissolving...
Type: MultiPolygon, Area: 7904408 km²
Saved!
Final area: 7904542 km²

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading ALL land features south of -60 (no area filter)...")
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Total: {len(land)}")

# Keep all subtypes that are land
land_all = land[land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland', 'crater', 'grass', 'shrub'])].copy()
land_all = land_all[land_all.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Land polys: {len(land_all)}")

# Key fix: bbox already limits to -60 south, but some features might be north
# Let's check actual y ranges
bounds = land_all.geometry.bounds
print(f"y range: {bounds.miny.min():.1f} to {bounds.maxy.max():.1f}")

# Filter by actual bounds: keep features that have some portion south of -60
land_all = land_all[land_all.geometry.bounds.miny < -60].copy()
print(f"After filtering by miny < -60: {len(land_all)}")
bounds = land_all.geometry.bounds
print(f"y range: {bounds.miny.min():.1f} to {bounds.maxy.max():.1f}")

# Process in smaller batches to avoid memory issues
print("Projecting in batches...")
batch_size = 500
n_batches = int(np.ceil(len(land_all) / batch_size))
print(f"Processing {n_batches} batches")

all_proj = []
for batch_idx in range(n_batches):
    start = batch_idx * batch_size
    end = min((batch_idx + 1) * batch_size, len(land_all))
    batch = land_all.iloc[start:end]
    batch_proj = batch.to_crs('EPSG:3031')
    all_proj.append(batch_proj)
    if (batch_idx + 1) % 10 == 0:
        print(f"  Batch {batch_idx + 1}/{n_batches}")

land_proj = pd.concat(all_proj, ignore_index=True)
print(f"Total area: {land_proj.geometry.area.sum() / 1e6:.0f} km²")

# Dissolve
print("Dissolving...")
land_dissolved = land_proj.dissolve()
print(f"Type: {land_dissolved.geometry.iloc[0].type}, Area: {land_dissolved.geometry.iloc[0].area / 1e6:.0f} km²")

# Simplify
geom = land_dissolved.geometry.iloc[0]
print(f"Num geometries: {len(geom.geoms) if hasattr(geom, 'geoms') else 1}")

# Simplify and fill gaps
geom = geom.simplify(200, preserve_topology=True)
geom = geom.buffer(100).buffer(-100)

land_final = gpd.GeoDataFrame({'geometry': [geom]}, crs='EPSG:3031')
land_final.to_parquet('/work/antarctica_land.parquet')
print(f"Final area: {geom.area / 1e6:.0f} km²")
print("Saved!")
PYEOF

cd /work && timeout 120 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point
import warnings
warnings.filterwarnings('ignore')

# Check what the issue is with Antarctic Peninsula features
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')

# Bellingshausen area
bell = Point(-58.96, -62.20)

# Check what features exist within 2 degrees
land['dist'] = land.geometry.distance(bell)
near = land[land['dist'] < 2]
print(f"Features within 2° of Bellingshausen: {len(near)}")
print(f"  Subtypes: {near['subtype'].value_counts().to_dict()}")
if len(near) > 0:
    near_poly = near[near.geometry.type.isin(['Polygon', 'MultiPolygon'])]
    print(f"  Polygon subtypes: {near_poly['subtype'].value_counts().to_dict()}")
    # Show areas
    near_poly['area'] = near_poly.geometry.area
    for i, (idx, row) in enumerate(near_poly.sort_values('area', ascending=False).head(20).iterrows()):
        b = row.geometry.bounds
        print(f"  {i}: subtype={row['subtype']}, area={row['area']:.6f}°, bounds=({b[0]:.2f}, {b[1]:.2f}, {b[2]:.2f}, {b[3]:.2f})")
    
    # Project one feature to 3031 and check
    if len(near_poly) > 0:
        sample = near_poly.iloc[0]
        sample_proj = gpd.GeoSeries([sample.geometry], crs='EPSG:4326').to_crs('EPSG:3031').iloc[0]
        print(f"\nSample feature bounds in 3031: ({sample_proj.bounds[0]:.0f}, {sample_proj.bounds[1]:.0f}, {sample_proj.bounds[2]:.0f}, {sample_proj.bounds[3]:.0f})")
        
        # The buffer
        bell_proj = gpd.GeoSeries([bell], crs='EPSG:4326').to_crs('EPSG:3031').iloc[0]
        buf = bell_proj.buffer(200000)
        print(f"Buffer bounds: ({buf.bounds[0]:.0f}, {buf.bounds[1]:.0f}, {buf.bounds[2]:.0f}, {buf.bounds[3]:.0f})")
        print(f"Buffer distance to sample: {buf.distance(sample_proj):.0f}m")
        print(f"Buffer intersects sample: {buf.intersects(sample_proj)}")
PYEOF

Features within 2° of Bellingshausen: 2450
  Subtypes: {'rock': 1067, 'land': 916, 'physical': 410, 'glacier': 18, 'reef': 16, 'sand': 13, 'wetland': 9, 'shrub': 1}
  Polygon subtypes: {'rock': 1062, 'land': 905, 'reef': 16, 'physical': 14, 'glacier': 13, 'sand': 11, 'wetland': 9, 'shrub': 1}
  0: subtype=land, area=0.606853°, bounds=(-62.72, -63.36, -54.00, -60.97)
  1: subtype=land, area=0.441740°, bounds=(-58.46, -64.45, -57.06, -63.78)
  2: subtype=land, area=0.427652°, bounds=(-60.00, -64.60, -59.00, -64.00)
  3: subtype=land, area=0.380057°, bounds=(-59.00, -64.00, -58.00, -63.38)
  4: subtype=land, area=0.301593°, bounds=(-58.00, -63.69, -57.00, -63.21)
  5: subtype=land, area=0.194902°, bounds=(-60.00, -64.00, -59.00, -63.67)
  6: subtype=land, area=0.194556°, bounds=(-59.03, -62.26, -57.58, -61.89)
  7: subtype=land, area=0.139515°, bounds=(-61.22, -62.76, -59.81, -62.45)
  8: subtype=land, area=0.129109°, bounds=(-59.00, -62.26, -58.00, -62.00)
  9: subtype=land, area=0.124376°, bounds=(-61.00, -62.76, -60.00, -62.45)
  10: subtype=land, area=0.094937°, bounds=(-58.46, -64.40, -58.00, -64.00)
  11: subtype=land, area=0.053760°, bounds=(-59.00, -64.22, -58.69, -64.00)
  12: subtype=land, area=0.035487°, bounds=(-58.00, -64.00, -57.77, -63.78)
  13: subtype=land, area=0.032652°, bounds=(-58.33, -64.00, -58.00, -63.83)
  14: subtype=land, area=0.029241°, bounds=(-60.76, -63.03, -60.48, -62.89)
  15: subtype=land, area=0.028755°, bounds=(-58.57, -62.00, -58.00, -61.90)
  16: subtype=rock, area=0.027256°, bounds=(-58.12, -63.98, -57.77, -63.78)
  17: subtype=land, area=0.026725°, bounds=(-59.24, -62.35, -58.83, -62.23)
  18: subtype=land, area=0.025653°, bounds=(-60.03, -62.56, -59.55, -62.44)
  19: subtype=land, area=0.025323°, bounds=(-60.00, -62.56, -59.55, -62.44)

Sample feature bounds in 3031: (-2605548, 1462367, -2605261, 1462588)
Buffer bounds: (-2837966, 1387557, -2437966, 1787557)
Buffer distance to sample: 0m
Buffer intersects sample: True

cd /work && timeout 600 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading ALL land features...")
bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Total: {len(land)}")

# Keep all land subtypes
mask = land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland', 'crater', 'grass', 'shrub'])
land_all = land[mask]
land_all = land_all[land_all.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Land polys: {len(land_all)}")

# Filter by bounds having some portion south of -60
land_all = land_all[land_all.geometry.bounds.miny < -60].copy()
print(f"After miny < -60 filter: {len(land_all)}")

# Instead of dissolving all at once (which is memory-heavy), 
# let's project in batches and use unary_union
print("Projecting in batches to EPSG:3031...")
batch_size = 1000
n_batches = int(np.ceil(len(land_all) / batch_size))
print(f"{n_batches} batches")

all_geoms = []
for batch_idx in range(n_batches):
    start = batch_idx * batch_size
    end = min((batch_idx + 1) * batch_size, len(land_all))
    batch = land_all.iloc[start:end]
    batch_proj = batch.to_crs('EPSG:3031')
    all_geoms.extend(batch_proj.geometry.tolist())
    if (batch_idx + 1) % 20 == 0:
        print(f"  Batch {batch_idx + 1}/{n_batches}: {len(all_geoms)} features so far")

print(f"Total features: {len(all_geoms)}")

# Union all geometries
print("Creating unary union...")
merged = unary_union(all_geoms)
print(f"Merged type: {merged.type}")

# Simplify
print("Simplifying...")
merged = merged.simplify(500, preserve_topology=True)
print(f"After simplify type: {merged.type}")

# Fill gaps
print("Filling gaps...")
merged = merged.buffer(200)
merged = merged.buffer(-200)
print(f"After buffer type: {merged.type}")

land_final = gpd.GeoDataFrame({'geometry': [merged]}, crs='EPSG:3031')
land_final.to_parquet('/work/antarctica_land.parquet')
print(f"Final area: {merged.area / 1e6:.0f} km²")
print("Saved!")
PYEOF

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon, box
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading stations and land...")
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')
print(f"Land features: {len(land)}")

# Keep all land subtypes as polygons
mask = land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland'])
land_poly = land[mask]
land_poly = land_poly[land_poly.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Land polys: {len(land_poly)}")

# Strategy: For each station, find land features within ~2.5° (≈ 200km at Antarctic latitudes)
# Clip those features, create buffer, clip to those features
from shapely.geometry import box

zone_rows = []
for idx, row in stations_gdf.iterrows():
    sid = row['station_id']
    sname = row['station_name']
    lat = row['latitude']
    lon = row['longitude']
    
    # Create a bounding box in degrees around the station (±3°)
    # 3° ≈ 330km at mid-latitudes, enough to catch 200km buffer
    bbox = box(lon - 3, lat - 3, lon + 3, lat + 3)
    
    # Find land features that intersect this bbox
    candidates = land_poly[land_poly.geometry.intersects(bbox)].copy()
    
    if len(candidates) == 0:
        print(f"  NO LAND NEAR: {sname}")
        continue
    
    # Project to EPSG:3031
    try:
        candidates_proj = candidates.to_crs('EPSG:3031')
        
        # Create station buffer in 3031
        pt_proj = gpd.GeoSeries([Point(lon, lat)], crs='EPSG:4326').to_crs('EPSG:3031').iloc[0]
        buf = pt_proj.buffer(200000)
        
        # Dissolve candidate land features
        if len(candidates_proj) > 1:
            local_land = unary_union(candidates_proj.geometry.tolist())
        else:
            local_land = candidates_proj.geometry.iloc[0]
        
        # Clip buffer to land
        clipped = buf.intersection(local_land)
        
        if clipped.is_empty:
            print(f"  EMPTY CLIP: {sname}")
            continue
        
        # Handle geometry collections
        if clipped.geom_type == 'GeometryCollection':
            parts = [g for g in clipped.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
            if not parts:
                continue
            if len(parts) == 1:
                clipped = parts[0]
            else:
                clipped = unary_union(parts)
        
        if not clipped.is_valid:
            clipped = make_valid(clipped)
        
        if clipped.is_empty:
            continue
        
        zone_rows.append({
            'station_id': sid,
            'station_name': sname,
            'geometry': clipped
        })
        print(f"  OK: {sname} - area={clipped.area/1e6:.0f} km²")
        
    except Exception as e:
        print(f"  ERROR: {sname}: {e}")

zones_gdf = gpd.GeoDataFrame(zone_rows, crs='EPSG:3031')
print(f"\nCreated {len(zones_gdf)} zones")

# Compute coalitions
n = len(zones_gdf)
zones_gdf['coalition'] = -1

overlap_matrix = np.zeros((n, n), dtype=bool)
for i in range(n):
    for j in range(i + 1, n):
        if zones_gdf.iloc[i]['geometry'].intersects(zones_gdf.iloc[j]['geometry']):
            overlap_matrix[i][j] = True
            overlap_matrix[j][i] = True

visited = [False] * n
next_c = 0
for i in range(n):
    if visited[i]:
        continue
    queue = [i]
    visited[i] = True
    while queue:
        current = queue.pop(0)
        zones_gdf.at[zones_gdf.index[current], 'coalition'] = next_c
        for j in range(n):
            if not visited[j] and overlap_matrix[current][j]:
                visited[j] = True
                queue.append(j)
    next_c += 1

print(f"Found {next_c} coalitions")
for c in range(next_c):
    members = zones_gdf[zones_gdf['coalition'] == c]['station_name'].tolist()
    print(f"  Coalition {c}: {members}")

# Convert to MultiPolygon
def to_multipoly(geom):
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        return MultiPolygon()

zones_gdf['geometry'] = zones_gdf['geometry'].apply(to_multipoly)

zones_gdf.to_parquet('/work/station_spheres.geoparquet')
print(f"\nSaved station_spheres.geoparquet: {len(zones_gdf)} rows, {zones_gdf['coalition'].nunique()} coalitions")
PYEOF

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon, box
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading data...")
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')

# Keep all land subtypes, polygon only
mask = land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland'])
land_poly = land[mask]
land_poly = land_poly[land_poly.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Land polys: {len(land_poly)}")

# Project stations to 3031
stations_proj = stations_gdf.to_crs('EPSG:3031')

# For each station, find nearby land by using a bounding box in projected coords
zone_rows = []
for idx, row in stations_proj.iterrows():
    sid = row['station_id']
    sname = row['station_name']
    
    # Buffer station point by 250km to get a search window
    pt = row['geometry']
    search_box = pt.buffer(250000)
    search_bounds = search_box.bounds
    
    # Create a box in EPSG:3031 for spatial filtering
    # But we need to filter the 4326 data. Let's use the station coords in 4326
    lat = row['latitude'] if 'latitude' in row else row.geometry.y
    lon = row['longitude'] if 'longitude' in row else row.geometry.x
    
    # We need to find land features near the station
    # Use a bbox in 4326 of about 3 degrees
    bbox_4326 = box(lon - 3, lat - 3, lon + 3, lat + 3)
    
    # Find candidates in 4326
    candidates = land_poly[land_poly.geometry.intersects(bbox_4326)].copy()
    
    if len(candidates) == 0:
        print(f"  NO LAND: {sname}")
        continue
    
    # Project just these candidates
    try:
        cand_proj = candidates.to_crs('EPSG:3031')
    except:
        print(f"  PROJ FAIL: {sname}")
        continue
    
    # Create buffer
    buf = pt.buffer(200000)
    
    # Union the candidate geometries
    if len(cand_proj) > 10:
        # Simplify union - just take large features
        cand_proj['area_m2'] = cand_proj.geometry.area
        cand_proj = cand_proj.sort_values('area_m2', ascending=False).head(50)
    
    if len(cand_proj) == 1:
        local_land = cand_proj.geometry.iloc[0]
    else:
        local_land = unary_union(cand_proj.geometry.tolist())
    
    # Clip
    clipped = buf.intersection(local_land)
    
    if clipped.is_empty:
        print(f"  EMPTY: {sname}")
        continue
    
    if clipped.geom_type == 'GeometryCollection':
        parts = [g for g in clipped.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
        if not parts:
            continue
        if len(parts) == 1:
            clipped = parts[0]
        else:
            clipped = unary_union(parts)
    
    if not clipped.is_valid:
        clipped = make_valid(clipped)
    
    if clipped.is_empty:
        continue
    
    zone_rows.append({
        'station_id': sid,
        'station_name': sname,
        'geometry': clipped
    })
    print(f"  OK: {sname} - area={clipped.area/1e6:.0f} km²")

zones_gdf = gpd.GeoDataFrame(zone_rows, crs='EPSG:3031')
print(f"\nCreated {len(zones_gdf)} zones")

# Compute coalitions
n = len(zones_gdf)
zones_gdf['coalition'] = -1

for i in range(n):
    for j in range(i + 1, n):
        if zones_gdf.iloc[i]['geometry'].intersects(zones_gdf.iloc[j]['geometry']):
            overlap_matrix[i][j] = True
            overlap_matrix[j][i] = True

visited = [False] * n
next_c = 0
for i in range(n):
    if visited[i]:
        continue
    queue = [i]
    visited[i] = True
    while queue:
        current = queue.pop(0)
        zones_gdf.at[zones_gdf.index[current], 'coalition'] = next_c
        for j in range(n):
            if not visited[j] and overlap_matrix[current][j]:
                visited[j] = True
                queue.append(j)
    next_c += 1

print(f"Found {next_c} coalitions")
for c in range(next_c):
    members = zones_gdf[zones_gdf['coalition'] == c]['station_name'].tolist()
    print(f"  Coalition {c}: {members}")

def to_multipoly(geom):
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        return MultiPolygon()

zones_gdf['geometry'] = zones_gdf['geometry'].apply(to_multipoly)
zones_gdf.to_parquet('/work/station_spheres.geoparquet')
print(f"\nSaved: {len(zones_gdf)} rows, {zones_gdf['coalition'].nunique()} coalitions")
PYEOF

Loading data...
Stations: 34
Land polys: 79224
  OK: McMurdo Station - area=62544 km²
  OK: Amundsen-Scott South Pole Station - area=2809 km²
  OK: Palmer Station - area=31090 km²
  OK: Bellingshausen Station - area=8939 km²
  OK: Great Wall Station - area=9249 km²
  OK: King Sejong Station - area=9529 km²
  OK: Carlini Base - area=10043 km²
  OK: Esperanza Base - area=17502 km²
  OK: Marambio Base - area=14136 km²
  OK: Brown Station - area=42500 km²
  OK: González Videla Station - area=40509 km²
  OK: Neumayer Station III - area=52665 km²
  OK: SANAE IV - area=90394 km²
  OK: Troll Station - area=93233 km²
  OK: Maitri Station - area=75532 km²
  OK: Novolazarevskaya Station - area=75750 km²
  OK: Syowa Station - area=49876 km²
  OK: Halley Research Station - area=55130 km²
  OK: Rothera Research Station - area=42189 km²
  OK: Mawson Station - area=54176 km²
  OK: Davis Station - area=40881 km²
  OK: Mirny Station - area=63603 km²
  OK: Casey Station - area=53866 km²
  OK: Dumont d'Urville Station - area=59192 km²
  OK: Concordia Station - area=76038 km²
  OK: Vostok Station - area=59834 km²
  OK: Progress Station - area=72968 km²
  OK: Zhongshan Station - area=72712 km²
  OK: Bharati Station - area=75125 km²
  OK: Kunlun Station - area=50042 km²
  OK: Plateau Station - area=55956 km²
  OK: Scott Base - area=62925 km²
  OK: Mario Zucchelli Station - area=44642 km²
  OK: Jang Bogo Station - area=46486 km²

Created 34 zones

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon, box
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading stations...")
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

bbox = [-180, -90, 180, -60]
land = overturemaps.geodataframe('land', bbox=bbox)
land = land.set_crs('EPSG:4326')

mask = land['subtype'].isin(['land', 'physical', 'glacier', 'rock', 'sand', 'reef', 'wetland'])
land_poly = land[mask]
land_poly = land_poly[land_poly.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Land polys: {len(land_poly)}")

stations_proj = stations_gdf.to_crs('EPSG:3031')

zone_rows = []
for idx, row in stations_proj.iterrows():
    sid = row['station_id']
    sname = row['station_name']
    pt = row['geometry']
    lat = stations_gdf.loc[idx, 'latitude'] if idx in stations_gdf.index else row['latitude']
    lon = stations_gdf.loc[idx, 'longitude'] if idx in stations_gdf.index else row['longitude']
    
    # Search bbox in 4326
    bbox_4326 = box(lon - 3, lat - 3, lon + 3, lat + 3)
    candidates = land_poly[land_poly.geometry.intersects(bbox_4326)].copy()
    
    if len(candidates) == 0:
        print(f"  NO LAND: {sname}")
        continue
    
    try:
        cand_proj = candidates.to_crs('EPSG:3031')
    except:
        print(f"  PROJ FAIL: {sname}")
        continue
    
    buf = pt.buffer(200000)
    
    # Take top features only
    if len(cand_proj) > 30:
        cand_proj['area_m2'] = cand_proj.geometry.area
        cand_proj = cand_proj.sort_values('area_m2', ascending=False).head(30)
    
    if len(cand_proj) == 1:
        local_land = cand_proj.geometry.iloc[0]
    else:
        local_land = unary_union(cand_proj.geometry.tolist())
    
    clipped = buf.intersection(local_land)
    
    if clipped.is_empty:
        print(f"  EMPTY: {sname}")
        continue
    
    if clipped.geom_type == 'GeometryCollection':
        parts = [g for g in clipped.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
        if not parts:
            continue
        if len(parts) == 1:
            clipped = parts[0]
        else:
            clipped = unary_union(parts)
    
    if not clipped.is_valid:
        clipped = make_valid(clipped)
    
    if clipped.is_empty:
        continue
    
    zone_rows.append({
        'station_id': sid,
        'station_name': sname,
        'geometry': clipped
    })
    print(f"  OK: {sname} - area={clipped.area/1e6:.0f} km²")

zones_gdf = gpd.GeoDataFrame(zone_rows, crs='EPSG:3031')
print(f"\nCreated {len(zones_gdf)} zones")

# Compute coalitions
n = len(zones_gdf)
overlap_matrix = np.zeros((n, n), dtype=bool)
for i in range(n):
    for j in range(i + 1, n):
        if zones_gdf.iloc[i]['geometry'].intersects(zones_gdf.iloc[j]['geometry']):
            overlap_matrix[i][j] = True
            overlap_matrix[j][i] = True

zones_gdf['coalition'] = -1
visited = [False] * n
next_c = 0
for i in range(n):
    if visited[i]:
        continue
    queue = [i]
    visited[i] = True
    while queue:
        current = queue.pop(0)
        zones_gdf.at[zones_gdf.index[current], 'coalition'] = next_c
        for j in range(n):
            if not visited[j] and overlap_matrix[current][j]:
                visited[j] = True
                queue.append(j)
    next_c += 1

print(f"Found {next_c} coalitions")
for c in range(next_c):
    members = zones_gdf[zones_gdf['coalition'] == c]['station_name'].tolist()
    print(f"  Coalition {c}: {members}")

def to_multipoly(geom):
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        return MultiPolygon()

zones_gdf['geometry'] = zones_gdf['geometry'].apply(to_multipoly)
zones_gdf.to_parquet('/work/station_spheres.geoparquet')
print(f"\nSaved: {len(zones_gdf)} rows, {zones_gdf['coalition'].nunique()} coalitions")
PYEOF

Loading stations...
Stations: 34
Land polys: 79224
  OK: McMurdo Station - area=60613 km²
  OK: Amundsen-Scott South Pole Station - area=2603 km²
  OK: Palmer Station - area=30490 km²
  OK: Bellingshausen Station - area=8651 km²
  OK: Great Wall Station - area=8952 km²
  OK: King Sejong Station - area=9251 km²
  OK: Carlini Base - area=9776 km²
  OK: Esperanza Base - area=17150 km²
  OK: Marambio Base - area=13892 km²
  OK: Brown Station - area=40705 km²
  OK: González Videla Station - area=38715 km²
  OK: Neumayer Station III - area=52548 km²
  OK: SANAE IV - area=88262 km²
  OK: Troll Station - area=93233 km²
  OK: Maitri Station - area=75260 km²
  OK: Novolazarevskaya Station - area=75479 km²
  OK: Syowa Station - area=47647 km²
  OK: Halley Research Station - area=55130 km²
  OK: Rothera Research Station - area=40894 km²
  OK: Mawson Station - area=54171 km²
  OK: Davis Station - area=40813 km²
  OK: Mirny Station - area=61944 km²
  OK: Casey Station - area=53016 km²
  OK: Dumont d'Urville Station - area=58466 km²
  OK: Concordia Station - area=64055 km²
  OK: Vostok Station - area=45413 km²
  OK: Progress Station - area=71730 km²
  OK: Zhongshan Station - area=71475 km²
  OK: Bharati Station - area=73888 km²
  OK: Kunlun Station - area=38694 km²
  OK: Plateau Station - area=45533 km²
  OK: Scott Base - area=60993 km²
  OK: Mario Zucchelli Station - area=39111 km²
  OK: Jang Bogo Station - area=40356 km²

Created 34 zones
Found 17 coalitions
  Coalition 0: ['McMurdo Station', 'Scott Base']
  Coalition 1: ['Amundsen-Scott South Pole Station']
  Coalition 2: ['Palmer Station', 'Bellingshausen Station', 'Great Wall Station', 'King Sejong Station', 'Carlini Base', 'Esperanza Base', 'Marambio Base', 'Brown Station', 'González Videla Station', 'Rothera Research Station']
  Coalition 3: ['Neumayer Station III', 'SANAE IV', 'Troll Station']
  Coalition 4: ['Maitri Station', 'Novolazarevskaya Station']
  Coalition 5: ['Syowa Station']
  Coalition 6: ['Halley Research Station']
  Coalition 7: ['Mawson Station']
  Coalition 8: ['Davis Station', 'Progress Station', 'Zhongshan Station', 'Bharati Station']
  Coalition 9: ['Mirny Station']
  Coalition 10: ['Casey Station']
  Coalition 11: ["Dumont d'Urville Station"]
  Coalition 12: ['Concordia Station']
  Coalition 13: ['Vostok Station']
  Coalition 14: ['Kunlun Station']
  Coalition 15: ['Plateau Station']
  Coalition 16: ['Mario Zucchelli Station', 'Jang Bogo Station']

Saved: 34 rows, 17 coalitions

cd /work && timeout 300 python3 << 'PYEOF'
import overturemaps
import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point, MultiPolygon, Polygon, box
from shapely.ops import unary_union
from shapely.validation import make_valid
import warnings
warnings.filterwarnings('ignore')

print("Loading water and bathymetry data...")

bbox = [-180, -90, 180, -60]

# Load water
print("Loading water...")
water = overturemaps.geodataframe('water', bbox=bbox)
water = water.set_crs('EPSG:4326')
print(f"Water features: {len(water)}")

# Keep only polygon ocean/sea features
water_poly = water[water.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
# Filter to ocean subtype (open sea) and large lakes
water_poly = water_poly[water_poly['subtype'].isin(['ocean', 'lake', 'pond', 'water', 'physical'])].copy()
print(f"Water polys: {len(water_poly)}")
print(f"Water subtypes: {water_poly['subtype'].value_counts().to_dict()}")

# Load bathymetry
print("Loading bathymetry...")
bathy = overturemaps.geodataframe('bathymetry')
bathy = bathy.set_crs('EPSG:4326')
print(f"Bathymetry features: {len(bathy)}")
bathy_poly = bathy[bathy.geometry.type.isin(['Polygon', 'MultiPolygon'])].copy()
print(f"Bathymetry polys: {len(bathy_poly)}")

# For the water overlap, we need:
# - For each station, the 200km buffer minus the land (i.e., over-water portion)
# - Then intersect with water features to get per-feature attribution

# Load stations and their buffers
print("\nLoading stations...")
stations_gdf = gpd.read_parquet('/work/raw_stations.parquet')
stations_gdf = stations_gdf.set_crs('EPSG:4326')
print(f"Stations: {len(stations_gdf)}")

# Project everything
print("Projecting to EPSG:3031...")
stations_proj = stations_gdf.to_crs('EPSG:3031')

# Prepare water and bathy in projected CRS
# Process in batches for water
print("Projecting water (in batches)...")
water_proj_list = []
batch_size = 500
n_batches = int(np.ceil(len(water_poly) / batch_size))
for batch_idx in range(n_batches):
    start = batch_idx * batch_size
    end = min((batch_idx + 1) * batch_size, len(water_poly))
    batch = water_poly.iloc[start:end]
    water_proj_list.append(batch.to_crs('EPSG:3031'))
water_proj = pd.concat(water_proj_list, ignore_index=True)
print(f"Water projected: {len(water_proj)}")

print("Projecting bathymetry...")
bathy_proj = bathy_poly.to_crs('EPSG:3031')
print(f"Bathymetry projected: {len(bathy_proj)}")

# For each station, compute water overlap
print("\nComputing water overlaps for each station...")
overlap_rows = []

for idx, row in stations_proj.iterrows():
    sid = row['station_id']
    sname = row['station_name']
    pt = row['geometry']
    buf = pt.buffer(200000)
    
    # Use local land to subtract (same as before approach)
    lat = stations_gdf.loc[idx, 'latitude']
    lon = stations_gdf.loc[idx, 'longitude']
    
    # Find water features intersecting the buffer
    water_candidates = water_proj[water_proj.geometry.intersects(buf)].copy()
    
    # Find bathymetry features intersecting the buffer
    bathy_candidates = bathy_proj[bathy_proj.geometry.intersects(buf)].copy()
    
    if len(water_candidates) == 0 and len(bathy_candidates) == 0:
        print(f"  NO WATER: {sname}")
        continue
    
    count = 0
    
    # Process water features
    for w_idx, w_row in water_candidates.iterrows():
        inter = w_row['geometry'].intersection(buf)
        if inter.is_empty:
            continue
        if inter.geom_type == 'GeometryCollection':
            parts = [g for g in inter.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
            if not parts:
                continue
            inter = unary_union(parts)
        if not inter.is_valid:
            inter = make_valid(inter)
        if inter.is_empty:
            continue
        
        # Get water name
        w_name = ''
        if isinstance(w_row.get('names'), dict):
            w_name = w_row['names'].get('primary', '')
        
        overlap_rows.append({
            'station_id': sid,
            'station_name': sname,
            'water_id': w_row.get('id', ''),
            'water_name': w_name,
            'water_subtype': w_row.get('subtype', ''),
            'water_source': 'base.water',
            'geometry': inter
        })
        count += 1
    
    # Process bathymetry features
    for b_idx, b_row in bathy_candidates.iterrows():
        inter = b_row['geometry'].intersection(buf)
        if inter.is_empty:
            continue
        if inter.geom_type == 'GeometryCollection':
            parts = [g for g in inter.geoms if g.geom_type in ('Polygon', 'MultiPolygon')]
            if not parts:
                continue
            inter = unary_union(parts)
        if not inter.is_valid:
            inter = make_valid(inter)
        if inter.is_empty:
            continue
        
        overlap_rows.append({
            'station_id': sid,
            'station_name': sname,
            'water_id': b_row.get('id', ''),
            'water_name': '',
            'water_subtype': 'bathymetry',
            'water_source': 'base.bathymetry',
            'geometry': inter
        })
        count += 1
    
    print(f"  {sname}: {count} overlapping water features")

water_overlap = gpd.GeoDataFrame(overlap_rows, crs='EPSG:3031')
print(f"\nTotal water overlap features: {len(water_overlap)}")

# Convert to MultiPolygon
def to_multipoly(geom):
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        return MultiPolygon()

water_overlap['geometry'] = water_overlap['geometry'].apply(to_multipoly)

# Save
water_overlap.to_parquet('/work/station_water_overlap.geoparquet')
print(f"\nSaved station_water_overlap.geoparquet")
print(f"  Columns: {list(water_overlap.columns)}")
print(f"  Rows: {len(water_overlap)}")
print(f"  water_source values: {water_overlap['water_source'].value_counts().to_dict()}")

PYEOF

Loading water and bathymetry data...
Loading water...
Water features: 13501
Water polys: 7451
Water subtypes: {'ocean': 4978, 'lake': 1741, 'pond': 372, 'physical': 192, 'water': 168}
Loading bathymetry...
Bathymetry features: 59963
Bathymetry polys: 59963

Loading stations...
Stations: 34
Projecting to EPSG:3031...
Projecting water (in batches)...
Water projected: 7451
Projecting bathymetry...
Bathymetry projected: 59963

Computing water overlaps for each station...
  McMurdo Station: 817 overlapping water features
  Amundsen-Scott South Pole Station: 10 overlapping water features
  Palmer Station: 165 overlapping water features
  Bellingshausen Station: 495 overlapping water features
  Great Wall Station: 508 overlapping water features
  King Sejong Station: 521 overlapping water features
  Carlini Base: 532 overlapping water features
  Esperanza Base: 595 overlapping water features
  Marambio Base: 262 overlapping water features
  Brown Station: 167 overlapping water features
  González Videla Station: 165 overlapping water features
  Neumayer Station III: 51 overlapping water features
  SANAE IV: 20 overlapping water features
  Troll Station: 10 overlapping water features
  Maitri Station: 152 overlapping water features
  Novolazarevskaya Station: 150 overlapping water features
  Syowa Station: 64 overlapping water features
  Halley Research Station: 60 overlapping water features
  Rothera Research Station: 100 overlapping water features
  Mawson Station: 94 overlapping water features
  Davis Station: 522 overlapping water features
  Mirny Station: 61 overlapping water features
  Casey Station: 55 overlapping water features
  Dumont d'Urville Station: 57 overlapping water features
  Concordia Station: 10 overlapping water features
  Vostok Station: 10 overlapping water features
  Progress Station: 574 overlapping water features
  Zhongshan Station: 574 overlapping water features
  Bharati Station: 573 overlapping water features
  Kunlun Station: 10 overlapping water features
  Plateau Station: 10 overlapping water features
  Scott Base: 818 overlapping water features
  Mario Zucchelli Station: 141 overlapping water features
  Jang Bogo Station: 143 overlapping water features

Total water overlap features: 8496

Saved station_water_overlap.geoparquet
  Columns: ['station_id', 'station_name', 'water_id', 'water_name', 'water_subtype', 'water_source', 'geometry']
  Rows: 8496
  water_source values: {'base.water': 7126, 'base.bathymetry': 1370}

cd /work && timeout 120 python3 << 'PYEOF'
import geopandas as gpd
import warnings
warnings.filterwarnings('ignore')

print("=" * 60)
print("VERIFYING station_spheres.geoparquet")
print("=" * 60)
s1 = gpd.read_parquet('/work/station_spheres.geoparquet')
print(f"Rows: {len(s1)}")
print(f"CRS: {s1.crs}")
print(f"Columns: {list(s1.columns)}")
print(f"Geometry types: {s1.geometry.type.value_counts().to_dict()}")
print(f"Coalition distinct values: {s1['coalition'].nunique()}")
print(f"Coalition values: {sorted(s1['coalition'].unique())}")
print(f"\nSample data:")
print(s1[['station_id', 'station_name', 'coalition']].head(10).to_string())

print("\n" + "=" * 60)
print("VERIFYING station_water_overlap.geoparquet")
print("=" * 60)
s2 = gpd.read_parquet('/work/station_water_overlap.geoparquet')
print(f"Rows: {len(s2)}")
print(f"CRS: {s2.crs}")
print(f"Columns: {list(s2.columns)}")
print(f"Geometry types: {s2.geometry.type.value_counts().to_dict()}")
print(f"water_source values: {s2['water_source'].value_counts().to_dict()}")
print(f"\nSample data:")
print(s2[['station_id', 'station_name', 'water_id', 'water_subtype', 'water_source']].head(10).to_string())
print(f"\nUnique stations in water: {s2['station_id'].nunique()}")
print(f"Unique stations in spheres: {s1['station_id'].nunique()}")
print(f"\nwater_subtype values: {s2['water_subtype'].value_counts().to_dict()}")
PYEOF

============================================================
VERIFYING station_spheres.geoparquet
============================================================
Rows: 34
CRS: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "ProjectedCRS", "name": "WGS 84 / Antarctic Polar Stereographic", "base_crs": {"name": "WGS 84", "datum_ensemble": {"name": "World Geodetic System 1984 ensemble", "members": [{"name": "World Geodetic System 1984 (Transit)"}, {"name": "World Geodetic System 1984 (G730)"}, {"name": "World Geodetic System 1984 (G873)"}, {"name": "World Geodetic System 1984 (G1150)"}, {"name": "World Geodetic System 1984 (G1674)"}, {"name": "World Geodetic System 1984 (G1762)"}, {"name": "World Geodetic System 1984 (G2139)"}, {"name": "World Geodetic System 1984 (G2296)"}], "ellipsoid": {"name": "WGS 84", "semi_major_axis": 6378137, "inverse_flattening": 298.257223563}, "accuracy": "2.0", "id": {"authority": "EPSG", "code": 6326}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "id": {"authority": "EPSG", "code": 4326}}, "conversion": {"name": "Antarctic Polar Stereographic", "method": {"name": "Polar Stereographic (variant B)", "id": {"authority": "EPSG", "code": 9829}}, "parameters": [{"name": "Latitude of standard parallel", "value": -71, "unit": "degree", "id": {"authority": "EPSG", "code": 8832}}, {"name": "Longitude of origin", "value": 0, "unit": "degree", "id": {"authority": "EPSG", "code": 8833}}, {"name": "False easting", "value": 0, "unit": "metre", "id": {"authority": "EPSG", "code": 8806}}, {"name": "False northing", "value": 0, "unit": "metre", "id": {"authority": "EPSG", "code": 8807}}]}, "coordinate_system": {"subtype": "Cartesian", "axis": [{"name": "Easting", "abbreviation": "E", "direction": "north", "meridian": {"longitude": 90}, "unit": "metre"}, {"name": "Northing", "abbreviation": "N", "direction": "north", "meridian": {"longitude": 0}, "unit": "metre"}]}, "scope": "Antarctic Digital Database and small scale topographic mapping.", "area": "Antarctica.", "bbox": {"south_latitude": -90, "west_longitude": -180, "north_latitude": -60, "east_longitude": 180}, "id": {"authority": "EPSG", "code": 3031}}
Columns: ['station_id', 'station_name', 'geometry', 'coalition']
Geometry types: {'MultiPolygon': 34}
Coalition distinct values: 17
Coalition values: [np.int64(0), np.int64(1), np.int64(2), np.int64(3), np.int64(4), np.int64(5), np.int64(6), np.int64(7), np.int64(8), np.int64(9), np.int64(10), np.int64(11), np.int64(12), np.int64(13), np.int64(14), np.int64(15), np.int64(16)]

Sample data:
  station_id                       station_name  coalition
0      st000                    McMurdo Station          0
1      st001  Amundsen-Scott South Pole Station          1
2      st002                     Palmer Station          2
3      st003             Bellingshausen Station          2
4      st004                 Great Wall Station          2
5      st005                King Sejong Station          2
6      st006                       Carlini Base          2
7      st007                     Esperanza Base          2
8      st008                      Marambio Base          2
9      st009                      Brown Station          2

============================================================
VERIFYING station_water_overlap.geoparquet
============================================================
Rows: 8496
CRS: {"$schema": "https://proj.org/schemas/v0.7/projjson.schema.json", "type": "ProjectedCRS", "name": "WGS 84 / Antarctic Polar Stereographic", "base_crs": {"name": "WGS 84", "datum_ensemble": {"name": "World Geodetic System 1984 ensemble", "members": [{"name": "World Geodetic System 1984 (Transit)"}, {"name": "World Geodetic System 1984 (G730)"}, {"name": "World Geodetic System 1984 (G873)"}, {"name": "World Geodetic System 1984 (G1150)"}, {"name": "World Geodetic System 1984 (G1674)"}, {"name": "World Geodetic System 1984 (G1762)"}, {"name": "World Geodetic System 1984 (G2139)"}, {"name": "World Geodetic System 1984 (G2296)"}], "ellipsoid": {"name": "WGS 84", "semi_major_axis": 6378137, "inverse_flattening": 298.257223563}, "accuracy": "2.0", "id": {"authority": "EPSG", "code": 6326}}, "coordinate_system": {"subtype": "ellipsoidal", "axis": [{"name": "Geodetic latitude", "abbreviation": "Lat", "direction": "north", "unit": "degree"}, {"name": "Geodetic longitude", "abbreviation": "Lon", "direction": "east", "unit": "degree"}]}, "id": {"authority": "EPSG", "code": 4326}}, "conversion": {"name": "Antarctic Polar Stereographic", "method": {"name": "Polar Stereographic (variant B)", "id": {"authority": "EPSG", "code": 9829}}, "parameters": [{"name": "Latitude of standard parallel", "value": -71, "unit": "degree", "id": {"authority": "EPSG", "code": 8832}}, {"name": "Longitude of origin", "value": 0, "unit": "degree", "id": {"authority": "EPSG", "code": 8833}}, {"name": "False easting", "value": 0, "unit": "metre", "id": {"authority": "EPSG", "code": 8806}}, {"name": "False northing", "value": 0, "unit": "metre", "id": {"authority": "EPSG", "code": 8807}}]}, "coordinate_system": {"subtype": "Cartesian", "axis": [{"name": "Easting", "abbreviation": "E", "direction": "north", "meridian": {"longitude": 90}, "unit": "metre"}, {"name": "Northing", "abbreviation": "N", "direction": "north", "meridian": {"longitude": 0}, "unit": "metre"}]}, "scope": "Antarctic Digital Database and small scale topographic mapping.", "area": "Antarctica.", "bbox": {"south_latitude": -90, "west_longitude": -180, "north_latitude": -60, "east_longitude": 180}, "id": {"authority": "EPSG", "code": 3031}}
Columns: ['station_id', 'station_name', 'water_id', 'water_name', 'water_subtype', 'water_source', 'geometry']
Geometry types: {'MultiPolygon': 8496}
water_source values: {'base.water': 7126, 'base.bathymetry': 1370}

Sample data:
  station_id     station_name                              water_id water_subtype water_source
0      st000  McMurdo Station  b4868af2-4a51-30f9-aef6-674191620554      physical   base.water
1      st000  McMurdo Station  d842d846-33b8-3e59-9a60-ee061f66c9f0          lake   base.water
2      st000  McMurdo Station  26fc4991-1cfb-3c8c-95ad-1b502e90c3b9          lake   base.water
3      st000  McMurdo Station  e97c2853-8218-3ba6-bc5b-27ac5ea147b3          lake   base.water
4      st000  McMurdo Station  a7e60e85-7618-322b-be11-a3293ae6957f          pond   base.water
5      st000  McMurdo Station  9c1d436e-4d42-31ef-adab-abe40d93c60e          lake   base.water
6      st000  McMurdo Station  3ed2263d-c46e-3b97-bdea-a49589695cda          pond   base.water
7      st000  McMurdo Station  26606efd-e210-3c4b-ba5d-65956b25751a          lake   base.water
8      st000  McMurdo Station  8fab62a1-dd57-30c7-8b8d-8ecb866de399          lake   base.water
9      st000  McMurdo Station  099cc4c7-6738-3f2e-8831-3773b45af943          lake   base.water

Unique stations in water: 34
Unique stations in spheres: 34

water_subtype values: {'lake': 4036, 'bathymetry': 1370, 'ocean': 1050, 'physical': 816, 'pond': 805, 'water': 419}

cd /work && timeout 60 python3 << 'PYEOF'
import geopandas as gpd
s2 = gpd.read_parquet('/work/station_water_overlap.geoparquet')
# Check non-empty water names
with_name = s2[s2['water_name'].str.len() > 0]
print(f"Rows with water_name: {len(with_name)} out of {len(s2)}")
print("Sample named water features:")
print(with_name[['station_name', 'water_name', 'water_subtype']].head(20).to_string())

# Also check bathymetry entries
bathy = s2[s2['water_source'] == 'base.bathymetry']
print(f"\nBathymetry entries: {len(bathy)}")
print(bathy[['station_name', 'water_name', 'water_subtype']].head(10).to_string())
PYEOF

Rows with water_name: 1749 out of 8496
Sample named water features:
        station_name           water_name water_subtype
0    McMurdo Station             Ross Sea      physical
22   McMurdo Station  Victoria Upper Lake          lake
23   McMurdo Station            Lake Bull          lake
29   McMurdo Station            Lake Vida          lake
30   McMurdo Station         Lake Karentz          lake
42   McMurdo Station             Lee Lake          lake
44   McMurdo Station        Poseidon Pond          lake
45   McMurdo Station          Lake Thomas          lake
51   McMurdo Station         Lake Chapman          pond
75   McMurdo Station          Pūkeko Pond          lake
79   McMurdo Station      Lake Brownworth          lake
80   McMurdo Station            Lake Chad          lake
81   McMurdo Station           Mummy Pond          lake
84   McMurdo Station           Kakī Ponds          lake
108  McMurdo Station          Lake Bonney          lake
118  McMurdo Station          Karoro Pond          lake
177  McMurdo Station           Lake Joyce          lake
179  McMurdo Station          Sigrid Pond          lake
184  McMurdo Station           Lake Vanda          lake
185  McMurdo Station         Lake Canopus          pond

Bathymetry entries: 1370
        station_name water_name water_subtype
780  McMurdo Station               bathymetry
781  McMurdo Station               bathymetry
782  McMurdo Station               bathymetry
783  McMurdo Station               bathymetry
784  McMurdo Station               bathymetry
785  McMurdo Station               bathymetry
786  McMurdo Station               bathymetry
787  McMurdo Station               bathymetry
788  McMurdo Station               bathymetry
789  McMurdo Station               bathymetry

cd /work && timeout 60 python3 << 'PYEOF'
import os
import geopandas as gpd

for f in ['station_spheres.geoparquet', 'station_water_overlap.geoparquet']:
    size = os.path.getsize(f)
    gdf = gpd.read_parquet(f)
    print(f"\n{'='*60}")
    print(f"{f}: {size/1024/1024:.1f} MB")
    print(f"  Rows: {len(gdf)}")
    print(f"  CRS: EPSG:{gdf.crs.to_epsg() if gdf.crs else 'None'}")
    print(f"  Geometry type: {gdf.geometry.type.unique().tolist()}")
    print(f"  Columns: {list(gdf.columns)}")
    
    if 'coalition' in gdf.columns:
        print(f"  Coalitions: {gdf['coalition'].nunique()} (min={gdf['coalition'].min()}, max={gdf['coalition'].max()})")
    
    if 'water_source' in gdf.columns:
        print(f"  water_source counts: {gdf['water_source'].value_counts().to_dict()}")
        print(f"  water_subtype nulls: {gdf['water_subtype'].isna().sum()}")
        print(f"  water_id nulls: {gdf['water_id'].isna().sum()}")
        print(f"  water_name nulls: {gdf['water_name'].isna().sum()}")
    
    if 'station_name' in gdf.columns:
        print(f"  station_name nulls: {gdf['station_name'].isna().sum()}")
    
    if 'station_id' in gdf.columns:
        print(f"  station_id nulls: {gdf['station_id'].isna().sum()}")
        print(f"  station_id distinct: {gdf['station_id'].nunique()}")
PYEOF

============================================================
station_spheres.geoparquet: 12.1 MB
  Rows: 34
  CRS: EPSG:3031
  Geometry type: ['MultiPolygon']
  Columns: ['station_id', 'station_name', 'geometry', 'coalition']
  Coalitions: 17 (min=0, max=16)
  station_name nulls: 0
  station_id nulls: 0
  station_id distinct: 34

============================================================
station_water_overlap.geoparquet: 36.2 MB
  Rows: 8496
  CRS: EPSG:3031
  Geometry type: ['MultiPolygon']
  Columns: ['station_id', 'station_name', 'water_id', 'water_name', 'water_subtype', 'water_source', 'geometry']
  water_source counts: {'base.water': 7126, 'base.bathymetry': 1370}
  water_subtype nulls: 0
  water_id nulls: 0
  water_name nulls: 0
  station_name nulls: 0
  station_id nulls: 0
  station_id distinct: 34