xopr.opr_tools

  1"""
  2Helper functions for working with OPR radar datasets.
  3
  4This module provides utility functions for common operations on radar datasets,
  5including merging frames from the same flight segment, generating citations, and
  6finding crossover points between flight lines.
  7"""
  8
  9import warnings
 10from typing import Iterable, Union
 11
 12import geopandas as gpd
 13import xarray as xr
 14
 15from xopr.util import get_ror_display_name, merge_dicts_no_conflicts
 16
 17
 18def merge_frames(frames: Iterable[xr.Dataset]) -> Union[list[xr.Dataset], xr.Dataset]:
 19    """
 20    Merge radar frames by segment.
 21
 22    Frames from the same segment are concatenated along slow_time and sorted.
 23
 24    Parameters
 25    ----------
 26    frames : Iterable[xr.Dataset]
 27        Radar frame Datasets to merge.
 28
 29    Returns
 30    -------
 31    xr.Dataset or list[xr.Dataset]
 32        Single merged Dataset if all frames are from one segment, otherwise
 33        list of Datasets (one per segment).
 34    """
 35    segments = {}
 36
 37    mergable_keys = ['source_url', 'collection', 'data_product', 'granule']
 38
 39    for frame in frames:
 40        # Get segment path from frame attributes
 41        granule = frame.attrs.get('granule')
 42        date, segment_id, frame_id = granule.split('_')
 43        segment_path = f"{date}_{segment_id}"
 44
 45        if not segment_path:
 46            warnings.warn("Frame missing 'granule' attribute or it was not in the expected format, skipping.", UserWarning)
 47            continue
 48
 49        if segment_path not in segments:
 50            segments[segment_path] = []
 51
 52        segments[segment_path].append(frame)
 53
 54    # Merge frames for each segment
 55    merged_segments = []
 56    for segment_path, segment_frames in segments.items():
 57        merged_segment = xr.concat(segment_frames, dim='slow_time', combine_attrs=merge_dicts_no_conflicts).sortby('slow_time')
 58
 59        for k in mergable_keys:
 60            if k not in merged_segment.attrs:
 61                merged_segment.attrs[k] = set([v for v in {f.attrs.get(k) for f in segment_frames} if v is not None])
 62
 63        merged_segments.append(merged_segment)
 64
 65    if len(merged_segments) == 1:
 66        return merged_segments[0]
 67    else:
 68        return merged_segments
 69
 70def generate_citation(ds : xr.Dataset) -> str:
 71        """
 72        Generate citation and acknowledgment text for dataset.
 73
 74        Parameters
 75        ----------
 76        ds : xr.Dataset
 77            Dataset with doi, ror, and funder_text attributes.
 78
 79        Returns
 80        -------
 81        str
 82            Formatted citation string including data citation, DOI, funding
 83            acknowledgments, and OPR toolbox citation.
 84        """
 85
 86        citation_string = ""
 87        any_citation_info = False
 88
 89        citation_string += "== Data Citation ==\n"
 90
 91        if 'ror' in ds.attrs and ds.attrs['ror']:
 92            any_citation_info = True
 93            if isinstance(ds.attrs['ror'], (set, list)):
 94                institution_name = ', '.join([get_ror_display_name(ror) for ror in ds.attrs['ror']])
 95            else:
 96                institution_name = get_ror_display_name(ds.attrs['ror'])
 97
 98            citation_string += f"This data was collected by {institution_name}.\n"
 99
