1"""
  2Geospatial operations for radar datasets
  3
  4This module provides utilities for working with geospatial data in the context
  5of polar radar analysis. It includes functions for:
  6
  7- Loading and filtering Antarctic regional boundaries from the MEaSUREs dataset
  8- Projecting radar dataset coordinates between different coordinate reference systems
  9- Reprojecting GeoJSON geometries for spatial analysis
 10
 11The module handles common polar coordinate system transformations, particularly
 12between WGS84 (EPSG:4326) and Antarctic Polar Stereographic (EPSG:3031), and
 13provides tools for merging and simplifying regional boundaries.
 14
 15Functions
 16---------
 17get_antarctic_regions : Load and filter Antarctic regional boundaries
 18get_greenland_regions : Load and filter Greenland regional boundaries
 19project_dataset : Project dataset coordinates to a target CRS
 20project_geojson : Reproject geometries between coordinate systems
 21
 22See Also
 23--------
 24MEaSUREs Antarctic Boundaries: https://nsidc.org/data/nsidc-0709/versions/2
 25EPSG:3031 (Antarctic Polar Stereographic): https://epsg.io/3031
 26"""
 27
 28import geopandas as gpd
 29import shapely
 30import shapely.ops
 31from pyproj import Transformer
 32
 33
 34def get_antarctic_regions(
 35    name=None,
 36    region=None,
 37    subregion=None,
 38    type=None,
 39    merge_regions=True,
 40    merge_in_projection="EPSG:3031",
 41    simplify_tolerance=None,
 42    regions=None, # Deprecated alias for "region"
 43    subregions=None # Deprecated alias for "subregion"
 44):
 45    """
 46    Load and filter Antarctic regional boundaries from the MEaSUREs dataset.
 47
 48    The data product is derived from:
 49
 50    > Maps of Antarctic ice shelves, the Antarctic coastline and Antarctic basins. The maps are assembled from 2008-2009 ice-front data from ALOS PALSAR and ENVISAT ASAR data acquired during International Polar Year, 2007-2009 (IPY), the InSAR-based grounding line data (MEaSUREs Antarctic Grounding Line from Differential Satellite Radar Interferometry), augmented with other grounding line sources, the Antarctic ice velocity map (MEaSUREs InSAR-Based Antarctica Ice Velocity Map), and the Bedmap-2 DEM.
 51    >
 52    > This data set is part of the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) program.
 53    >
 54    > Mouginot, J., B. Scheuchl, and E. Rignot. 2017. MEaSUREs Antarctic Boundaries for IPY 2007-2009 from Satellite Radar, Version 2. [Indicate subset used]. Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center. doi: http://dx.doi.org/10.5067/AXE4121732AD. [Date Accessed].
 55    >
 56    > https://nsidc.org/data/nsidc-0709/versions/2
 57
 58    See these two notebooks for examples of how to use this function:
 59    - https://docs.englacial.org/xopr/search-and-scaling/
 60    - https://docs.englacial.org/xopr/crossovers/
 61
 62    Parameters
 63    ----------
 64    name : str or list, optional
 65        NAME field value(s) to filter by
 66    region : str or list, optional
 67        REGION field value(s) to filter by
 68    subregion : str or list, optional
 69        SUBREGION field value(s) to filter by
 70    type : str or list, optional
 71        TYPE field value(s) to filter by
 72    merge_regions : bool, default True
 73        If True, return a single merged geometry; if False, return list of geometries
 74    measures_boundaries_url : str, default "https://storage.googleapis.com/opr_stac/reference_geometry/measures_boundaries_4326.geojson"
 75        URL to the GeoJSON file containing Antarctic region boundaries
 76
 77    Returns
 78    -------
 79    list or dict
 80        If merge_regions=False: List of GeoJSON geometry dicts
 81        If merge_regions=True: Single GeoJSON geometry dict of merged regions
 82
 83    Examples
 84    --------
 85    # Get George VI ice shelf
 86    >>> george_vi = get_antarctic_regions(name="George_VI", type="FL")
 87
 88    # Get all ice shelves, merged into one geometry
 89    >>> all_shelves = get_antarctic_regions(type="FL", merge_regions=True)
 90
 91    # Get multiple regions by name
 92    >>> regions = get_antarctic_regions(name=["George_VI", "LarsenC"])
 93    """
 94
 95    if regions is not None:
 96        print("Warning: 'regions' parameter is deprecated; use 'region' instead.")
 97        if region is not None:
 98            raise ValueError("Cannot use both 'regions' and 'region' parameters.")
 99        region = regions
