ISD to Kernel support

CHANGELOG.md

@@ -42,6 +42,7 @@ release.
 ### Added
 
 - Added getTargetStatesRanged so a request can be made with start, end, and range of ETs instead of a list [#115](https://github.com/DOI-USGS/SpiceQL/pull/115)
+- Added ck_writer.py and additional support for ISD to kernel generation [#116](https://github.com/DOI-USGS/SpiceQL/pull/116)
 
 ### 1.2.7
 

SpiceQL/include/io.h

@@ -6,6 +6,7 @@
  *
  **/
 
+#include <cstddef>
 #include <optional>
 #include <string>
 #include <vector>
@@ -207,3 +208,28 @@ namespace SpiceQL {
   void writeTextKernel(std::string fileName, std::string type, nlohmann::json &keywords, std::string comment = "");
 
   }
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+/** Returns 0 on success, -1 on error. On error, call writeCkFromBuffersLastError() for the message. */
+int writeCkFromBuffers(
+  const char* path,
+  const double* quats,
+  size_t n_quats,
+  const double* times,
+  size_t n_times,
+  int bodyCode,
+  const char* referenceFrame,
+  const char* segmentId,
+  const char* sclk,
+  const char* lsk,
+  const double* av,
+  size_t n_av,
+  const char* comment
+);
+/** Last error message from writeCkFromBuffers (valid until next writeCk call). */
+const char* writeCkFromBuffersLastError(void);
+#ifdef __cplusplus
+}
+#endif

SpiceQL/src/api.cpp

@@ -33,6 +33,7 @@ namespace SpiceQL {
     vector<string> default_KernelQualities = {"smithed", "reconstructed"};
 
     json aliasMap = {
+      {"A15_METRIC", "apollo"},
       {"AMICA", "amica"},
       {"CHANDRAYAAN-1_M3", "m3"},
       {"CHANDRAYAAN-1_MRFFR", "mrffr"},
@@ -71,6 +72,7 @@ namespace SpiceQL {
       {"MRO_MARCI_UV", "marci"},
       {"MRO_CTX", "ctx"},
       {"MRO_HIRISE", "hirise"},
+      {"MRO_HIRISE_LOOK_DIRECTION", "hirise"},
       {"MRO_CRISM_VNIR", "crism"},
       {"NEAR EARTH ASTEROID RENDEZVOUS", ""},
       {"NH_LORRI", "lorri"},

SpiceQL/src/io.cpp

@@ -1,5 +1,7 @@
 #include <iostream>
 #include <fstream>
+#include <memory>
+#include <sstream>
 
 #include <fmt/format.h>
 #include <ghc/fs_std.hpp>
@@ -10,6 +12,8 @@
 #include "io.h"
 #include "utils.h"
 
+#include <spdlog/spdlog.h>
+
 
 using namespace std;
 
@@ -96,7 +100,6 @@ namespace SpiceQL {
     for(auto &et : times) {
       double sclkdp;
       checkNaifErrors();
-      sce2c_c(bodyCode/1000, et, &sclkdp);
       checkNaifErrors();
       et = sclkdp;
     }
@@ -141,6 +144,22 @@ namespace SpiceQL {
                  vector<vector<double>> stateVelocities,
                  string segmentComment) {
 
+    if (stateTimes.empty() || statePositions.empty()) {
+      throw runtime_error("writeSpk: stateTimes and statePositions must be non-empty.");
+    }
+
+    // NAIF spkw13_c requires segment start time < end time. Single-epoch (e.g. from ISD) has start == end.
+    if (stateTimes.size() == 1) {
+      stateTimes.push_back(stateTimes.front() + 1E-6);
+      statePositions.push_back(statePositions.front());
+      if (!stateVelocities.empty()) {
+        stateVelocities.push_back(stateVelocities.front());
+      }
+    } else if (stateTimes.front() >= stateTimes.back()) {
+      throw runtime_error(
+          "writeSpk: segment start time must be less than end time (got start == end or reversed order).");
+    }
+
     vector<vector<double>> states;
 
     if (stateVelocities.empty()) {
@@ -361,4 +380,102 @@ namespace SpiceQL {
     SPDLOG_TRACE("Text kernel written to {}", fileName);
   }
 
+  namespace {
+    thread_local std::string g_writeCkFromBuffersLastError;
+
+    std::vector<std::string> split(const std::string& s, char delim) {
+      std::vector<std::string> out;
+      std::istringstream ss(s);
+      std::string part;
+      while (std::getline(ss, part, delim)) {
+        auto start = part.find_first_not_of(" \t");
+        if (start == std::string::npos) continue;
+        auto end = part.find_last_not_of(" \t");
+        out.push_back(part.substr(start, end == std::string::npos ? part.size() : end - start + 1));
+      }
+      return out;
+    }
+  }
+
+  extern "C" int writeCkFromBuffers(
+    const char* path,
+    const double* quats,
+    size_t n_quats,
+    const double* times,
+    size_t n_times,
+    int bodyCode,
+    const char* referenceFrame,
+    const char* segmentId,
+    const char* sclk,
+    const char* lsk,
+    const double* av,
+    size_t n_av,
+    const char* comment
+  ) {
+    g_writeCkFromBuffersLastError.clear();
+    if (path == nullptr || quats == nullptr || times == nullptr || n_quats == 0 || n_times == 0) {
+      g_writeCkFromBuffersLastError = "writeCkFromBuffers: path/quats/times must be non-null and non-empty.";
+      return -1;
+    }
+    try {
+      std::string commentStr(comment ? comment : "");
+      if (commentStr.empty()) commentStr = "CK Kernel";
+
+      // sclk: single path or comma-separated list; must keep Kernel objects alive or destructors unload them
+      std::vector<std::unique_ptr<Kernel>> sclkKernels;
+      if (sclk && *sclk) {
+        for (const std::string& sclkPath : split(sclk, ','))
+          sclkKernels.push_back(std::make_unique<Kernel>(sclkPath));
+      }
+      Kernel lskKernel(lsk ? lsk : "");
+
+      int clockId = bodyCode / 1000;
+      std::vector<double> sclkTimes(n_times);
+      for (size_t i = 0; i < n_times; i++) {
+        double sclkdp;
+        checkNaifErrors();
+        sce2c_c(clockId, times[i], &sclkdp);
+        checkNaifErrors();
+        sclkTimes[i] = sclkdp;
+      }
+      checkNaifErrors();
+
+      SpiceInt handle;
+      ckopn_c(path, "CK", (SpiceInt)commentStr.size(), &handle);
+      checkNaifErrors();
+      ckw03_c(
+        handle,
+        sclkTimes[0],
+        sclkTimes[n_times - 1],
+        bodyCode,
+        referenceFrame ? referenceFrame : "",
+        (av != nullptr && n_av > 0) ? SPICETRUE : SPICEFALSE,
+        segmentId ? segmentId : "",
+        (SpiceInt)n_times, sclkTimes.data(),
+        quats,
+        (av != nullptr && n_av > 0) ? av : nullptr,
+        (SpiceInt)n_times,
+        sclkTimes.data()
+      );
+      checkNaifErrors();
+
+      ckcls_c(handle);
+      checkNaifErrors();
+      writeComment(path, commentStr);
+      return 0;
+    } catch (const std::exception& e) {
+      g_writeCkFromBuffersLastError = e.what();
+      reset_c();
+      return -1;
+    } catch (...) {
+      g_writeCkFromBuffersLastError = "writeCkFromBuffers: unknown exception";
+      reset_c();
+      return -1;
+    }
+  }
+
+  extern "C" const char* writeCkFromBuffersLastError(void) {
+    return g_writeCkFromBuffersLastError.c_str();
+  }
+
 }

bindings/python/CMakeLists.txt

@@ -44,6 +44,9 @@ configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/__init__.py
                ${PYSPICEQL_OUTPUT_DIR}/__init__.py
                COPYONLY)
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/ck_writer.py
+               ${PYSPICEQL_OUTPUT_DIR}/ck_writer.py
+               COPYONLY)
 
 # Setup to run setup tools on install 
 install(CODE "execute_process(COMMAND pip install . WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})")

bindings/python/ck_writer.py

@@ -0,0 +1,136 @@
+import ctypes
+import numpy as np
+import os
+import sys
+
+def _find_lib():
+    # Load library
+    try:
+        from . import _pyspiceql as mod
+
+        lib = ctypes.CDLL(mod.__file__)
+        if hasattr(lib, "writeCkFromBuffers"):
+            return lib, mod.__file__, []
+    except (ImportError, AttributeError, OSError) as e:
+        print(f"DEBUG: Internal import failed: {e}")
+        pass
+
+    # Manual search and load library
+    try:
+        current_dir = os.path.dirname(os.path.abspath(__file__))
+        # Look for the .so or .dylib file in the same folder as ck_writer.py
+        for f in os.listdir(current_dir):
+            if f.startswith("_pyspiceql") and (f.endswith(".so") or f.endswith(".dylib")):
+                path = os.path.join(current_dir, f)
+                lib = ctypes.CDLL(path)
+                if hasattr(lib, "writeCkFromBuffers"):
+                    return lib, path, []
+    except Exception as e:
+        print(f"DEBUG: Manual directory search failed: {e}")
+
+    return None, None, []
+
+
+_ck_lib, _ck_lib_path, _ck_lib_search_dirs = _find_lib()
+
+if _ck_lib is not None:
+    try:
+        _ck_lib.writeCkFromBuffers.argtypes = [
+            ctypes.c_char_p,
+            ctypes.POINTER(ctypes.c_double),
+            ctypes.c_size_t,
+            ctypes.POINTER(ctypes.c_double),
+            ctypes.c_size_t,
+            ctypes.c_int,
+            ctypes.c_char_p,
+            ctypes.c_char_p,
+            ctypes.c_char_p,
+            ctypes.c_char_p,
+            ctypes.POINTER(ctypes.c_double),
+            ctypes.c_size_t,
+            ctypes.c_char_p,
+        ]
+        _ck_lib.writeCkFromBuffers.restype = ctypes.c_int
+        if hasattr(_ck_lib, "writeCkFromBuffersLastError"):
+            _ck_lib.writeCkFromBuffersLastError.restype = ctypes.c_char_p
+            _ck_lib.writeCkFromBuffersLastError.argtypes = []
+        _ck_lib_ok = True
+    except AttributeError:
+        _ck_lib_ok = False
+else:
+    _ck_lib_ok = False
+
+
+def write_ck(
+    path,
+    quats,
+    times,
+    body_code,
+    reference_frame,
+    segment_id,
+    sclk,
+    lsk,
+    angular_velocities=None,
+    comment="",
+):
+    if not _ck_lib_ok:
+        dirs = _ck_lib_search_dirs
+        hint = ""
+        if dirs:
+            hint = " Searched: " + ", ".join(dirs) + ". "
+        raise RuntimeError(
+            "writeCkFromBuffers not found. Rebuild SpiceQL and ensure libSpiceQL is on the library path."
+            + hint
+        )
+
+    quats = np.ascontiguousarray(np.asarray(quats, dtype=np.float64))
+    times = np.ascontiguousarray(np.asarray(times, dtype=np.float64))
+    n = quats.shape[0]
+    if quats.ndim != 2 or quats.shape[1] != 4:
+        raise ValueError("quats must have shape (n, 4)")
+    if times.shape[0] != n:
+        raise ValueError("times length must match quats rows")
+
+    if angular_velocities is not None and len(angular_velocities) != 0:
+        angular_velocities = np.ascontiguousarray(np.asarray(angular_velocities, dtype=np.float64))
+        if angular_velocities.shape[0] != n or angular_velocities.shape[1] != 3:
+            raise ValueError("angular_velocities must have shape (n, 3)")
+        av_ptr = angular_velocities.ctypes.data_as(ctypes.POINTER(ctypes.c_double))
+        n_av = n
+    else:
+        av_ptr = None
+        n_av = 0
+
+    def _b(s):
+        """Convert value to bytes"""
+        if s is None:
+            return b""
+        if isinstance(s, str):
+            return s.encode("utf-8")
+        return str(s).encode("utf-8")
+
+    # Accept str or list of str (multiple SCLK kernels) for sclk, comma deliminated
+    sclk_arg = ",".join(sclk) if isinstance(sclk, (list, tuple)) else sclk
+
+    rc = _ck_lib.writeCkFromBuffers(
+        _b(path),
+        quats.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+        n,
+        times.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
+        n,
+        int(body_code),
+        _b(reference_frame),
+        _b(segment_id),
+        _b(sclk_arg),
+        _b(lsk),
+        av_ptr,
+        n_av,
+        _b(comment),
+    )
+    if rc != 0:
+        err = "unknown error"
+        if hasattr(_ck_lib, "writeCkFromBuffersLastError"):
+            msg = _ck_lib.writeCkFromBuffersLastError()
+            if msg:
+                err = msg.decode("utf-8") if isinstance(msg, bytes) else msg
+        raise RuntimeError(f"writeCkFromBuffers failed: {err}")

bindings/python/io.i

@@ -4,4 +4,5 @@
   #include "io.h"
 %}
 
+%ignore writeCkFromBuffers;
 %include "io.h"