DAY 12 #1

Generated difference images and 3D images for different interpolation methods.
Started writing the report.

comparison_2D_data_interpolation.py

@@ -1,10 +1,13 @@
+import os
+
 import numpy as np
 from matplotlib import pyplot as plt
 from matplotlib.cm import ScalarMappable
 from matplotlib.colors import Normalize
 
 from functions_2D import interpolate_discretized_data, plot_contours
 from helpers.save_figure_position import ShowLastPos
+from mayavi import mlab
 
 """
 This script implements a comparative test for 2D interpolation techniques
@@ -13,6 +16,8 @@
 
 block_stride = 1
 height_levels = 50
+dim_x = 10 * 1e3
+dim_y = 10 * 1e3
 
 
 def mean_square_error(original_datagrid, interpolated_datagrid):
@@ -34,31 +39,75 @@ def test_epsilons_for_method(xo, yo, original_datagrid, x, y, datagrid, epsilons
     plt.show()
 
 
-def test_methods(xo, yo, original_datagrid, x, y, datagrid, plot):
+def test_methods(plot):
     def plot_result():
-        fig = plt.figure(figsize=(16, 9))
-        axes = fig.subplots(1, 4, subplot_kw={"aspect": "equal"})
-        plt.tight_layout()
+        # Set up Figure
+        plt.rc("text", usetex=True)
+        plt.rc("font", family="serif")
+        size = 1
+        fig = plt.figure(figsize=(16 * size, 9 * size))
+        fig.tight_layout()
+        axes = fig.subplots(
+            1,
+            4,
+            sharex="all",
+            sharey="all",
+            subplot_kw={
+                "adjustable": "box",
+                "aspect": "equal",
+                "xticks": [i for i in np.linspace(0, dim_x, 6)],
+                "yticks": [i for i in np.linspace(0, dim_y, 6)],
+                "xticklabels": [f"{i:d} km" for i in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
+                "yticklabels": [f"{i:d} km" for i in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
+                "xlim": (0, dim_x),
+                "ylim": (0, dim_y),
+            },
+            gridspec_kw={"hspace": 0.3},
+        )
+
         axes[0].pcolormesh(xo, yo, original_datagrid, cmap=cmap, norm=norm)
+        axes[0].set_title("Original Raster")
         axes[1].pcolormesh(x, y, datagrid, cmap=cmap, norm=norm)
+        axes[1].set_title("Discretized Raster")
         axes[2].pcolormesh(xo, yo, interpolated_datagrid, cmap=cmap, norm=norm)
-        plt.colorbar(colors, ticks=levels, ax=axes[0:3].ravel().tolist())
+        axes[2].set_title(f"Interpolated Raster, using\n{method_name}")
+        plt.colorbar(colors, ticks=levels, ax=axes[0:3].ravel().tolist(), shrink=0.4, aspect=15)
+        mlab.figure(bgcolor=(1, 1, 1))
+        mlab.surf(
+            yo,
+            xo,
+            np.rot90(interpolated_datagrid.T),
+            warp_scale=5,
+            colormap="terrain",
+            vmin=vmin,
+            vmax=vmax,
+        )
+        filename = method_name.replace("\n", " ").replace(" ", "_")
+        mlab.savefig(f"images/differences/{filename}_3D.png", magnification=10)
+        # Use ImageMagick to remove background from image.
+        os.system(
+            f"convert images/differences/{filename}_3D.png -transparent white images/differences/{filename}_3D.png"
+        )
 
         diff = np.nan_to_num(interpolated_datagrid - original_datagrid)
         cmap_diff = plt.get_cmap("coolwarm")
         norm_diff = Normalize(-np.max(abs(diff)), np.max(abs(diff)))
         colors_diff = ScalarMappable(cmap=cmap_diff, norm=norm_diff)
         axes[3].pcolormesh(xo, yo, diff, cmap=cmap_diff, norm=norm_diff)
-        plt.colorbar(colors_diff, ax=axes[3])
-        plt.show()
-
-    for method in [
-        "linear",  # The spline interpolator don't work well
-        "cubic",
-        "rbf_linear",
-        "rbf_thin_plate_spline",
-        "rbf_cubic",
-        "rbf_quintic",
+        axes[3].set_title("Difference of Interpolated\nand Original Raster")
+        plt.colorbar(colors_diff, ax=axes[3], fraction=0.05, pad=0.1)
+        plt.savefig(f"images/differences/{filename})_2D.png")
+
+    for method, method_name in [
+        [
+            "linear",
+            "Delaunay Triangulation\n and Barycentric Interpolation",
+        ],  # The spline interpolators don't work well
+        ["cubic", "Delaunay Triangulation\n and Clough-Tocher scheme"],
+        # ["rbf_linear", "Radial Basis\nLinear Function"],
+        ["rbf_thin_plate_spline", "Radial Basis\nThin Plate Spline"],
+        # ["rbf_cubic", "Radial Basis\nCubic Function"],
+        # ["rbf_quintic", "Radial Basis\nQuintic Function"],
     ]:
         interpolated_datagrid = interpolate_discretized_data(
             x,
@@ -71,13 +120,15 @@ def plot_result():
         print(f"Method: {method}, Mean Error: {mean_square_error(original_datagrid, interpolated_datagrid)}.")
         if plot:
             plot_result()
-    for epsilon, method in [
-        [3, "rbf_multiquadric"],  # Still not sure what value for epsilon is the best
-        [10, "rbf_multiquadric"],
-        [30, "rbf_multiquadric"],
-        [100, "rbf_multiquadric"],
-        [10, "rbf_inverse_multiquadric"],  # Inverse methods don't work well in sparse places.
-        [3.5, "rbf_inverse_quadratic"],
+    for epsilon, method, method_name in [
+        # Still not sure what value for epsilon is the best
+        # [3, "rbf_multiquadric", "Radial Basis\nMultiquadric Function"],
+        # [10, "rbf_multiquadric", "Radial Basis\nMultiquadric Function"],
+        # [30, "rbf_multiquadric", "Radial Basis\nMultiquadric Function"],
+        # [100, "rbf_multiquadric", "Radial Basis\nMultiquadric Function"],
+        # Inverse methods don't work well in sparse places.
+        # [10, "rbf_inverse_multiquadric", "Radial Basis\nInverse Multiquadric Function"],
+        # [3.5, "rbf_inverse_quadratic", "Radial Basis\nInverse Quadratic Function"],
     ]:
         interpolated_datagrid = interpolate_discretized_data(
             x,
@@ -95,23 +146,25 @@ def plot_result():
             plot_result()
 
 
-# for tile in ["NN17"]:
-for tile in ["NN17", "NO33", "NO44", "NO51"]:
+for tile in ["NO44"]:
+    # for tile in ["NN17", "NO33", "NO44", "NO51"]:
     print(f"Tile name: {tile}.")
     original_datagrid = np.loadtxt(f"data/{tile}.asc", skiprows=5)[::-1, :]
     # Each tile is of dimension 10km x 10km, sampled by 50m, thus we have 200 x 200 samples
-    xo = np.linspace(0, 10, original_datagrid.shape[1])
-    yo = np.linspace(0, 10, original_datagrid.shape[0])
-    maximum = int(np.ceil(np.max(original_datagrid)))
-    levels = [i for i in range(0, maximum + height_levels, height_levels)]
-    maximum = levels[-1]
+    xo = np.linspace(0, dim_x, original_datagrid.shape[1])
+    yo = np.linspace(0, dim_y, original_datagrid.shape[0])
+    minimum = int(np.floor(np.min(original_datagrid) / height_levels) * height_levels)
+    maximum = int(np.ceil(np.max(original_datagrid) / height_levels) * height_levels)
+    levels = np.arange(minimum, maximum + 1e-10, height_levels, dtype=int)
 
+    # Set up Colormaps
     cmap = plt.get_cmap("terrain")
-    norm = Normalize(
-        min(np.min(original_datagrid) * 1.1, -10), maximum * 1.1
-    )  # Leave extra 10% for interpolation overshoot
+    vmin = -0.25 * maximum * 1.1
+    vmax = maximum * 1.1
+    norm = Normalize(vmin, vmax)  # Leave extra 10% for interpolation overshoot
     colors = ScalarMappable(norm=norm, cmap=cmap)
+
     datagrid = height_levels * np.floor(original_datagrid[::block_stride, ::block_stride] / height_levels)
-    x = np.linspace(0, 10, datagrid.shape[1])
-    y = np.linspace(0, 10, datagrid.shape[0])
-    test_methods(xo, yo, original_datagrid, x, y, datagrid, plot=False)
+    x = np.linspace(0, dim_x, datagrid.shape[1])
+    y = np.linspace(0, dim_y, datagrid.shape[0])
+    test_methods(plot=True)

contour_example.py

@@ -10,15 +10,16 @@
 
 # Set up Figure
 plt.rc("text", usetex=True)
-plt.rc("font", family="serif")
+plt.rc("font", family="serif", size=25)
 size = 0.7
 fig = plt.figure(figsize=(16 * size, 9 * size))
 fig.tight_layout()
 axes = fig.subplots(1, 3, sharex="all", sharey="all", subplot_kw={"frame_on": False, "xticks": [], "yticks": []})
 
 # Load and plot the original contour
 contour = np.loadtxt("data/RealContour1A.txt")
-axes[0].plot(contour[:, 0], contour[:, 1], color=plt.get_cmap("tab10")(0))
+axes[0].plot(contour[:, 0], contour[:, 1], "o", markersize=1, color=plt.get_cmap("tab10")(0))
+# axes[0].plot(contour[:, 0], contour[:, 1], color=plt.get_cmap("tab10")(0))
 axes[0].set_title("Terraced contour")
 
 cubic_spline = smooth_contour(contour, collinearity_tol=1e-2)
@@ -35,5 +36,5 @@
 # plt.plot(line[:, 0], line[:, 1], "gx", markersize=10)
 # plt.plot(isolated_points[:, 0], isolated_points[:, 1], "ro")
 
-plt.savefig("images/2D_Contour.png")
+# plt.savefig("images/2D_Contour.png")
 plt.show()

function_example.py

@@ -11,7 +11,7 @@
 
 # Set up Figure
 plt.rc("text", usetex=True)
-plt.rc("font", family="serif")
+plt.rc("font", family="serif", size=16)
 size = 0.5
 fig = plt.figure(figsize=(16 * size, 9 * size))
 fig.tight_layout()
@@ -20,7 +20,12 @@
 
 # Load and plot the original function
 line = np.loadtxt("data/Simple_Line1A.txt")
-plt.plot(line[:, 0], line[:, 1])
+plt.plot(line[:100, 0], line[:100, 1], "o", markersize=2)
+plt.legend(["Original line"], frameon=False)
+plt.savefig("images/1D_Function.png")
+plt.show()
+
+fig = plt.figure(figsize=(16 * size, 9 * size))
 
 # Isolate points from collinear segments
 isolated_points = isolate_from_collinear(line)
@@ -31,6 +36,7 @@
 )
 x = np.linspace(line[0, 0], line[-1, 0], 1000)
 y = cubic_spline(x)
+plt.plot(line[:, 0], line[:, 1], "o", markersize=2)
 plt.plot(x, y)
 
 # BEZIER CURVE
@@ -41,7 +47,7 @@
 
 # Plot line points and isolated points
 # plt.plot(line[:, 0], line[:, 1], "gx", markersize=10)
-# plt.plot(isolated_points[:, 0], isolated_points[:, 1], 'ro')
+plt.plot(isolated_points[:, 0], isolated_points[:, 1], "ro", markersize=3)
 
-plt.savefig("images/1D_Function.png")
+plt.savefig("images/1D_Function_Interpolated.png")
 plt.show()

functions_1D.py

@@ -71,7 +71,7 @@ def add_point(points: np.ndarray, point: np.ndarray):
             i += 1
             continue
 
-        if check_half_plane(
+        if points.shape[1] == 2 and check_half_plane(
             path[[firstpoint, lastpoint]], path[[firstpoint - 1, min(lastpoint + 1, path.shape[0] - 1)]]
         ):
             points = add_point(points, 2 / 3 * path[firstpoint] + 1 / 3 * path[lastpoint])

interpolating_2D_contours.py

@@ -13,6 +13,7 @@
 height_levels = 50
 # NN17 - Fort William, NO44 - north of Dundee, NO51 - in St Andrews, NO33 - in Dundee, NH52 - Drumnadrochit
 tile = "NO44"
+tile_name = "South of Forfar"
 dim_x = 10 * 1e3  # Dimensions of loaded data in m
 dim_y = 10 * 1e3  # Dimensions of loaded data in m
 
@@ -31,10 +32,10 @@
     subplot_kw={
         "adjustable": "box",
         "aspect": "equal",
-        "xticks": [x for x in np.linspace(0, dim_x, 6)],
-        "yticks": [y for y in np.linspace(0, dim_y, 6)],
-        "xticklabels": [f"{x:d} km" for x in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
-        "yticklabels": [f"{y:d} km" for y in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
+        "xticks": [i for i in np.linspace(0, dim_x, 6)],
+        "yticks": [i for i in np.linspace(0, dim_y, 6)],
+        "xticklabels": [f"{i:d} km" for i in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
+        "yticklabels": [f"{i:d} km" for i in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
         "xlim": (0, dim_x),
         "ylim": (0, dim_y),
     },
@@ -102,7 +103,7 @@ def plot_contours_wrap(x, y, datagrid, axes, plot_title, levels=None, discretize
     axes[0].set_title(plot_title)
 
 
-plot_contours_wrap(x, y, datagrid, axes[:, 0], "Original Data", levels=levels)
+plot_contours_wrap(x, y, datagrid, axes[:, 0], plot_title=f"{tile_name}\n50m Resolution Raster", levels=levels)
 
 # Lower Spatial Resolution
 low_spatial_res_datagrid = datagrid[::block_stride, ::block_stride]

interpolating_2D_data.py

@@ -37,10 +37,10 @@
     subplot_kw={
         "adjustable": "box",
         "aspect": "equal",
-        "xticks": [x for x in np.linspace(0, dim_x, 6)],
-        "yticks": [y for y in np.linspace(0, dim_y, 6)],
-        "xticklabels": [f"{x:d} km" for x in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
-        "yticklabels": [f"{y:d} km" for y in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
+        "xticks": [i for i in np.linspace(0, dim_x, 6)],
+        "yticks": [i for i in np.linspace(0, dim_y, 6)],
+        "xticklabels": [f"{i:d} km" for i in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
+        "yticklabels": [f"{i:d} km" for i in np.linspace(0, dim_x / 1e3, 6, dtype=int)],
         "xlim": (0, dim_x),
         "ylim": (0, dim_y),
     },

report/Report.aux

@@ -0,0 +1,34 @@
+\relax 
+\@writefile{tdo}{\contentsline {todo}{Add abstract}{1}{}\protected@file@percent }
+\pgfsyspdfmark {pgfid1}{3729359}{34316829}
+\pgfsyspdfmark {pgfid4}{38001205}{34329117}
+\pgfsyspdfmark {pgfid5}{39852597}{34105203}
+\@writefile{toc}{\contentsline {chapter}{\numberline {1}Introduction}{1}{}\protected@file@percent }
+\@writefile{lof}{\addvspace {10\p@ }}
+\@writefile{lot}{\addvspace {10\p@ }}
+\@writefile{tdo}{\contentsline {todo}{Add introduction.}{1}{}\protected@file@percent }
+\pgfsyspdfmark {pgfid6}{3729359}{40394457}
+\pgfsyspdfmark {pgfid9}{38001205}{40406745}
+\pgfsyspdfmark {pgfid10}{39852597}{40182831}
+\@writefile{toc}{\contentsline {chapter}{\numberline {2}Interpolation in 1D}{2}{}\protected@file@percent }
+\@writefile{lof}{\addvspace {10\p@ }}
+\@writefile{lot}{\addvspace {10\p@ }}
+\@writefile{tdo}{\contentsline {todo}{Possibly add a bit of theory to splines.}{2}{}\protected@file@percent }
+\pgfsyspdfmark {pgfid11}{3729359}{38035161}
+\pgfsyspdfmark {pgfid14}{38001205}{38047449}
+\pgfsyspdfmark {pgfid15}{39852597}{37823535}
+\@writefile{lof}{\contentsline {figure}{\numberline {2.1}{\ignorespaces An example of discretized 1D data.\relax }}{2}{}\protected@file@percent }
+\providecommand*\caption@xref[2]{\@setref\relax\@undefined{#1}}
+\newlabel{fig:1D_fun}{{2.1}{2}}
+\@writefile{lof}{\contentsline {figure}{\numberline {2.2}{\ignorespaces An example of curves interpolated from discretized 1D data.\relax }}{3}{}\protected@file@percent }
+\newlabel{fig:1D_fun_interpolated}{{2.2}{3}}
+\@writefile{toc}{\contentsline {chapter}{\numberline {3}Interpolation in 2D}{4}{}\protected@file@percent }
+\@writefile{lof}{\addvspace {10\p@ }}
+\@writefile{lot}{\addvspace {10\p@ }}
+\@writefile{toc}{\contentsline {section}{\numberline {3.1}Parametric Curves}{4}{}\protected@file@percent }
+\@writefile{lof}{\contentsline {figure}{\numberline {3.1}{\ignorespaces An example of parametric curves interpolated from discretized parametric 1D curve in 2D.\relax }}{4}{}\protected@file@percent }
+\newlabel{fig:2D_parametric_curve}{{3.1}{4}}
+\@writefile{toc}{\contentsline {section}{\numberline {3.2}Countour lines}{4}{}\protected@file@percent }
+\newlabel{fig:2D_contour}{{\caption@xref {fig:2D_contour}{ on input line 87}}{5}}
+\@writefile{lof}{\contentsline {figure}{\numberline {3.2}{\ignorespaces Contour interpolation from data with various resolution.\relax }}{5}{}\protected@file@percent }
+\gdef \@abspage@last{7}