rearrange tropo_delay imports for raider-base

tools/RAiDER/delay.py

@@ -25,11 +25,9 @@
 from RAiDER.logger import logger
 from RAiDER.losreader import build_ray
 from RAiDER.utilFcns import (
-    lla2ecef,  writeResultsToXarray,
-    rio_profile, transformPoints,
+    writeResultsToXarray, transformPoints,
 )
 
-
 ###############################################################################
 def tropo_delay(
         dt,
@@ -62,12 +60,12 @@ def tropo_delay(
     crs = CRS(out_proj)
 
     # Load CRS from weather model file
-    wm_proj = rio_profile(f"netcdf:{weather_model_file}:t")["crs"]
-    if wm_proj is None:
-       logger.warning("WARNING: I can't find a CRS in the weather model file, so I will assume you are using WGS84")
-       wm_proj = CRS.from_epsg(4326)
-    else:
-        wm_proj = CRS.from_wkt(wm_proj.to_wkt())
+    with xarray.load_dataset(weather_model_file) as ds:
+        try:
+            wm_proj = CRS.from_wkt(ds['t'].attrs['crs'].to_wkt())
+        except KeyError:
+            logger.warning("WARNING: I can't find a CRS in the weather model file, so I will assume you are using WGS84")
+            wm_proj = CRS.from_epsg(4326)
 
     # get heights
     if height_levels is None:
@@ -215,6 +213,7 @@ def _build_cube_ray(
     """
     # Get model heights in an array
     # Assumption: All interpolators here are on the same grid
+    T = Transformer.from_crs(4326, 4978, always_xy=True)
     model_zs = interpolators[0].grid[2]
 
     # Create a regular 2D grid
@@ -246,7 +245,7 @@ def _build_cube_ray(
         else:
             llh = [xx, yy, np.full(yy.shape, ht)]
 
-        xyz = np.stack(lla2ecef(llh[1], llh[0], np.full(yy.shape, ht)), axis=-1)
+        xyz = np.stack(T.transform(llh[1], llh[0], np.full(yy.shape, ht)), axis=-1)
 
         # Step 2 - get LOS vectors for targets
         LOS = los.getLookVectors(ht, llh, xyz, yy)
@@ -307,3 +306,114 @@ def _build_cube_ray(
 
     if output_created_here:
         return outputArrs
+
+
+def writeResultsToXarray(dt, xpts, ypts, zpts, crs, wetDelay, hydroDelay, weather_model_file, out_type):
+    '''
+    write a 1-D array to a NETCDF5 file
+    '''
+       # Modify this as needed for NISAR / other projects
+    ds = xarray.Dataset(
+        data_vars=dict(
+            wet=(["z", "y", "x"],
+                 wetDelay,
+                 {"units" : "m",
+                  "description": f"wet {out_type} delay",
+                  # 'crs': crs.to_epsg(),
+                  "grid_mapping": "crs",
+
+                 }),
+            hydro=(["z", "y", "x"],
+                   hydroDelay,
+                   {"units": "m",
+                    # 'crs': crs.to_epsg(),
+                    "description": f"hydrostatic {out_type} delay",
+                    "grid_mapping": "crs",
+                   }),
+        ),
+        coords=dict(
+            x=(["x"], xpts),
+            y=(["y"], ypts),
+            z=(["z"], zpts),
+        ),
+        attrs=dict(
+            Conventions="CF-1.7",
+            title="RAiDER geo cube",
+            source=os.path.basename(weather_model_file),
+            history=str(datetime.utcnow()) + " RAiDER",
+            description=f"RAiDER geo cube - {out_type}",
+            reference_time=str(dt),
+        ),
+    )
+
+    # Write projection system mapping
+    ds["crs"] = int(-2147483647) # dummy placeholder
+    for k, v in crs.to_cf().items():
+        ds.crs.attrs[k] = v
+
+    # Write z-axis information
+    ds.z.attrs["axis"] = "Z"
+    ds.z.attrs["units"] = "m"
+    ds.z.attrs["description"] = "height above ellipsoid"
+
+    # If in degrees
+    if crs.axis_info[0].unit_name == "degree":
+        ds.y.attrs["units"] = "degrees_north"
+        ds.y.attrs["standard_name"] = "latitude"
+        ds.y.attrs["long_name"] = "latitude"
+
+        ds.x.attrs["units"] = "degrees_east"
+        ds.x.attrs["standard_name"] = "longitude"
+        ds.x.attrs["long_name"] = "longitude"
+
+    else:
+        ds.y.attrs["axis"] = "Y"
+        ds.y.attrs["standard_name"] = "projection_y_coordinate"
+        ds.y.attrs["long_name"] = "y-coordinate in projected coordinate system"
+        ds.y.attrs["units"] = "m"
+
+        ds.x.attrs["axis"] = "X"
+        ds.x.attrs["standard_name"] = "projection_x_coordinate"
+        ds.x.attrs["long_name"] = "x-coordinate in projected coordinate system"
+        ds.x.attrs["units"] = "m"
+
+    return ds
+
+
+
+def transformPoints(lats: np.ndarray, lons: np.ndarray, hgts: np.ndarray, old_proj: CRS, new_proj: CRS) -> np.ndarray:
+    '''
+    Transform lat/lon/hgt data to an array of points in a new
+    projection
+
+    Args:
+        lats: ndarray   - WGS-84 latitude (EPSG: 4326)
+        lons: ndarray   - ditto for longitude
+        hgts: ndarray   - Ellipsoidal height in meters
+        old_proj: CRS   - the original projection of the points
+        new_proj: CRS   - the new projection in which to return the points
+
+    Returns:
+        ndarray: the array of query points in the weather model coordinate system (YX)
+    '''
+    # Flags for flipping inputs or outputs
+    if not isinstance(new_proj, CRS):
+        new_proj = CRS.from_epsg(new_proj.lstrip('EPSG:'))
+    if not isinstance(old_proj, CRS):
+        old_proj = CRS.from_epsg(old_proj.lstrip('EPSG:'))
+
+    t = Transformer.from_crs(old_proj, new_proj)
+
+    in_flip = old_proj.axis_info[0].direction
+    out_flip = new_proj.axis_info[0].direction
+
+    if in_flip == 'east':
+        res = t.transform(lons, lats, hgts)
+    else:
+        res = t.transform(lats, lons, hgts)
+
+    if out_flip == 'east':
+        return np.stack((res[1], res[0], res[2]), axis=-1).T
+    else:
+        return np.stack(res, axis=-1).T
+

tools/RAiDER/delayFcns.py

@@ -5,10 +5,15 @@
 # RESERVED. United States Government Sponsorship acknowledged.
 #
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-import multiprocessing as mp
+try:
+    import multiprocessing as mp
+except ImportError:
+    mp = None
+
 import xarray
 
 import numpy as np
+
 from scipy.interpolate import RegularGridInterpolator as Interpolator
 
 from RAiDER.utilFcns import transformPoints
@@ -62,6 +67,9 @@ def make_shared_raw(inarr):
     Make numpy view array of mp.Array
     """
     # Create flat shared array
+    if mp is None:
+        raise ImportError('multiprocessing is not available')
+
     shared_arr = mp.RawArray('d', inarr.size)
     # Create a numpy view of it
     shared_arr_np = np.ndarray(inarr.shape, dtype=np.float64,

tools/RAiDER/losreader.py

@@ -10,15 +10,17 @@
 import datetime
 import shelve
 
-import xml.etree.ElementTree as ET
 import numpy as np
+try:
+    import xml.etree.ElementTree as ET
+except ImportError:
+    ET = None
+
 from abc import ABC
-from scipy.interpolate import interp1d
 
 from RAiDER.utilFcns import (
-    cosd, sind, rio_open, enu2ecef, lla2ecef, ecef2enu, ecef2lla
+    cosd, sind, rio_open, lla2ecef, ecef2lla
 )
-from RAiDER.constants import _ZREF
 
 
 class LOS(ABC):
@@ -478,6 +480,8 @@ def read_ESA_Orbit_file(filename):
     x, y, z: Nt x 1 ndarrays    - x/y/z positions of the sensor at the times t
     vx, vy, vz: Nt x 1 ndarrays - x/y/z velocities of the sensor at the times t
     '''
+    if ET is None:
+        raise ImportError('read_ESA_Orbit_file: cannot import xml.etree.ElementTree')
     tree = ET.parse(filename)
     root = tree.getroot()
     data_block = root[1]