1import geopandas as gpd
  2import shapely
  3import shapely.ops
  4from pyproj import Transformer
  5
  6
  7def get_antarctic_regions(
  8    name=None,
  9    regions=None,
 10    subregions=None,
 11    type=None,
 12    merge_regions=True,
 13    measures_boundaries_url = "https://storage.googleapis.com/opr_stac/reference_geometry/measures_boundaries_4326.geojson",
 14    merge_in_projection="EPSG:3031",
 15    simplify_tolerance=None
 16):
 17    """
 18    Load and filter Antarctic regional boundaries from the MEASURES dataset.
 19
 20    Parameters
 21    ----------
 22    name : str or list, optional
 23        NAME field value(s) to filter by
 24    regions : str or list, optional
 25        REGIONS field value(s) to filter by
 26    subregions : str or list, optional
 27        SUBREGION field value(s) to filter by
 28    type : str or list, optional
 29        TYPE field value(s) to filter by
 30    merge_regions : bool, default True
 31        If True, return a single merged geometry; if False, return list of geometries
 32    measures_boundaries_url : str, default "https://storage.googleapis.com/opr_stac/reference_geometry/measures_boundaries_4326.geojson"
 33        URL to the GeoJSON file containing Antarctic region boundaries
 34
 35    Returns
 36    -------
 37    list or dict
 38        If merge_regions=False: List of GeoJSON geometry dicts
 39        If merge_regions=True: Single GeoJSON geometry dict of merged regions
 40
 41    Examples
 42    --------
 43    # Get George VI ice shelf
 44    >>> george_vi = get_antarctic_regions(name="George_VI", type="FL")
 45
 46    # Get all ice shelves, merged into one geometry
 47    >>> all_shelves = get_antarctic_regions(type="FL", merge_regions=True)
 48
 49    # Get multiple regions by name
 50    >>> regions = get_antarctic_regions(name=["George_VI", "LarsenC"])
 51    """
 52
 53
 54    # Load the boundaries GeoJSON from the reference URL
 55    filtered = gpd.read_file(measures_boundaries_url)
 56
 57    # Apply filters based on provided parameters
 58    if name is not None:
 59        if isinstance(name, str):
 60            name = [name]
 61        filtered = filtered[filtered['NAME'].isin(name)]
 62
 63    if regions is not None:
 64        if isinstance(regions, str):
 65            regions = [regions]
 66        filtered = filtered[filtered['Regions'].isin(regions)]
 67
 68    if subregions is not None:
 69        if isinstance(subregions, str):
 70            subregions = [subregions]
 71        filtered = filtered[filtered['Subregions'].isin(subregions)]
 72
 73    if type is not None:
 74        if isinstance(type, str):
 75            type = [type]
 76        filtered = filtered[filtered['TYPE'].isin(type)]
 77
 78    if len(filtered) == 0:
 79        return [] if not merge_regions else None
 80
 81    if merge_regions:
 82
 83        if merge_in_projection:
 84            filtered = filtered.to_crs(merge_in_projection)
 85
 86        # Check for invalid regions and attempt to fix them
 87        invalid_geometries = filtered[~filtered.is_valid]
 88        if len(invalid_geometries) > 0:
 89            filtered = filtered.make_valid()
 90            print(f"Warning: {len(invalid_geometries)} invalid geometries were fixed before merging.")
 91            if merge_in_projection != "EPSG:3031":
 92                print("Consider using merge_in_projection='EPSG:3031' to reproject before merging.")
 93            print(f"Invalid geometry regions were: {', '.join(invalid_geometries['NAME'])}")
 94
 95        merged = shapely.ops.unary_union(filtered.geometry)  # geopandas < 1.0 compat
 96
 97        if simplify_tolerance is None and (merge_in_projection == "EPSG:3031"): # Set a reasonable default based on the size
 98            area_km2 = merged.area / 1e6
 99            if area_km2 < 10000:
100                simplify_tolerance = 0
101            elif area_km2 < 100000:
102                print(f"Area is {area_km2:.1f} km^2, automatically applying 100m simplification tolerance")
103                print(f"To disable simplification, set simplify_tolerance=0")
104                simplify_tolerance = 100
105            else:
106                print(f"Area is {area_km2:.1f} km^2, automatically applying 1km simplification tolerance")
107                print(f"To disable simplification, set simplify_tolerance=0")
108                simplify_tolerance = 1000
109
110        if simplify_tolerance and (simplify_tolerance > 0):
111            merged = shapely.buffer(merged, simplify_tolerance).simplify(tolerance=simplify_tolerance)
112
113        if merge_in_projection:
114            merged = project_geojson(merged, source_crs=merge_in_projection, target_crs="EPSG:4326")
115
116        return merged
117    else:
118        return filtered
119
120def project_dataset(ds, target_crs):
121    """
122    Project dataset coordinates from WGS84 to a target coordinate reference system.
123
124    Takes longitude and latitude coordinates from a dataset and projects them to
125    the specified target CRS, adding 'x' and 'y' coordinate arrays to the dataset.
126
127    Parameters
128    ----------
129    ds : xarray.Dataset
130        Input dataset containing 'Longitude' and 'Latitude' coordinates
131    target_crs : cartopy.crs.CRS or str
132        Target coordinate reference system. Can be a cartopy CRS object or
133        a string representation (e.g., "EPSG:3031")
134
135    Returns
136    -------
137    xarray.Dataset
138        Dataset with added 'x' and 'y' coordinate arrays in the target CRS
139
140    Examples
141    --------
142    >>> import cartopy.crs as ccrs
143    >>> projected_ds = project_dataset(ds, ccrs.SouthPolarStereo())
144    >>> projected_ds = project_dataset(ds, "EPSG:3031")
145    """
146    if hasattr(target_crs, 'to_epsg') and target_crs.to_epsg():
147        target_crs_str = f"EPSG:{target_crs.to_epsg()}"
148    elif isinstance(target_crs, str):
149        target_crs_str = target_crs
150    else:
151        target_crs_str = target_crs.to_proj4_string()
152
153    transformer = Transformer.from_crs("EPSG:4326", target_crs_str, always_xy=True)
154    projected_coords = transformer.transform(ds['Longitude'].values, ds['Latitude'].values)
155
156    ds = ds.assign_coords({
157        'x': (('slow_time'), projected_coords[0]),
158        'y': (('slow_time'), projected_coords[1])
159    })
160    return ds
161
162def project_geojson(geometry, source_crs="EPSG:4326", target_crs="EPSG:3031"):
163    """
164    Project a geometry from one coordinate reference system to another.
165
166    Uses pyproj.Transformer to reproject geometries between different
167    coordinate reference systems. Commonly used for projecting geometries
168    from WGS84 (lat/lon) to polar stereographic projections.
169
170    Parameters
171    ----------
172    geometry : shapely.geometry.base.BaseGeometry
173        Input geometry to be projected
174    source_crs : str, default "EPSG:4326"
175        Source coordinate reference system (default is WGS84)
176    target_crs : str, default "EPSG:3031"
177        Target coordinate reference system (default is Antarctic Polar Stereographic)
178
179    Returns
180    -------
181    shapely.geometry.base.BaseGeometry
182        Projected geometry in the target coordinate reference system
183
184    Examples
185    --------
186    >>> from shapely.geometry import Point
187    >>> point = Point(-70, -75)  # lon, lat in WGS84
188    >>> projected = project_geojson(point, "EPSG:4326", "EPSG:3031")
189    """
190    transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)
191    projected_geometry = shapely.ops.transform(transformer.transform, geometry)
192    return projected_geometry