100        if 'doi' in ds.attrs and ds.attrs['doi']:
101            any_citation_info = True
102            citation_string += f"Please cite the dataset DOI: https://doi.org/{ds.attrs['doi']}\n"
103
104        if 'funder_text' in ds.attrs and ds.attrs['funder_text']:
105            any_citation_info = True
106            citation_string += f"Please include the following funder acknowledgment:\n{ds.attrs['funder_text']}\n"
107
108        if not any_citation_info:
109            citation_string += "No specific citation information was retrieved for this dataset. By default, please cite:\n"
110            citation_string += "CReSIS. 2024. REPLACE_WITH_RADAR_NAME Data, Lawrence, Kansas, USA. Digital Media. http://data.cresis.ku.edu/."
111
112        # Add general OPR Toolbox citation
113        citation_string += "\n== Processing Citation ==\n"
114        citation_string += "Data was processed using the Open Polar Radar (OPR) Toolbox: https://doi.org/10.5281/zenodo.5683959\n"
115        citation_string += "Please cite the OPR Toolbox as:\n"
116        citation_string += "Open Polar Radar. (2024). opr (Version 3.0.1) [Computer software]. https://gitlab.com/openpolarradar/opr/. https://doi.org/10.5281/zenodo.5683959\n"
117        citation_string += "And include the following acknowledgment:\n"
118        citation_string += "We acknowledge the use of software from Open Polar Radar generated with support from the University of Kansas, NASA grants 80NSSC20K1242 and 80NSSC21K0753, and NSF grants OPP-2027615, OPP-2019719, OPP-1739003, IIS-1838230, RISE-2126503, RISE-2127606, and RISE-2126468.\n"
119
120        return citation_string if citation_string else "No citation information available for this dataset."
121
122
123def find_intersections(gdf : gpd.GeoDataFrame, remove_self_intersections: bool = True,
124                        remove_adjacent_intersections: bool = True,
125                        calculate_crossing_angles: bool = False) -> gpd.GeoDataFrame:
126    """
127    Find crossover points between radar flight lines.
128
129    Parameters
130    ----------
131    gdf : gpd.GeoDataFrame
132        GeoDataFrame of STAC items with flight line geometries.
133    remove_self_intersections : bool, optional
134        If True, exclude geometry intersecting with itself.
135    remove_adjacent_intersections : bool, optional
136        If True, exclude consecutive frames in same flight line.
137    calculate_crossing_angles : bool, optional
138        If True, calculate and include crossing angles in degrees.
139
140    Returns
141    -------
142    gpd.GeoDataFrame
143        GeoDataFrame with columns for intersecting pairs (id_1, id_2),
144        intersection_geometry, and optionally crossing_angle.
145    """
146
147    tmp_df = gdf.reset_index()
148    tmp_df['geom'] = tmp_df.geometry
149    intersections = gpd.sjoin(tmp_df, tmp_df, how='inner', predicate='intersects', lsuffix='1', rsuffix='2')
150
151    if remove_self_intersections:
152        intersections = intersections[intersections['id_1'] != intersections['id_2']]
153
154    intersections['intersection_geometry'] = intersections.apply(lambda row: row['geom_1'].intersection(row['geom_2']), axis=1)
155    intersections.set_geometry('intersection_geometry', inplace=True, crs=gdf.crs)
156    intersections = intersections.drop_duplicates(subset=['intersection_geometry'])
157    intersections = intersections.explode(index_parts=True).reset_index(drop=True)
158
159    intersections_tmp = intersections[['id_1', 'id_2', 'intersection_geometry', 'collection_1', 'collection_2', 'geom_1', 'geom_2']].copy()
160
161    for k in ['opr:date', 'opr:segment', 'opr:frame']:
162        intersections_tmp[f'{k}_1'] = intersections['properties_1'].apply(lambda x: x[k])
163        intersections_tmp[f'{k}_2'] = intersections['properties_2'].apply(lambda x: x[k])
164
165    intersections = intersections_tmp
166
167    if remove_adjacent_intersections:
168        intersections = intersections[
169            (intersections['opr:date_1'] != intersections['opr:date_2']) |
170            (intersections['opr:segment_1'] != intersections['opr:segment_2']) |
171            ((intersections['opr:frame_1'] != (intersections['opr:frame_2'] + 1)) &
172            (intersections['opr:frame_1'] != (intersections['opr:frame_2'] - 1)))
173        ]
174
175    if calculate_crossing_angles:
176        intersections['crossing_angle'] = intersections.apply(
177            lambda row: _calculate_crossing_angle(row['geom_1'], row['geom_2'], row['intersection_geometry']),
178            axis=1
179        )
180
181    return intersections
182
183def _calculate_crossing_angle(line1, line2, intersection_point):
184    """
185    Calculate the crossing angle between two lines at their intersection point.
186
187    Parameters
188    ----------
189    line1 : shapely.geometry.LineString
190        The first line.
191    line2 : shapely.geometry.LineString
192        The second line.
193    intersection_point : shapely.geometry.Point
194        The point of intersection between the two lines.
195
196    Returns
197    -------
198    float
199        The crossing angle in degrees.
200    """
201    import numpy as np
202
203    def get_line_angle(line, point):
204        # Get the nearest point on the line to the intersection point
205        nearest_point = line.interpolate(line.project(point))
206        # Get a small segment of the line around the nearest point to calculate the angle
207        buffer_distance = 1e-6  # Small distance to create a segment
208        start_point = line.interpolate(max(0, line.project(nearest_point) - buffer_distance))
209        end_point = line.interpolate(min(line.length, line.project(nearest_point) + buffer_distance))
210
211        # Calculate the angle of the segment
212        delta_x = end_point.x - start_point.x
213        delta_y = end_point.y - start_point.y
214        angle = np.arctan2(delta_y, delta_x)
215        return angle
216
217    angle1 = get_line_angle(line1, intersection_point)
218    angle2 = get_line_angle(line2, intersection_point)
219
220    # Calculate the absolute difference in angles and convert to degrees
221    angle_diff = np.abs(angle1 - angle2)
222    crossing_angle = np.degrees(angle_diff)
223
224    # Ensure the angle is between 0 and 180 degrees
225    if crossing_angle > 180:
226        crossing_angle = 360 - crossing_angle
227
228    if crossing_angle > 90:
229        crossing_angle = 180 - crossing_angle
230
231    return crossing_angle