305eb64a-9a20-4e4e-9c1e-8cc67e861cd6

Committed 09 Sep 2025 04:56PM UTC coverage: 51.086% (-43.4%) from 94.444%

Build # 305eb64a-9a20-4e4e-9c1e-8cc67e861cd6

Build Type

push

circleci

Committed by

web-flow

Commit Message

Merge pull request #766 from dbekaert/dev

v0.6.0

Run Details

133 of 261 new or added lines in 18 files covered. (50.96%)

3199 of 6262 relevant lines covered (51.09%)

0.51 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

77.49

/tools/RAiDER/llreader.py

# noqa: D100
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# Author: Jeremy Maurer, Raymond Hogenson & David Bekaert
# Copyright 2019, by the California Institute of Technology. ALL RIGHTS
# RESERVED. United States Government Sponsorship acknowledged.
#
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
import os
from pathlib import Path
from typing import Optional, Union

import numpy as np
import pyproj
import xarray as xr


try:
    import pandas as pd
except ImportError:
    pd = None

from pyproj import CRS

from RAiDER.logger import logger
from RAiDER.types import BB, RIO


class AOI:
    """
    This instantiates a generic AOI class object.

    Attributes:
       _bounding_box    - S N W E bounding box
       _proj            - pyproj-compatible CRS
       _type            - Type of AOI
    """

    def __init__(self, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:
        self._output_directory = output_directory
        self._bounding_box = None
        self._proj = CRS.from_epsg(4326)
        self._geotransform = None
        self._cube_spacing_m = cube_spacing_in_m
    

    def __repr__(self):
        return f'AOI: {self.__class__.__name__}({self._bounding_box}, {self._type})'

    def type(self):
        return self._type

    def bounds(self):
        return list(self._bounding_box).copy()

    def geotransform(self):
        return self._geotransform

    def projection(self):
        return self._proj

    def get_output_spacing(self, crs=4326):
        """Return the output spacing in desired units."""
        output_spacing_deg = self._output_spacing
        if not isinstance(crs, CRS):
            crs = CRS.from_epsg(crs)

        ## convert it to meters users wants a projected coordinate system
        if all(axis_info.unit_name == 'degree' for axis_info in crs.axis_info):
            output_spacing = output_spacing_deg
        else:
            output_spacing = output_spacing_deg * 1e5

        return output_spacing

    def set_output_spacing(self, ll_res=None) -> None:
        """Calculate the spacing for the output grid and weather model.

        Use the requested spacing if exists or the weather model grid itself

        Returns:
            None. Sets self._output_spacing
        """
        assert ll_res or self._cube_spacing_m, 'Must pass lat/lon resolution if _cube_spacing_m is None'

        out_spacing = self._cube_spacing_m / 1e5 if self._cube_spacing_m else ll_res

        logger.debug(f'Output cube spacing: {out_spacing} degrees')
        self._output_spacing = out_spacing

    def add_buffer(self, ll_res, digits=2) -> None:
        """
        Add a fixed buffer to the AOI, accounting for the cube spacing.

        Ensures cube is slighly larger than requested area.
        The AOI will always be in EPSG:4326
        Args:
            ll_res          - weather model lat/lon resolution
            digits          - number of decimal digits to include in the output

        Returns:
            None. Updates self._bounding_box

        Example:
        >>> from RAiDER.models.hrrr import HRRR
        >>> from RAiDER.llreader import BoundingBox
        >>> wm = HRRR()
        >>> aoi = BoundingBox([37, 38, -92, -91])
        >>> aoi.add_buffer(buffer = 1.5 * wm.getLLRes())
        >>> aoi.bounds()
         [36.93, 38.07, -92.07, -90.93]
        """
        from RAiDER.utilFcns import clip_bbox

        ## add an extra buffer around the user specified region
        S, N, W, E = self.bounds()
        buffer = 1.5 * ll_res
        S, N = np.max([S - buffer, -90]), np.min([N + buffer, 90])
        W, E = W - buffer, E + buffer  # TODO: handle dateline crossings

        ## clip the buffered region to a multiple of the spacing
        self.set_output_spacing(ll_res)
        S, N, W, E = clip_bbox([S, N, W, E], self._output_spacing)

        if np.max([np.abs(W), np.abs(E)]) > 180:
            logger.warning('Bounds extend past +/- 180. Results may be incorrect.')

        self._bounding_box = [np.round(a, digits) for a in (S, N, W, E)]


    def calc_buffer_ray(self, direction, lookDir='right', incAngle=30, maxZ=80, digits=2):
        """
        Calculate the buffer for ray tracing. This only needs to be done in the east-west
        direction due to satellite orbits, and only needs extended on the side closest to
        the sensor.

        Args:
            lookDir (str)      - Sensor look direction, can be "right" or "left"
            losAngle (float)   - Incidence angle in degrees
            maxZ (float)       - maximum integration elevation in km
        """
        direction = direction.lower()
        # for isce object
        try:
            lookDir = lookDir.name.lower()
        except AttributeError:
            lookDir = lookDir.lower()

        assert direction in 'asc desc'.split(), f'Incorrect orbital direction: {direction}. Choose asc or desc.'
        assert lookDir in 'right left'.split(), f'Incorrect look direction: {lookDir}. Choose right or left.'

        S, N, W, E = self.bounds()

        # use a small look angle to calculate near range
        lat_max = np.max([np.abs(S), np.abs(N)])
        near = maxZ * np.tan(np.deg2rad(incAngle))
        buffer = near / (np.cos(np.deg2rad(lat_max)) * 100)

        # buffer on the side nearest the sensor
        if (lookDir == 'right' and direction == 'asc') or (lookDir == 'left' and direction == 'desc'):
            W = W - buffer
        else:
            E = E + buffer

        bounds = [np.round(a, digits) for a in (S, N, W, E)]
        if np.max([np.abs(W), np.abs(E)]) > 180:
            logger.warning('Bounds extend past +/- 180. Results may be incorrect.')
        return bounds

    def set_output_directory(self, output_directory) -> None:
        self._output_directory = output_directory

    def set_output_xygrid(self, dst_crs: Union[int, str]=4326) -> None:
        """Define the locations where the delays will be returned."""
        from RAiDER.utilFcns import transform_bbox

        try:
            out_proj = CRS.from_epsg(dst_crs.replace('EPSG:', ''))
        except AttributeError:
            try:
                out_proj = CRS.from_epsg(dst_crs)
            except pyproj.exceptions.CRSError:
                out_proj = dst_crs

        out_snwe = transform_bbox(self.bounds(), src_crs=4326, dest_crs=out_proj)
        logger.debug(f'Output SNWE: {out_snwe}')

        # Build the output grid
        out_spacing = self.get_output_spacing(out_proj)
        self.xpts = np.arange(out_snwe[2], out_snwe[3] + out_spacing, out_spacing)
        self.ypts = np.arange(out_snwe[1], out_snwe[0] - out_spacing, -out_spacing)


class StationFile(AOI):
    """Use a .csv file containing at least Lat, Lon, and optionally Hgt_m columns."""

    def __init__(self, station_file, demFile=None, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:
        super().__init__(cube_spacing_in_m, output_directory)
        self._filename = station_file
        self._demfile = demFile
        self._bounding_box = bounds_from_csv(station_file)
        self._type = 'station_file'

    def readLL(self) -> tuple[np.ndarray, np.ndarray]:
        """Read the station lat/lons from the csv file."""
        df = pd.read_csv(self._filename).drop_duplicates(subset=['Lat', 'Lon'])
        return df['Lat'].to_numpy(), df['Lon'].to_numpy()

    def readZ(self):
        """Read the station heights from the file, or download a DEM if not present."""
        df = pd.read_csv(self._filename).drop_duplicates(subset=['Lat', 'Lon'])
        if 'Hgt_m' in df.columns:
            return df['Hgt_m'].values
        else:
            # Download the DEM
            from RAiDER.dem import download_dem
            from RAiDER.interpolator import interpolateDEM

            demFile = (
                os.path.join(self._output_directory, 'GLO30_fullres_dem.tif')
                if self._demfile is None
                else self._demfile
            )

            download_dem(
                self._bounding_box,
                writeDEM=True,
                dem_path=Path(demFile),
            )

            ## interpolate the DEM to the query points
            z_out0 = interpolateDEM(demFile, self.readLL())
            if np.isnan(z_out0).all():
                raise Exception('DEM interpolation failed. Check DEM bounds and station coords.')
            z_out = np.diag(z_out0)  # the diagonal is the actual stations coordinates

            # write the elevations to the file
            df['Hgt_m'] = z_out
            df.to_csv(self._filename, index=False)
            self._demfile = None
            return z_out


class RasterRDR(AOI):
    """Use a 2-band raster file containing lat/lon coordinates."""

    def __init__(self, lat_file, lon_file=None, *, hgt_file=None, dem_file=None, convention='isce', cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:
        super().__init__(cube_spacing_in_m, output_directory)
        self._type = 'radar_rasters'
        self._latfile = lat_file
        self._lonfile = lon_file

        if (self._latfile is None) and (self._lonfile is None):
            raise ValueError('You need to specify a 2-band file or two single-band files')

        if not os.path.exists(self._latfile):
            raise ValueError(f'{self._latfile} cannot be found!')

        try:
            bpg = bounds_from_latlon_rasters(lat_file, lon_file)
            self._bounding_box, self._proj, self._geotransform = bpg
        except Exception as e:
            raise ValueError(f'Could not read lat/lon rasters: {e}')

        # keep track of the height file, dem and convention
        self._hgtfile = hgt_file
        self._demfile = dem_file
        self._convention = convention

    def readLL(self) -> tuple[np.ndarray, Optional[np.ndarray]]:
        # allow for 2-band lat/lon raster
        from RAiDER.utilFcns import rio_open
        lats, _ = rio_open(Path(self._latfile))

        if self._lonfile is None:
            return lats, None
        else:
            lons, _ = rio_open(Path(self._lonfile))
            return lats, lons

    def readZ(self) -> np.ndarray:
        """Read the heights from the raster file, or download a DEM if not present."""
        from RAiDER.utilFcns import rio_open
        if self._hgtfile is not None and os.path.exists(self._hgtfile):
            logger.info('Using existing heights at: %s', self._hgtfile)
            hgts, _ = rio_open(self._hgtfile)
            return hgts

        else:
            # Download the DEM
            from RAiDER.dem import download_dem
            from RAiDER.interpolator import interpolateDEM

            demFile = (
                os.path.join(self._output_directory, 'GLO30_fullres_dem.tif')
                if self._demfile is None
                else self._demfile
            )

            download_dem(
                self._bounding_box,
                writeDEM=True,
                dem_path=Path(demFile),
            )
            z_out = interpolateDEM(demFile, self.readLL())

            return z_out


class BoundingBox(AOI):
    """Parse a bounding box AOI."""

    def __init__(self, bbox, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:
        super().__init__(cube_spacing_in_m, output_directory)
        self._bounding_box = bbox
        self._type = 'bounding_box'


class GeocodedFile(AOI):
    """Parse a Geocoded file for coordinates."""

    p: RIO.Profile
    _bounding_box: BB.SNWE
    _is_dem: bool

    def __init__(self, path: Path, is_dem=False, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:
        super().__init__(cube_spacing_in_m, output_directory)

        from RAiDER.utilFcns import rio_extents, rio_profile, rio_stats

        self._filename = path
        self.p = rio_profile(path)
        self._bounding_box = rio_extents(self.p)
        self._is_dem = is_dem
        _, self._proj, self._geotransform = rio_stats(path)
        self._type = 'geocoded_file'
        try:
            self.crs = self.p['crs']
        except KeyError:
            self.crs = None

    def readLL(self) -> tuple[np.ndarray, np.ndarray]:
        # ll_bounds are SNWE
        S, N, W, E = self._bounding_box
        w, h = self.p['width'], self.p['height']
        px = (E - W) / w
        py = (N - S) / h
        x = np.array([W + (t * px) for t in range(w)])
        y = np.array([S + (t * py) for t in range(h)])
        X, Y = np.meshgrid(x, y)
        return Y, X  # lats, lons

    def readZ(self):
        """Download a DEM for the file."""
        from RAiDER.dem import download_dem
        from RAiDER.interpolator import interpolateDEM

        demFile = self._filename if self._is_dem else 'GLO30_fullres_dem.tif'
        bbox = self._bounding_box
        _, _ = download_dem(bbox, writeDEM=True, dem_path=Path(demFile))
        z_out = interpolateDEM(demFile, self.readLL())

        return z_out


class Geocube(AOI):
    """Pull lat/lon/height from a georeferenced data cube."""

    def __init__(self, path_cube, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:
        from RAiDER.utilFcns import rio_stats
        super().__init__(cube_spacing_in_m, output_directory)
        self.path = path_cube
        self._type = 'Geocube'
        self._bounding_box = self.get_extent()
        _, self._proj, self._geotransform = rio_stats(path_cube)

    def get_extent(self):
        with xr.open_dataset(self.path) as ds:
            S, N = ds['latitude'].min().item(), ds['latitude'].max().item()
            W, E = ds['longitude'].min().item(), ds['longitude'].max().item()
        return [S, N, W, E]

    ## untested
    def readLL(self) -> tuple[np.ndarray, np.ndarray]:
        with xr.open_dataset(self.path) as ds:
            lats = ds['latitutde'].data()
            lons = ds['longitude'].data()
        Lats, Lons = np.meshgrid(lats, lons)
        return Lats, Lons

    def readZ(self):
        with xr.open_dataset(self.path) as ds:
            heights = ds['heights'].data
        return heights


def bounds_from_latlon_rasters(lat_filestr: str, lon_filestr: str) -> tuple[BB.SNWE, CRS, RIO.GDAL]:
    """
    Parse lat/lon/height inputs and return
    the appropriate outputs.
    """
    from RAiDER.utilFcns import get_file_and_band, rio_stats

    latinfo = get_file_and_band(lat_filestr)
    loninfo = get_file_and_band(lon_filestr)
    lat_stats, lat_proj, lat_gt = rio_stats(latinfo[0], band=latinfo[1])
    lon_stats, lon_proj, lon_gt = rio_stats(loninfo[0], band=loninfo[1])

    assert lat_proj == lon_proj, 'Projection information for Latitude and Longitude files does not match'
    assert lat_gt == lon_gt, 'Affine transform for Latitude and Longitude files does not match'

    # TODO - handle dateline crossing here
    snwe = (lat_stats.min, lat_stats.max,
            lon_stats.min, lon_stats.max)

    if lat_proj is None:
        logger.debug('Assuming lat/lon files are in EPSG:4326')
        lat_proj = CRS.from_epsg(4326)

    return snwe, lat_proj, lat_gt


def bounds_from_csv(station_file):
    """
    station_file should be a comma-delimited file with at least "Lat"
    and "Lon" columns, which should be EPSG: 4326 projection (i.e WGS84).
    """
    stats = pd.read_csv(station_file).drop_duplicates(subset=['Lat', 'Lon'])
    snwe = [stats['Lat'].min(), stats['Lat'].max(), stats['Lon'].min(), stats['Lon'].max()]
    return snwe

1	# noqa: D100
2	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3	#
4	# Author: Jeremy Maurer, Raymond Hogenson & David Bekaert
5	# Copyright 2019, by the California Institute of Technology. ALL RIGHTS
6	# RESERVED. United States Government Sponsorship acknowledged.
7	#
8	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
9	import os	1✔
10	from pathlib import Path	1✔
11	from typing import Optional, Union	1✔
12
13	import numpy as np	1✔
14	import pyproj	1✔
15	import xarray as xr	1✔
16
17
18	try:	1✔
19	import pandas as pd	1✔
20	except ImportError:	×
21	pd = None	×
22
23	from pyproj import CRS	1✔
24
25	from RAiDER.logger import logger	1✔
26	from RAiDER.types import BB, RIO	1✔
27
28
29	class AOI:	1✔
30	"""
31	This instantiates a generic AOI class object.
32
33	Attributes:
34	_bounding_box - S N W E bounding box
35	_proj - pyproj-compatible CRS
36	_type - Type of AOI
37	"""
38
39	def __init__(self, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:	1✔
40	self._output_directory = output_directory	1✔
41	self._bounding_box = None	1✔
42	self._proj = CRS.from_epsg(4326)	1✔
43	self._geotransform = None	1✔
44	self._cube_spacing_m = cube_spacing_in_m	1✔
45
46
47	def __repr__(self):	1✔
48	return f'AOI: {self.__class__.__name__}({self._bounding_box}, {self._type})'	×
49
50	def type(self):	1✔
51	return self._type	1✔
52
53	def bounds(self):	1✔
54	return list(self._bounding_box).copy()	1✔
55
56	def geotransform(self):	1✔
57	return self._geotransform	×
58
59	def projection(self):	1✔
60	return self._proj	1✔
61
62	def get_output_spacing(self, crs=4326):	1✔
63	"""Return the output spacing in desired units."""
64	output_spacing_deg = self._output_spacing	1✔
65	if not isinstance(crs, CRS):	1✔
66	crs = CRS.from_epsg(crs)	1✔
67
68	## convert it to meters users wants a projected coordinate system
69	if all(axis_info.unit_name == 'degree' for axis_info in crs.axis_info):	1✔
70	output_spacing = output_spacing_deg	1✔
71	else:
72	output_spacing = output_spacing_deg * 1e5	1✔
73
74	return output_spacing	1✔
75
76	def set_output_spacing(self, ll_res=None) -> None:	1✔
77	"""Calculate the spacing for the output grid and weather model.
78
79	Use the requested spacing if exists or the weather model grid itself
80
81	Returns:
82	None. Sets self._output_spacing
83	"""
84	assert ll_res or self._cube_spacing_m, 'Must pass lat/lon resolution if _cube_spacing_m is None'	1✔
85
86	out_spacing = self._cube_spacing_m / 1e5 if self._cube_spacing_m else ll_res	1✔
87
88	logger.debug(f'Output cube spacing: {out_spacing} degrees')	1✔
89	self._output_spacing = out_spacing	1✔
90
91	def add_buffer(self, ll_res, digits=2) -> None:	1✔
92	"""
93	Add a fixed buffer to the AOI, accounting for the cube spacing.
94
95	Ensures cube is slighly larger than requested area.
96	The AOI will always be in EPSG:4326
97	Args:
98	ll_res - weather model lat/lon resolution
99	digits - number of decimal digits to include in the output
100
101	Returns:
102	None. Updates self._bounding_box
103
104	Example:
105	>>> from RAiDER.models.hrrr import HRRR
106	>>> from RAiDER.llreader import BoundingBox
107	>>> wm = HRRR()
108	>>> aoi = BoundingBox([37, 38, -92, -91])
109	>>> aoi.add_buffer(buffer = 1.5 * wm.getLLRes())
110	>>> aoi.bounds()
111	[36.93, 38.07, -92.07, -90.93]
112	"""
113	from RAiDER.utilFcns import clip_bbox	1✔
114
115	## add an extra buffer around the user specified region
116	S, N, W, E = self.bounds()	1✔
117	buffer = 1.5 * ll_res	1✔
118	S, N = np.max([S - buffer, -90]), np.min([N + buffer, 90])	1✔
119	W, E = W - buffer, E + buffer # TODO: handle dateline crossings	1✔
120
121	## clip the buffered region to a multiple of the spacing
122	self.set_output_spacing(ll_res)	1✔
123	S, N, W, E = clip_bbox([S, N, W, E], self._output_spacing)	1✔
124
125	if np.max([np.abs(W), np.abs(E)]) > 180:	1✔
126	logger.warning('Bounds extend past +/- 180. Results may be incorrect.')	×
127
128	self._bounding_box = [np.round(a, digits) for a in (S, N, W, E)]	1✔
129
130
131	def calc_buffer_ray(self, direction, lookDir='right', incAngle=30, maxZ=80, digits=2):	1✔
132	"""
133	Calculate the buffer for ray tracing. This only needs to be done in the east-west
134	direction due to satellite orbits, and only needs extended on the side closest to
135	the sensor.
136
137	Args:
138	lookDir (str) - Sensor look direction, can be "right" or "left"
139	losAngle (float) - Incidence angle in degrees
140	maxZ (float) - maximum integration elevation in km
141	"""
142	direction = direction.lower()	1✔
143	# for isce object
144	try:	1✔
145	lookDir = lookDir.name.lower()	1✔
146	except AttributeError:	×
147	lookDir = lookDir.lower()	×
148
149	assert direction in 'asc desc'.split(), f'Incorrect orbital direction: {direction}. Choose asc or desc.'	1✔
150	assert lookDir in 'right left'.split(), f'Incorrect look direction: {lookDir}. Choose right or left.'	1✔
151
152	S, N, W, E = self.bounds()	1✔
153
154	# use a small look angle to calculate near range
155	lat_max = np.max([np.abs(S), np.abs(N)])	1✔
156	near = maxZ * np.tan(np.deg2rad(incAngle))	1✔
157	buffer = near / (np.cos(np.deg2rad(lat_max)) * 100)	1✔
158
159	# buffer on the side nearest the sensor
160	if (lookDir == 'right' and direction == 'asc') or (lookDir == 'left' and direction == 'desc'):	1✔
161	W = W - buffer	1✔
162	else:
163	E = E + buffer	1✔
164
165	bounds = [np.round(a, digits) for a in (S, N, W, E)]	1✔
166	if np.max([np.abs(W), np.abs(E)]) > 180:	1✔
167	logger.warning('Bounds extend past +/- 180. Results may be incorrect.')	×
168	return bounds	1✔
169
170	def set_output_directory(self, output_directory) -> None:	1✔
171	self._output_directory = output_directory	1✔
172
173	def set_output_xygrid(self, dst_crs: Union[int, str]=4326) -> None:	1✔
174	"""Define the locations where the delays will be returned."""
175	from RAiDER.utilFcns import transform_bbox	1✔
176
177	try:	1✔
178	out_proj = CRS.from_epsg(dst_crs.replace('EPSG:', ''))	1✔
179	except AttributeError:	1✔
180	try:	1✔
181	out_proj = CRS.from_epsg(dst_crs)	1✔
182	except pyproj.exceptions.CRSError:	1✔
183	out_proj = dst_crs	1✔
184
185	out_snwe = transform_bbox(self.bounds(), src_crs=4326, dest_crs=out_proj)	1✔
186	logger.debug(f'Output SNWE: {out_snwe}')	1✔
187
188	# Build the output grid
189	out_spacing = self.get_output_spacing(out_proj)	1✔
190	self.xpts = np.arange(out_snwe[2], out_snwe[3] + out_spacing, out_spacing)	1✔
191	self.ypts = np.arange(out_snwe[1], out_snwe[0] - out_spacing, -out_spacing)	1✔
192
193
194	class StationFile(AOI):	1✔
195	"""Use a .csv file containing at least Lat, Lon, and optionally Hgt_m columns."""
196
197	def __init__(self, station_file, demFile=None, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:	1✔
198	super().__init__(cube_spacing_in_m, output_directory)	1✔
199	self._filename = station_file	1✔
200	self._demfile = demFile	1✔
201	self._bounding_box = bounds_from_csv(station_file)	1✔
202	self._type = 'station_file'	1✔
203
204	def readLL(self) -> tuple[np.ndarray, np.ndarray]:	1✔
205	"""Read the station lat/lons from the csv file."""
206	df = pd.read_csv(self._filename).drop_duplicates(subset=['Lat', 'Lon'])	1✔
207	return df['Lat'].to_numpy(), df['Lon'].to_numpy()	1✔
208
209	def readZ(self):	1✔
210	"""Read the station heights from the file, or download a DEM if not present."""
211	df = pd.read_csv(self._filename).drop_duplicates(subset=['Lat', 'Lon'])	1✔
212	if 'Hgt_m' in df.columns:	1✔
213	return df['Hgt_m'].values	1✔
214	else:
215	# Download the DEM
216	from RAiDER.dem import download_dem	×
217	from RAiDER.interpolator import interpolateDEM	×
218
219	demFile = (	×
220	os.path.join(self._output_directory, 'GLO30_fullres_dem.tif')
221	if self._demfile is None
222	else self._demfile
223	)
224
225	download_dem(	×
226	self._bounding_box,
227	writeDEM=True,
228	dem_path=Path(demFile),
229	)
230
231	## interpolate the DEM to the query points
232	z_out0 = interpolateDEM(demFile, self.readLL())	×
233	if np.isnan(z_out0).all():	×
234	raise Exception('DEM interpolation failed. Check DEM bounds and station coords.')	×
235	z_out = np.diag(z_out0) # the diagonal is the actual stations coordinates	×
236
237	# write the elevations to the file
238	df['Hgt_m'] = z_out	×
239	df.to_csv(self._filename, index=False)	×
NEW 240	self._demfile = None	×
241	return z_out	×
242
243
244	class RasterRDR(AOI):	1✔
245	"""Use a 2-band raster file containing lat/lon coordinates."""
246
247	def __init__(self, lat_file, lon_file=None, *, hgt_file=None, dem_file=None, convention='isce', cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:	1✔
248	super().__init__(cube_spacing_in_m, output_directory)	1✔
249	self._type = 'radar_rasters'	1✔
250	self._latfile = lat_file	1✔
251	self._lonfile = lon_file	1✔
252
253	if (self._latfile is None) and (self._lonfile is None):	1✔
254	raise ValueError('You need to specify a 2-band file or two single-band files')	1✔
255
256	if not os.path.exists(self._latfile):	1✔
257	raise ValueError(f'{self._latfile} cannot be found!')	1✔
258
259	try:	1✔
260	bpg = bounds_from_latlon_rasters(lat_file, lon_file)	1✔
261	self._bounding_box, self._proj, self._geotransform = bpg	1✔
262	except Exception as e:	1✔
263	raise ValueError(f'Could not read lat/lon rasters: {e}')	1✔
264
265	# keep track of the height file, dem and convention
266	self._hgtfile = hgt_file	1✔
267	self._demfile = dem_file	1✔
268	self._convention = convention	1✔
269
270	def readLL(self) -> tuple[np.ndarray, Optional[np.ndarray]]:	1✔
271	# allow for 2-band lat/lon raster
272	from RAiDER.utilFcns import rio_open	1✔
273	lats, _ = rio_open(Path(self._latfile))	1✔
274
275	if self._lonfile is None:	1✔
276	return lats, None	×
277	else:
278	lons, _ = rio_open(Path(self._lonfile))	1✔
279	return lats, lons	1✔
280
281	def readZ(self) -> np.ndarray:	1✔
282	"""Read the heights from the raster file, or download a DEM if not present."""
283	from RAiDER.utilFcns import rio_open	×
284	if self._hgtfile is not None and os.path.exists(self._hgtfile):	×
285	logger.info('Using existing heights at: %s', self._hgtfile)	×
286	hgts, _ = rio_open(self._hgtfile)	×
287	return hgts	×
288
289	else:
290	# Download the DEM
291	from RAiDER.dem import download_dem	×
292	from RAiDER.interpolator import interpolateDEM	×
293
294	demFile = (	×
295	os.path.join(self._output_directory, 'GLO30_fullres_dem.tif')
296	if self._demfile is None
297	else self._demfile
298	)
299
300	download_dem(	×
301	self._bounding_box,
302	writeDEM=True,
303	dem_path=Path(demFile),
304	)
305	z_out = interpolateDEM(demFile, self.readLL())	×
306
307	return z_out	×
308
309
310	class BoundingBox(AOI):	1✔
311	"""Parse a bounding box AOI."""
312
313	def __init__(self, bbox, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:	1✔
314	super().__init__(cube_spacing_in_m, output_directory)	1✔
315	self._bounding_box = bbox	1✔
316	self._type = 'bounding_box'	1✔
317
318
319	class GeocodedFile(AOI):	1✔
320	"""Parse a Geocoded file for coordinates."""
321
322	p: RIO.Profile	1✔
323	_bounding_box: BB.SNWE	1✔
324	_is_dem: bool	1✔
325
326	def __init__(self, path: Path, is_dem=False, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:	1✔
327	super().__init__(cube_spacing_in_m, output_directory)	1✔
328
329	from RAiDER.utilFcns import rio_extents, rio_profile, rio_stats	1✔
330
331	self._filename = path	1✔
332	self.p = rio_profile(path)	1✔
333	self._bounding_box = rio_extents(self.p)	1✔
334	self._is_dem = is_dem	1✔
335	_, self._proj, self._geotransform = rio_stats(path)	1✔
336	self._type = 'geocoded_file'	1✔
337	try:	1✔
338	self.crs = self.p['crs']	1✔
339	except KeyError:	×
340	self.crs = None	×
341
342	def readLL(self) -> tuple[np.ndarray, np.ndarray]:	1✔
343	# ll_bounds are SNWE
344	S, N, W, E = self._bounding_box	1✔
345	w, h = self.p['width'], self.p['height']	1✔
346	px = (E - W) / w	1✔
347	py = (N - S) / h	1✔
348	x = np.array([W + (t * px) for t in range(w)])	1✔
349	y = np.array([S + (t * py) for t in range(h)])	1✔
350	X, Y = np.meshgrid(x, y)	1✔
351	return Y, X # lats, lons	1✔
352
353	def readZ(self):	1✔
354	"""Download a DEM for the file."""
355	from RAiDER.dem import download_dem	1✔
356	from RAiDER.interpolator import interpolateDEM	1✔
357
358	demFile = self._filename if self._is_dem else 'GLO30_fullres_dem.tif'	1✔
359	bbox = self._bounding_box	1✔
360	_, _ = download_dem(bbox, writeDEM=True, dem_path=Path(demFile))	1✔
361	z_out = interpolateDEM(demFile, self.readLL())	1✔
362
363	return z_out	1✔
364
365
366	class Geocube(AOI):	1✔
367	"""Pull lat/lon/height from a georeferenced data cube."""
368
369	def __init__(self, path_cube, cube_spacing_in_m: Optional[float]=None, output_directory=os.getcwd()) -> None:	1✔
370	from RAiDER.utilFcns import rio_stats	×
NEW 371	super().__init__(cube_spacing_in_m, output_directory)	×
372	self.path = path_cube	×
373	self._type = 'Geocube'	×
374	self._bounding_box = self.get_extent()	×
375	_, self._proj, self._geotransform = rio_stats(path_cube)	×
376
377	def get_extent(self):	1✔
378	with xr.open_dataset(self.path) as ds:	×
NEW 379	S, N = ds['latitude'].min().item(), ds['latitude'].max().item()	×
NEW 380	W, E = ds['longitude'].min().item(), ds['longitude'].max().item()	×
381	return [S, N, W, E]	×
382
383	## untested
384	def readLL(self) -> tuple[np.ndarray, np.ndarray]:	1✔
385	with xr.open_dataset(self.path) as ds:	×
NEW 386	lats = ds['latitutde'].data()	×
NEW 387	lons = ds['longitude'].data()	×
388	Lats, Lons = np.meshgrid(lats, lons)	×
389	return Lats, Lons	×
390
391	def readZ(self):	1✔
392	with xr.open_dataset(self.path) as ds:	×
NEW 393	heights = ds['heights'].data	×
394	return heights	×
395
396
397	def bounds_from_latlon_rasters(lat_filestr: str, lon_filestr: str) -> tuple[BB.SNWE, CRS, RIO.GDAL]:	1✔
398	"""
399	Parse lat/lon/height inputs and return
400	the appropriate outputs.
401	"""
402	from RAiDER.utilFcns import get_file_and_band, rio_stats	1✔
403
404	latinfo = get_file_and_band(lat_filestr)	1✔
405	loninfo = get_file_and_band(lon_filestr)	1✔
406	lat_stats, lat_proj, lat_gt = rio_stats(latinfo[0], band=latinfo[1])	1✔
407	lon_stats, lon_proj, lon_gt = rio_stats(loninfo[0], band=loninfo[1])	1✔
408
409	assert lat_proj == lon_proj, 'Projection information for Latitude and Longitude files does not match'	1✔
410	assert lat_gt == lon_gt, 'Affine transform for Latitude and Longitude files does not match'	1✔
411
412	# TODO - handle dateline crossing here
413	snwe = (lat_stats.min, lat_stats.max,	1✔
414	lon_stats.min, lon_stats.max)
415
416	if lat_proj is None:	1✔
417	logger.debug('Assuming lat/lon files are in EPSG:4326')	1✔
418	lat_proj = CRS.from_epsg(4326)	1✔
419
420	return snwe, lat_proj, lat_gt	1✔
421
422
423	def bounds_from_csv(station_file):	1✔
424	"""
425	station_file should be a comma-delimited file with at least "Lat"
426	and "Lon" columns, which should be EPSG: 4326 projection (i.e WGS84).
427	"""
428	stats = pd.read_csv(station_file).drop_duplicates(subset=['Lat', 'Lon'])	1✔
429	snwe = [stats['Lat'].min(), stats['Lat'].max(), stats['Lon'].min(), stats['Lon'].max()]	1✔
430	return snwe	1✔

dbekaert / RAiDER / 305eb64a-9a20-4e4e-9c1e-8cc67e861cd6

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous