Implement bilinear interpolation

linear-interpolation/src/CMakeLists.txt

@@ -16,7 +16,13 @@ add_executable(LinearInterpolate main.cpp
         readdata.cpp
         readdata.h
         vecutils.cpp
-        vecutils.h)
+        vecutils.h
+        interpolate.cpp
+        interpolate.h
+        UVGrid.cpp
+        UVGrid.h
+        point.h
+        point.cpp)
 
 execute_process(
     COMMAND nc-config --includedir

linear-interpolation/src/UVGrid.cpp

@@ -0,0 +1,28 @@
+#include <mdspan>
+#include <ranges>
+
+#include "UVGrid.h"
+#include "readdata.h"
+
+using namespace std;
+
+UVGrid::UVGrid() {
+    auto [us, sizeU] = readHydrodynamicU();
+    auto [vs, sizeV] = readHydrodynamicV();
+    // Assuming sizeU == sizeV
+    uvData = views::zip(us, vs) | ranges::to<vector>();
+    uvMatrix = mdspan(uvData.data(), sizeU);
+    tie(times, lats, lons) = readGrid();
+}
+
+double UVGrid::lonStep() const {
+    return lons[1] - lons[0];
+}
+
+double UVGrid::latStep() const {
+    return lats[1]-lats[0];
+}
+
+int UVGrid::timeStep() const {
+    return times[1]-times[0];
+}

linear-interpolation/src/UVGrid.h

@@ -0,0 +1,30 @@
+#ifndef LINEARINTERPOLATE_UVGRID_H
+#define LINEARINTERPOLATE_UVGRID_H
+
+#include <vector>
+#include "vecutils.h"
+#include "point.h"
+
+class UVGrid {
+private:
+    /**
+     * u == Eastward Current Velocity in the Water Column
+     * v == Northward Current Velocity in the Water Column
+     */
+    std::vector<point> uvData;
+public:
+    UVGrid();
+
+    // The step functions assume regular spacing
+    double lonStep() const;
+    double latStep() const;
+    int timeStep() const;
+
+    std::vector<int> times;
+    std::vector<double> lats;
+    std::vector<double> lons;
+
+    arr3d<point> uvMatrix;
+};
+
+#endif //LINEARINTERPOLATE_UVGRID_H

linear-interpolation/src/interpolate.cpp

@@ -0,0 +1,43 @@
+#include <cmath>
+
+#include <ranges>
+#include "interpolate.h"
+
+using namespace std;
+
+// Inspired by https://numerical.recipes/book.html Chapter 3.6
+point bilinearInterpolate(const UVGrid &uvGrid, std::tuple<int, double, double> timeLatLon) {
+    auto [time, lat, lon] = timeLatLon;
+
+    double latStep = uvGrid.latStep();
+    double lonStep = uvGrid.lonStep();
+    int timeStep = uvGrid.timeStep();
+
+    int latIndex = (lat-uvGrid.lats[0])/latStep;
+    int lonIndex = (lon-uvGrid.lons[0])/lonStep;
+    int timeIndex = (time-uvGrid.times[0])/timeStep;
+
+    double timeRatio = (static_cast<double>(time)-uvGrid.times[timeIndex])/timeStep;
+    double latRatio = (lat-uvGrid.lats[latIndex]) / latStep;
+    double lonRatio = (lon-uvGrid.lons[lonIndex]) / lonStep;
+
+    point point = {0, 0};
+    for(int time = 0; time <= 1; time++) {
+        for(int lat = 0; lat <= 1; lat++) {
+            for(int lon = 0; lon <= 1; lon++) {
+                auto vertex = uvGrid.uvMatrix[
+                        timeIndex + time < uvGrid.uvMatrix.extent(0) ? timeIndex + 1 : timeIndex,
+                        latIndex + lat < uvGrid.uvMatrix.extent(1) ? latIndex + 1 : latIndex,
+                        lonIndex + lon < uvGrid.uvMatrix.extent(2) ? lonIndex + 1 : lonIndex
+                ];
+
+                double timeRation = (1 - time)*(1-timeRatio) + time*timeRatio;
+                double latRation = (1 - lat)*(1-latRatio) + lat*latRatio;
+                double lonRation = (1 - lon)*(1-lonRatio) + lon*lonRatio;
+                point += timeRation * latRation * lonRation * vertex;
+            }
+        }
+    }
+
+    return point;
+}

linear-interpolation/src/interpolate.h

@@ -0,0 +1,12 @@
+#ifndef LINEARINTERPOLATE_INTERPOLATE_H
+#define LINEARINTERPOLATE_INTERPOLATE_H
+
+#include <vector>
+
+#include "UVGrid.h"
+
+point bilinearInterpolate(const UVGrid &uvGrid, std::tuple<int, double, double> timeLatLong);
+
+std::vector<double> bilinearInterpolate();
+
+#endif //LINEARINTERPOLATE_INTERPOLATE_H

linear-interpolation/src/main.cpp

@@ -1,18 +1,13 @@
-#include "readdata.h"
-
-#include <vector>
+#include "interpolate.h"
 
 using namespace std;
 
 int main() {
-    auto [vec, size] = readHydrodynamicU();
+    UVGrid uvGrid;
+    auto p = bilinearInterpolate(uvGrid, {392400, 53, -14.5});
 
-    auto arr = std::mdspan(vec.data(), size);
-
-    print3DMatrixSlice(arr, 100);
-
-    auto [times, lats, longs] = readGrid();
-    printContentsOfVec(lats);
+    println("({}, {})", p.first, p.second);
+    p = bilinearInterpolate(uvGrid, {802400, 62, -14.5});
 
     return 0;
 }

linear-interpolation/src/point.cpp

@@ -0,0 +1,11 @@
+#include "point.h"
+
+point operator+(const point& p1, const point& p2) {
+    return {p1.first + p2.first, p1.second + p2.second};
+}
+
+point& operator+=(point& p1, const point& p2) {
+    p1.first += p2.first;
+    p1.second += p2.second;
+    return p1;
+}

linear-interpolation/src/point.h

@@ -0,0 +1,27 @@
+#ifndef LINEARINTERPOLATE_POINT_H
+#define LINEARINTERPOLATE_POINT_H
+
+#include <utility>
+
+using point = std::pair<double, double>; // {u, v}
+
+point operator+(const point& p1, const point& p2);
+
+point& operator+=(point& p1, const point& p2);
+
+template<typename Scalar>
+point operator*(const point& p, Scalar scalar) {
+    return {p.first * scalar, p.second * scalar};
+}
+
+template<typename Scalar>
+point operator*(Scalar scalar, const point& p) {
+    return {p.first * scalar, p.second * scalar};
+}
+
+template<typename Scalar>
+point operator/(const point& p, Scalar scalar) {
+    return {p.first / scalar, p.second / scalar};
+}
+
+#endif //LINEARINTERPOLATE_POINT_H

linear-interpolation/src/readdata.cpp

@@ -55,10 +55,18 @@ pair<vector<double>, std::dextents<std::size_t, 3>> readHydrodynamicU() {
     return get3DMat<double>(vars.find("uo")->second);
 }
 
-tuple<vector<double>, vector<double>, vector<double>> readGrid() {
+pair<vector<double>, std::dextents<std::size_t, 3>> readHydrodynamicV() {
+    netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read);
+
+    multimap< string, NcVar > vars = data.getVars();
+
+    return get3DMat<double>(vars.find("vo")->second);
+}
+
+tuple<vector<int>, vector<double>, vector<double>> readGrid() {
     netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read);
     multimap< string, NcVar > vars = data.getVars();
-    vector<double> time = getVarVector<double>(vars.find("times")->second);
+    vector<int> time = getVarVector<int>(vars.find("times")->second);
     vector<double> longitude = getVarVector<double>(vars.find("longitude")->second);
     vector<double> latitude = getVarVector<double>(vars.find("latitude")->second);
 

linear-interpolation/src/readdata.h

@@ -4,15 +4,21 @@
 #include "vecutils.h"
 
 /**
- * Reads the file hydrodynamic_U.h5
+ * reads the file hydrodynamic_U.h5
  * @return a pair of the data vector of the contents and its dimensions to be used with mdspan
  */
 std::pair<std::vector<double>, std::dextents<std::size_t, 3>> readHydrodynamicU();
 
+/**
+ * reads the file hydrodynamic_V.h5
+ * @return a pair of the data vector of the contents and its dimensions to be used with mdspan
+ */
+std::pair<std::vector<double>, std::dextents<std::size_t, 3>> readHydrodynamicV();
+
 /**
  * Reads the file grid.h5
  * @return a tuple of (times, latitude, longitude)
  */
-std::tuple<std::vector<double>, std::vector<double>, std::vector<double>> readGrid();
+std::tuple<std::vector<int>, std::vector<double>, std::vector<double>> readGrid();
 
 #endif //LINEARINTERPOLATE_READDATA_H

linear-interpolation/src/vecutils.cpp

@@ -12,3 +12,13 @@ void print3DMatrixSlice(const arr3d<double> &arr, int t) {
         println("");
     }
 }
+
+void print3DMatrixSlice(const arr3d<std::pair<double, double>> &arr, int t) {
+    for (int x = 0; x < arr.extent(1); x++) {
+        for (int y = 0; y < arr.extent(2); y++) {
+            auto [u,v] = arr[t,x,y];
+            print("({:>7.4f}, {:>7.4f}) ", u, v);
+        }
+        println("");
+    }
+}

linear-interpolation/src/vecutils.h

@@ -43,5 +43,6 @@ void print3DMatrixSlice(const arr3d<T> &arr, int t) {
 }
 
 void print3DMatrixSlice(const arr3d<double> &arr, int t);
+void print3DMatrixSlice(const arr3d<std::pair<double, double>> &arr, int t);
 
 #endif //LINEARINTERPOLATE_VECUTILS_H