100    if subregions is not None:
101        print("Warning: 'subregions' parameter is deprecated; use 'subregion' instead.")
102        if subregion is not None:
103            raise ValueError("Cannot use both 'subregions' and 'subregion' parameters.")
104        subregion = subregions
105
106    return _get_regions(
107        name=name,
108        region=region,
109        subregion=subregion,
110        type=type,
111        merge_regions=merge_regions,
112        regions_geojson_url="https://storage.googleapis.com/opr_stac/reference_geometry/measures_boundaries_4326.geojson",
113        merge_in_projection=merge_in_projection,
114        simplify_tolerance=simplify_tolerance
115    )
116
117def get_greenland_regions(
118    name=None,
119    region=None,
120    subregion=None,
121    type=None,
122    merge_regions=True,
123    merge_in_projection="EPSG:3413",
124    simplify_tolerance=None,
125):
126    """
127    Load and filter Greenland regional boundaries as defined by the GrIMP project
128
129    The data product is derived from:
130
131    > We divide Greenland, including its peripheral glaciers and ice caps, into 260 basins grouped in seven regions:  southwest (SW), central west (CW), (iii) northwest (NW), north (NO), northeast (NE), central east (CE), and southeast (SE). These regions are selected based on ice flow regimes, climate, and the need to partition the ice sheet into zones comparable in size (200,000 km2 to 400,000 km2) and ice production (50 Gt/y to 100 Gt/y, or billion tons per year). Out of the 260 surveyed glaciers, 217 are marine-terminating, i.e., calving into icebergs and melting in contact with ocean waters, and 43 are land-terminating.The actual number of land-terminating glaciers is far larger than 43, but we lump them into larger units for simplification.
132    >
133    > Each glacier catchment is defined using a combination of ice flow direction and surface slope. In areas of fast flow (> 100 m), we use a composite velocity mosaic (Mouginot et al. 2017). In slowmoving areas, we use surface slope using the GIMP DEM (https://nsidc.org/data/nsidc- 0715/versions/1) (Howat et al. 2014) smoothed over 10 ice thicknesses to remove shortwavelength undulations.
134    >
135    > References:
136    >
137    > Mouginot J, Rignot E, Scheuchl B, Millan R (2017) Comprehensive annual ice sheet velocity mapping using landsat-8, sentinel-1, and radarsat-2 data. Remote Sensing 9(4).
138    >
139    > Howat IM, Negrete A, Smith BE (2014) The greenland ice mapping project (GIMP) land classification and surface elevation data sets. The Cryosphere 8(4):1509–1518.
140    >
141    > https://datadryad.org/dataset/doi:10.7280/D1WT11
142
143    Parameters
144    ----------
145    name : str or list, optional
146        NAME field value(s) to filter by
147    region : str or list, optional
148        REGION field value(s) to filter by
149    subregion : str or list, optional
150        SUBREGION field value(s) to filter by
151    type : str or list, optional
152        TYPE field value(s) to filter by
153    merge_regions : bool, default True
154        If True, return a single merged geometry; if False, return list of geometries
155    measures_boundaries_url : str, default "https://storage.googleapis.com/opr_stac/reference_geometry/measures_boundaries_4326.geojson"
156        URL to the GeoJSON file containing Antarctic region boundaries
157
158    Returns
159    -------
160    list or dict
161        If merge_regions=False: List of GeoJSON geometry dicts
162        If merge_regions=True: Single GeoJSON geometry dict of merged regions
163    """
164
165    greenland_boundaries_url = "https://storage.googleapis.com/opr_stac/reference_geometry/greenland_boundaries_4326.geojson"
166
167    return _get_regions(
168        name=name,
169        region=region,
170        subregion=subregion,
171        type=type,
172        merge_regions=merge_regions,
173        regions_geojson_url=greenland_boundaries_url,
174        merge_in_projection=merge_in_projection,
175        simplify_tolerance=simplify_tolerance
176    )
177
178
179def _get_regions(
180    name=None,
181    region=None,
182    subregion=None,
183    type=None,
184    merge_regions=True,
185    regions_geojson_url = None,
186    merge_in_projection=None,
187    simplify_tolerance=None
188):
189    if regions_geojson_url is None:
190        raise ValueError("regions_geojson_url must be provided")
191
192    if merge_in_projection is None:
193        raise ValueError("merge_in_projection must be provided")
194
195    # Load the boundaries GeoJSON from the reference URL
196    filtered = gpd.read_file(regions_geojson_url)
197
198    standard_columns = {
199        'NAME': ['NAME', 'NAMES'],
200        'REGION': ['REGION', 'Regions'],
201        'SUBREGION': ['SUBREGION', 'Subregions'],
202        'TYPE': ['TYPE', 'TYPES', 'GL_TYPE']
203    }
204
205    # Standardize column names
206    for standard_name, possible_names in standard_columns.items():
207        for col in possible_names:
208            if col in filtered.columns:
209                filtered = filtered.rename(columns={col: standard_name})
210                break
211
212    # Apply filters based on provided parameters
213    if name is not None:
214        if isinstance(name, str):
215            name = [name]
216        filtered = filtered[filtered['NAME'].isin(name)]
217
218    if region is not None:
219        if isinstance(region, str):
220            region = [region]
221        filtered = filtered[filtered['REGION'].isin(region)]
222
223    if subregion is not None:
224        if isinstance(subregion, str):
225            subregion = [subregion]
226        filtered = filtered[filtered['SUBREGION'].isin(subregion)]
227
228    if type is not None:
229        if isinstance(type, str):
230            type = [type]
231        filtered = filtered[filtered['TYPE'].isin(type)]
232
233    if len(filtered) == 0:
234        return [] if not merge_regions else None
235
236    if merge_regions:
237
238        if merge_in_projection:
239            filtered = filtered.to_crs(merge_in_projection)
240
241        # Check for invalid regions and attempt to fix them
242        invalid_geometries = filtered[~filtered.is_valid]
243        if len(invalid_geometries) > 0:
244            filtered = filtered.make_valid()
245            print(f"Warning: {len(invalid_geometries)} invalid geometries were fixed before merging.")
246            if merge_in_projection != "EPSG:3031":
247                print("Consider using merge_in_projection='EPSG:3031' to reproject before merging.")
248            print(f"Invalid geometry regions were: {', '.join(invalid_geometries['NAME'])}")
249
250        merged = shapely.ops.unary_union(filtered.geometry)  # geopandas < 1.0 compat
251
252        if simplify_tolerance is None and (merge_in_projection is not None): # Set a reasonable default based on the size
253            area_km2 = merged.area / 1e6
254            if area_km2 < 10000:
255                simplify_tolerance = 0
256            elif area_km2 < 100000:
257                print(f"Area is {area_km2:.1f} km^2, automatically applying 100m simplification tolerance")
258                print(f"To disable simplification, set simplify_tolerance=0")
259                simplify_tolerance = 100
260            else:
261                print(f"Area is {area_km2:.1f} km^2, automatically applying 1km simplification tolerance")
262                print(f"To disable simplification, set simplify_tolerance=0")
263                simplify_tolerance = 1000
264
265        if simplify_tolerance and (simplify_tolerance > 0):
266            merged = shapely.buffer(merged, simplify_tolerance).simplify(tolerance=simplify_tolerance)
267
268        if merge_in_projection:
269            merged = project_geojson(merged, source_crs=merge_in_projection, target_crs="EPSG:4326")
270
271        return merged
272    else:
273        return filtered
274
275def project_dataset(ds, target_crs):
276    """
277    Project dataset coordinates from WGS84 to a target coordinate reference system.
278
279    Takes longitude and latitude coordinates from a dataset and projects them to
280    the specified target CRS, adding 'x' and 'y' coordinate arrays to the dataset.
281
282    Parameters
283    ----------
284    ds : xarray.Dataset
285        Input dataset containing 'Longitude' and 'Latitude' coordinates
286    target_crs : cartopy.crs.CRS or str
287        Target coordinate reference system. Can be a cartopy CRS object or
288        a string representation (e.g., "EPSG:3031")
289
290    Returns
291    -------
292    xarray.Dataset
293        Dataset with added 'x' and 'y' coordinate arrays in the target CRS
294
295    Examples
296    --------
297    >>> import cartopy.crs as ccrs
298    >>> projected_ds = project_dataset(ds, ccrs.SouthPolarStereo(true_scale_latitude=-71))
299    >>> projected_ds = project_dataset(ds, "EPSG:3031")
300    """
301    if hasattr(target_crs, 'to_epsg') and target_crs.to_epsg():
302        target_crs_str = f"EPSG:{target_crs.to_epsg()}"
303    elif isinstance(target_crs, str):
304        target_crs_str = target_crs
305    else:
306        target_crs_str = target_crs.to_proj4_string()
307
308    transformer = Transformer.from_crs("EPSG:4326", target_crs_str, always_xy=True)
309    projected_coords = transformer.transform(ds['Longitude'].values, ds['Latitude'].values)
310
311    ds = ds.assign_coords({
312        'x': (('slow_time'), projected_coords[0]),
313        'y': (('slow_time'), projected_coords[1])
314    })
315    return ds
316
317def project_geojson(geometry, source_crs="EPSG:4326", target_crs=None):
318    """
319    Project a geometry from one coordinate reference system to another.
320
321    Uses pyproj.Transformer to reproject geometries between different
322    coordinate reference systems. Commonly used for projecting geometries
323    from WGS84 (lat/lon) to polar stereographic projections.
324
325    Parameters
326    ----------
327    geometry : shapely.geometry.base.BaseGeometry
328        Input geometry to be projected
329    source_crs : str, default "EPSG:4326"
330        Source coordinate reference system (default is WGS84)
331    target_crs : str, default None
332        Target coordinate reference system. If None, automatically selects
333        Antarctic Polar Stereographic (EPSG:3031) for southern hemisphere
334        and Arctic Polar Stereographic (EPSG:3413) for northern hemisphere based on
335        the geometry's centroid latitude.
336
337    Returns
338    -------
339    shapely.geometry.base.BaseGeometry
340        Projected geometry in the target coordinate reference system
341
342    Examples
343    --------
344    >>> from shapely.geometry import Point
345    >>> point = Point(-70, -75)  # lon, lat in WGS84
346    >>> projected = project_geojson(point, "EPSG:4326", "EPSG:3031")
347    """
348
349    if target_crs is None:
350        if geometry.centroid.y < 0:
351            target_crs = "EPSG:3031"  # Antarctic Polar Stereographic
352        else:
353            target_crs = "EPSG:3413"  # Arctic Polar Stereographic
354
355    transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)
356    projected_geometry = shapely.ops.transform(transformer.transform, geometry)
357    return projected_geometry