xopr.opr_tools

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