|
| 1 | +""" |
| 2 | +Diagnostics module for xarray-spatial. |
| 3 | +
|
| 4 | +Provides utilities to help users identify common pitfalls and issues |
| 5 | +with DataArrays before running xarray-spatial operations. |
| 6 | +""" |
| 7 | +from __future__ import annotations |
| 8 | + |
| 9 | +from dataclasses import dataclass, field |
| 10 | +from typing import List, Optional |
| 11 | + |
| 12 | +import xarray as xr |
| 13 | + |
| 14 | +from xrspatial.utils import ( |
| 15 | + _infer_coord_unit_type, |
| 16 | + _infer_vertical_unit_type, |
| 17 | + get_dataarray_resolution, |
| 18 | +) |
| 19 | + |
| 20 | + |
| 21 | +@dataclass |
| 22 | +class DiagnosticIssue: |
| 23 | + """Represents a single diagnostic issue found during analysis.""" |
| 24 | + code: str |
| 25 | + severity: str # 'warning' or 'error' |
| 26 | + message: str |
| 27 | + suggestion: str |
| 28 | + |
| 29 | + |
| 30 | +@dataclass |
| 31 | +class DiagnosticReport: |
| 32 | + """Results from diagnosing a DataArray.""" |
| 33 | + issues: List[DiagnosticIssue] = field(default_factory=list) |
| 34 | + horizontal_unit_type: Optional[str] = None |
| 35 | + vertical_unit_type: Optional[str] = None |
| 36 | + resolution: Optional[tuple] = None |
| 37 | + |
| 38 | + @property |
| 39 | + def has_issues(self) -> bool: |
| 40 | + return len(self.issues) > 0 |
| 41 | + |
| 42 | + @property |
| 43 | + def has_warnings(self) -> bool: |
| 44 | + return any(i.severity == 'warning' for i in self.issues) |
| 45 | + |
| 46 | + @property |
| 47 | + def has_errors(self) -> bool: |
| 48 | + return any(i.severity == 'error' for i in self.issues) |
| 49 | + |
| 50 | + def __str__(self) -> str: |
| 51 | + if not self.issues: |
| 52 | + return "No issues detected." |
| 53 | + |
| 54 | + lines = [] |
| 55 | + for issue in self.issues: |
| 56 | + lines.append(f"[{issue.severity.upper()}] {issue.code}: {issue.message}") |
| 57 | + lines.append(f" Suggestion: {issue.suggestion}") |
| 58 | + return "\n".join(lines) |
| 59 | + |
| 60 | + |
| 61 | +def _check_unit_mismatch(agg: xr.DataArray, report: DiagnosticReport) -> None: |
| 62 | + """ |
| 63 | + Check for horizontal vs vertical unit mismatch. |
| 64 | +
|
| 65 | + Detects the common case of coordinates in degrees (lon/lat) with |
| 66 | + elevation values in meters/feet. |
| 67 | + """ |
| 68 | + try: |
| 69 | + cellsize_x, cellsize_y = get_dataarray_resolution(agg) |
| 70 | + report.resolution = (cellsize_x, cellsize_y) |
| 71 | + except Exception: |
| 72 | + return |
| 73 | + |
| 74 | + if len(agg.dims) < 2: |
| 75 | + return |
| 76 | + |
| 77 | + dim_y, dim_x = agg.dims[-2], agg.dims[-1] |
| 78 | + coord_x = agg.coords.get(dim_x, None) |
| 79 | + coord_y = agg.coords.get(dim_y, None) |
| 80 | + |
| 81 | + if coord_x is None or coord_y is None: |
| 82 | + return |
| 83 | + |
| 84 | + horiz_x = _infer_coord_unit_type(coord_x, cellsize_x) |
| 85 | + horiz_y = _infer_coord_unit_type(coord_y, cellsize_y) |
| 86 | + vert = _infer_vertical_unit_type(agg) |
| 87 | + |
| 88 | + report.vertical_unit_type = vert |
| 89 | + |
| 90 | + horiz_types = {horiz_x, horiz_y} - {"unknown"} |
| 91 | + if horiz_types: |
| 92 | + report.horizontal_unit_type = next(iter(horiz_types)) |
| 93 | + |
| 94 | + if not horiz_types or vert == "unknown": |
| 95 | + return |
| 96 | + |
| 97 | + if "degrees" in horiz_types and vert == "elevation": |
| 98 | + report.issues.append(DiagnosticIssue( |
| 99 | + code="UNIT_MISMATCH", |
| 100 | + severity="warning", |
| 101 | + message=( |
| 102 | + "Input DataArray appears to have coordinates in degrees " |
| 103 | + "but elevation values in a linear unit (e.g. meters/feet)." |
| 104 | + ), |
| 105 | + suggestion=( |
| 106 | + "Slope/aspect/curvature operations expect horizontal distances " |
| 107 | + "in the same units as vertical. Consider reprojecting to a " |
| 108 | + "projected CRS with meter-based coordinates." |
| 109 | + ), |
| 110 | + )) |
| 111 | + |
| 112 | + |
| 113 | +def diagnose(agg: xr.DataArray, tool: Optional[str] = None) -> DiagnosticReport: |
| 114 | + """ |
| 115 | + Diagnose a DataArray for common xarray-spatial pitfalls. |
| 116 | +
|
| 117 | + Runs a series of heuristic checks to identify potential issues |
| 118 | + that could lead to incorrect results when using xarray-spatial |
| 119 | + functions. |
| 120 | +
|
| 121 | + Parameters |
| 122 | + ---------- |
| 123 | + agg : xr.DataArray |
| 124 | + The input DataArray to diagnose. |
| 125 | + tool : str, optional |
| 126 | + Name of the xarray-spatial tool you intend to use (e.g., 'slope', |
| 127 | + 'aspect', 'curvature'). When specified, only diagnostics relevant |
| 128 | + to that tool are run. If None, all diagnostics are run. |
| 129 | +
|
| 130 | + Returns |
| 131 | + ------- |
| 132 | + DiagnosticReport |
| 133 | + A report containing any issues found, along with inferred |
| 134 | + metadata about the DataArray. |
| 135 | +
|
| 136 | + Examples |
| 137 | + -------- |
| 138 | + >>> import xarray as xr |
| 139 | + >>> import numpy as np |
| 140 | + >>> from xrspatial.diagnostics import diagnose |
| 141 | + >>> # Create a DataArray with lon/lat coordinates but meter elevations |
| 142 | + >>> data = np.random.rand(100, 100) * 1000 + 500 |
| 143 | + >>> da = xr.DataArray( |
| 144 | + ... data, |
| 145 | + ... dims=['y', 'x'], |
| 146 | + ... coords={ |
| 147 | + ... 'y': np.linspace(40.0, 41.0, 100), |
| 148 | + ... 'x': np.linspace(-105.0, -104.0, 100), |
| 149 | + ... } |
| 150 | + ... ) |
| 151 | + >>> report = diagnose(da) |
| 152 | + >>> print(report) |
| 153 | + [WARNING] UNIT_MISMATCH: Input DataArray appears to have coordinates... |
| 154 | + >>> if report.has_warnings: |
| 155 | + ... print("Consider reprojecting your data!") |
| 156 | + """ |
| 157 | + report = DiagnosticReport() |
| 158 | + |
| 159 | + # Tools that are sensitive to unit mismatch |
| 160 | + unit_mismatch_tools = {'slope', 'aspect', 'curvature', 'hillshade'} |
| 161 | + |
| 162 | + # Run unit mismatch check if tool is relevant or no tool specified |
| 163 | + if tool is None or tool.lower() in unit_mismatch_tools: |
| 164 | + _check_unit_mismatch(agg, report) |
| 165 | + |
| 166 | + return report |
0 commit comments