diff --git a/creating_geometric_models/Face.py b/creating_geometric_models/Face.py
index d6cf998..3fd815a 100644
--- a/creating_geometric_models/Face.py
+++ b/creating_geometric_models/Face.py
@@ -1,22 +1,48 @@
-from collections import defaultdict
+import os
+
 import numpy as np
+import pyvista as pv
+from tqdm import tqdm, trange
 
 
 class Face:
     def __init__(self, vertex, connect):
-
         self.connect = connect
         self.vertex = vertex
         self.local_vid_lookup = {}
         self.ntriangles = self.connect.shape[0]
 
+    @classmethod
+    def from_file(cls, fname):
+        ext = os.path.splitext(fname)[1].lower()
+
+        if ext in [".ts"]:
+            with open(fname, "r") as fid:
+                return cls.from_ts(fid)
+
+        elif ext in [".vtk", ".vtu", ".xdmf", ".ply", ".stl", ".obj"]:
+            return cls.from_pyvista(fname)
+
+        else:
+            raise ValueError(f"Unsupported file format: {ext}")
+
+    @classmethod
+    def from_pyvista(self, fname):
+        mesh = pv.read(fname)
+        mesh = mesh.clean(tolerance=0.1)
+        if mesh.faces is None:
+            mesh = mesh.extract_surface()
+        faces = mesh.faces.reshape((-1, 4))[:, 1:]
+        myFace = Face(mesh.points, faces)
+
+        return myFace
+
     @classmethod
     def from_ts(self, fid):
         surface_name = "undefined"
         xyzl = []
         vid = {}
         trl = []
-        prev_vert = -1
         ivertex = 0
         line = fid.readline()
         if not line:
@@ -72,7 +98,9 @@ def proj(self, sProj):
 
         transformer = Transformer.from_crs("epsg:4326", sProj, always_xy=True)
         print("projecting the node coordinates")
-        self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(self.vertex[:, 0], self.vertex[:, 1])
+        self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(
+            self.vertex[:, 0], self.vertex[:, 1]
+        )
         print(self.vertex)
         print("done projecting")
 
@@ -82,7 +110,9 @@ def convert_projected_to_latlon(self, sProj):
 
         transformer = Transformer.from_crs(sProj, "epsg:4326", always_xy=True)
         print("convert the node coordinates to lat/lon")
-        self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(self.vertex[:, 0], self.vertex[:, 1])
+        self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(
+            self.vertex[:, 0], self.vertex[:, 1]
+        )
         print(self.vertex)
         print("done converting")
 
@@ -109,11 +139,15 @@ def enforce_min_depth(self, depth_min):
         mesh = trimesh.Trimesh(self.vertex, self.connect)
         assert depth_min > 0
         # list vertex of the face boundary
-        unique_edges = mesh.edges[trimesh.grouping.group_rows(mesh.edges_sorted, require_count=1)]
-        if len(unique_edges)==0:
-            raise ValueError("the surface mesh has no boundary edges. "
+        unique_edges = mesh.edges[
+            trimesh.grouping.group_rows(mesh.edges_sorted, require_count=1)
+        ]
+        if len(unique_edges) == 0:
+            raise ValueError(
+                "the surface mesh has no boundary edges. "
                 "This means you extracted an incorrect set of surfaces "
-                "(the surface mesh of a closed region).")
+                "(the surface mesh of a closed region)."
+            )
         boundary_vid = np.array(list(set(list(unique_edges.flatten()))))
         # mask shallow water region
         shallow = np.where(mesh.vertices[:, 2] > -depth_min)
@@ -142,17 +176,38 @@ def __writeTs(self, fname, write_full_vertex_array=True, append=False):
 
         mode = "a" if append else "w"
         with open(fname, mode) as fout:
-            fout.write("GOCAD TSURF 1\nHEADER {\nname:%s\nborder: true\nmesh: false\n*border*bstone: true\n}\nTFACE\n" % (fname))
-            for ivx in vertex_id_2_write:
-                fout.write("VRTX %s %s %s %s\n" % (ivx, self.vertex[ivx - 1, 0], self.vertex[ivx - 1, 1], self.vertex[ivx - 1, 2]))
+            fout.write(
+                "GOCAD TSURF 1\nHEADER {\nname:%s\nborder: true\nmesh: false\n*border*bstone: true\n}\nTFACE\n"
+                % (fname)
+            )
+            for ivx in tqdm(vertex_id_2_write):
+                fout.write(
+                    "VRTX %s %s %s %s\n"
+                    % (
+                        ivx,
+                        self.vertex[ivx - 1, 0],
+                        self.vertex[ivx - 1, 1],
+                        self.vertex[ivx - 1, 2],
+                    )
+                )
 
-            for i in range(self.ntriangles):
-                fout.write("TRGL %d %d %d\n" % (self.connect[i, 0] + 1, self.connect[i, 1] + 1, self.connect[i, 2] + 1))
+            for i in trange(self.ntriangles):
+                fout.write(
+                    "TRGL %d %d %d\n"
+                    % (
+                        self.connect[i, 0] + 1,
+                        self.connect[i, 1] + 1,
+                        self.connect[i, 2] + 1,
+                    )
+                )
             fout.write("END\n")
 
     def __computeNormal(self):
         "compute efficiently the normals"
-        normal = np.cross(self.vertex[self.connect[:, 1], :] - self.vertex[self.connect[:, 0], :], self.vertex[self.connect[:, 2], :] - self.vertex[self.connect[:, 0], :]).astype(float)
+        normal = np.cross(
+            self.vertex[self.connect[:, 1], :] - self.vertex[self.connect[:, 0], :],
+            self.vertex[self.connect[:, 2], :] - self.vertex[self.connect[:, 0], :],
+        ).astype(float)
         norm = np.linalg.norm(normal, axis=1)
         self.normal = np.divide(normal, norm[:, None])
 
@@ -161,11 +216,14 @@ def __writeStl(self, fname, append=False):
         mode = "a" if append else "w"
         with open(fname, mode) as fout:
             fout.write(f"solid {fname}\n")
-            for k in range(self.ntriangles):
+            for k in trange(self.ntriangles):
                 fout.write("facet normal %e %e %e\n" % tuple(self.normal[k, :]))
                 fout.write("outer loop\n")
                 for i in range(0, 3):
-                    fout.write("vertex %.10e %.10e %.10e\n" % tuple(self.vertex[self.connect[k, i], :]))
+                    fout.write(
+                        "vertex %.10e %.10e %.10e\n"
+                        % tuple(self.vertex[self.connect[k, i], :])
+                    )
                 fout.write("endloop\n")
                 fout.write("endfacet\n")
             fout.write(f"endsolid {fname}\n")
@@ -177,12 +235,20 @@ def __writebStl(self, fname):
         fout.seek(80)
         fout.write(struct.pack("<L", self.ntriangles))
         self.__computeNormal()
-        for k in range(self.ntriangles):
+        for k in trange(self.ntriangles):
             fout.write(struct.pack("<3f", *self.normal[k, :]))
             for i in range(0, 3):
                 fout.write(struct.pack("<3f", *self.vertex[self.connect[k, i], :]))
             fout.write(struct.pack("<H", 0))
 
+    def __writeVtk(self, fname):
+        import pyvista as pv
+
+        padding = np.full((self.ntriangles, 1), 3)
+        faces = np.hstack((padding, self.connect))
+        mesh = pv.PolyData(self.vertex, faces)
+        mesh.save(fname)
+
     def write(self, fname, write_full_vertex_array=True, append=False):
         import os
 
@@ -191,6 +257,8 @@ def write(self, fname, write_full_vertex_array=True, append=False):
             self.__writeTs(fname, write_full_vertex_array, append)
         elif ext == ".stl":
             self.__writeStl(fname, append)
+        elif ext == ".vtk":
+            self.__writeVtk(fname)
         elif ext == ".bstl":
             self.__writebStl(fname)
         else:
diff --git a/creating_geometric_models/Grid.py b/creating_geometric_models/Grid.py
index e3d998f..f0b3c2b 100644
--- a/creating_geometric_models/Grid.py
+++ b/creating_geometric_models/Grid.py
@@ -1,7 +1,7 @@
 import numpy as np
 import os
 from netCDF4 import Dataset
-
+from tqdm import tqdm
 
 class Grid:
     def __init__(self, fname, downsample):
@@ -53,9 +53,22 @@ def read_netcdf(self, fname, downsample):
         xvar = self.determine_netcdf_variables(["lon", "x"])
         yvar = self.determine_netcdf_variables(["lat", "y"])
         zvar = self.determine_netcdf_variables(["elevation", "z", "Band1"])
+
+        zvar_obj = fh.variables[zvar]
+        units = getattr(zvar_obj, "units", "").lower()
+
+        scale = 1.0  # default: meters
+
+        if units in ["m", "meter", "meters", "metre", "metres"]:
+            scale = 1.0
+        elif units in ["km", "kilometer", "kilometers"]:
+            scale = 1e3
+        else:
+            print(f"Warning: unknown elevation unit '{units}', assuming meters")
+
         self.x = fh.variables[xvar][0::downsample]
         self.y = fh.variables[yvar][0::downsample]
-        self.z = fh.variables[zvar][0::downsample, 0::downsample].astype(float)
+        self.z = fh.variables[zvar][0::downsample, 0::downsample].astype(float) * scale
 
         # transpose z if z was given as (lon, lat)
         dim_name_y = fh.variables[yvar].dimensions
@@ -64,7 +77,7 @@ def read_netcdf(self, fname, downsample):
             self.z = self.z.T
 
         if np.ma.is_masked(self.z):
-            self.z = np.ma.filled(self.z, float("nan")) * 1e3
+            self.z = np.ma.filled(self.z, float("nan"))
             self.is_sparse = True
         self.compute_nx_ny()
 
@@ -194,7 +207,7 @@ def isolate_hole(self, argHole):
         if argHole:
             x0, x1, y0, y1 = argHole
             print("tagging hole...")
-            for k in range(nconnect):
+            for k in tqdm(range(nconnect)):
                 coords = self.vertex[self.connect[k, :], 0:2]
                 xmin = min(coords[:, 0])
                 xmax = max(coords[:, 0])
@@ -208,12 +221,12 @@ def isolate_hole(self, argHole):
                 raise ValueError("no hole was tagged")
             print("done tagging hole")
 
-    def smooth(self, zrange):
+    def smooth(self, zrange, sig):
         "Smooth bathymetry in range +- zrange"
         "This is useful when preprocessing bathymetry data before intersection"
         import scipy.ndimage as ndimage
 
-        z_smooth = ndimage.gaussian_filter(self.z, sigma=(1, 1), order=0)
+        z_smooth = ndimage.gaussian_filter(self.z, sigma=(sig, sig), order=0)
         ids = np.where(abs(self.z) < zrange)
         self.z[ids] = z_smooth[ids]
 
diff --git a/creating_geometric_models/create_surface_from_rectilinear_grid.py b/creating_geometric_models/create_surface_from_rectilinear_grid.py
index 7fff090..76da909 100755
--- a/creating_geometric_models/create_surface_from_rectilinear_grid.py
+++ b/creating_geometric_models/create_surface_from_rectilinear_grid.py
@@ -9,8 +9,8 @@
 parser = argparse.ArgumentParser(description="create surface from a (possibly sparse, e.g. Slab2.0 dataset) structured dataset (e.g. netcdf)")
 parser.add_argument("input_file", help="netcdf file")
 parser.add_argument("output_file", help="output file (ext in stl, bstl, ts)")
-parser.add_argument("--subsample", nargs=1, type=int, metavar=("onesample_every"), default=[1], help="use only one value every onesample_every in both direction")
-parser.add_argument("--objectname", nargs=1, metavar=("objectname"), help="name of the surface in gocad")
+parser.add_argument("--downsample", type=int, metavar=("onesample_every"), default=1, help="use only one value every onesample_every in both direction")
+parser.add_argument("--objectname", metavar=("objectname"), help="name of the surface in gocad")
 parser.add_argument("--hole", nargs=4, metavar=(("x0"), ("x1"), ("y0"), ("y1")), help="isolate a hole in surface defined by x0<=x<=x1 and y0<=y<=y1 (stl and ts output only)", type=float)
 parser.add_argument("--crop", nargs=4, metavar=(("x0"), ("x1"), ("y0"), ("y1")), help="select only surfaces in x0<=x<=x1 and y0<=y<=y1", type=float)
 parser.add_argument(
@@ -21,23 +21,23 @@
  projname: name of the (projected) Coordinate Reference System (CRS) (e.g. EPSG:32646 for UTM46N)",
 )
 parser.add_argument("--translate", nargs=2, metavar=("x0", "y0"), default=([0, 0]), help="translates all nodes by (x0,y0)", type=float)
-parser.add_argument("--smooth", nargs=1, metavar=("zrange"), help="smooth zgrid for -zrange < z < zrange", type=float)
-parser.add_argument("--change_zero_elevation", nargs=1, metavar=("new_elevation"), help="change elevation of points with z=0, thus avoiding errors when interesecting with sea surface", type=float)
+parser.add_argument("--smooth", nargs=2, metavar=("zrange", "sigma"), help="smooth zgrid for -zrange < z < zrange", type=float)
+parser.add_argument("--change_zero_elevation", metavar=("new_elevation"), help="change elevation of points with z=0, thus avoiding errors when interesecting with sea surface", type=float)
 args = parser.parse_args()
 
 if not args.objectname:
     base = os.path.basename(args.input_file)
     args.objectname = os.path.splitext(base)[0]
 else:
-    args.objectname = args.objectname[0]
+    args.objectname = args.objectname
 
-structured_grid = Grid(args.input_file, args.subsample[0])
+structured_grid = Grid(args.input_file, args.downsample)
 structured_grid.crop(args.crop)
 
 if args.change_zero_elevation:
-    structured_grid.change_zero_elevation(args.change_zero_elevation[0])
+    structured_grid.change_zero_elevation(args.change_zero_elevation)
 if args.smooth:
-    structured_grid.smooth(args.smooth[0])
+    structured_grid.smooth(args.smooth[0], rgs.smooth[1])
 
 structured_grid.generate_vertex()
 
@@ -62,7 +62,7 @@
     for sid in range(nsolid):
         idtr = np.where(structured_grid.solid_id == sid)[0]
         aVid = np.unique(structured_grid.connect[idtr, :].flatten())
-        myFace = Face(vertex=None, connect=structured_grid.connect[idtr, :])
+        myFace = Face(vertex=structured_grid.vertex, connect=structured_grid.connect[idtr, :])
         if structured_grid.is_sparse:
             myFace.reindex(structured_grid.vid_lookup)
-        myFace.write(f"{basename}{sid}{ext}", structured_grid.vertex, write_full_vertex_array=False)
+        myFace.write(f"{basename}{sid}{ext}", write_full_vertex_array=False)
diff --git a/creating_geometric_models/generate_box.py b/creating_geometric_models/generate_box.py
index f01f7c1..fb1bc60 100755
--- a/creating_geometric_models/generate_box.py
+++ b/creating_geometric_models/generate_box.py
@@ -46,6 +46,13 @@
     default=[10e3],
 )
 
+parser.add_argument(
+    "--rotate",
+    metavar=("angle_in_degree"),
+    help="rotate box around it center and Uz",
+    type=float,
+)
+
 parser.add_argument(
     "--shrink",
     nargs=1,
@@ -107,12 +114,26 @@ def read_lat_lon_range_from_netcdf(fname):
     xwidth = xmax - xmin
     ywidth = ymax - ymin
     height = args.zdim[1] - zmin
+
     box_origin = [
         xmin + 0.5 * (1 - args.shrink[0]) * xwidth,
         ymin + 0.5 * (1 - args.shrink[0]) * ywidth,
         args.zdim[0],
     ]
     box_dimension = [args.shrink[0] * xwidth, args.shrink[0] * ywidth, height]
-    geom.add_box(box_origin, box_dimension, args.meshSize[0])
+    box = geom.add_box(box_origin, box_dimension, args.meshSize[0])
+
+    if args.rotate:
+        center = [
+            box_origin[0] + 0.5 * box_dimension[0],
+            box_origin[1] + 0.5 * box_dimension[1],
+            box_origin[2] + 0.5 * box_dimension[2],
+        ]
+
+        angle_rad = np.radians(args.rotate)
+        axis_vector = [0.0, 0.0, 1.0]
+
+        geom.rotate(box, center, angle_rad, axis_vector)
+
     mesh = geom.generate_mesh(dim=2)
     mesh.write(args.output_file)
diff --git a/creating_geometric_models/projTs.py b/creating_geometric_models/project_surface.py
similarity index 69%
rename from creating_geometric_models/projTs.py
rename to creating_geometric_models/project_surface.py
index 678ceed..40007c9 100755
--- a/creating_geometric_models/projTs.py
+++ b/creating_geometric_models/project_surface.py
@@ -3,14 +3,17 @@
 
 parser = argparse.ArgumentParser(description="project ts file")
 parser.add_argument("ts_file", help="ts filename")
-parser.add_argument("--proj", nargs=1, metavar=("projname"), help="transform vertex array to projected system.\
- projname: name of the (projected) Coordinate Reference System (CRS) (e.g. EPSG:32646 for UTM46N)")
+parser.add_argument("--inputCRS", metavar=("projname"), help="input CRS", default = "epsg:4326")
+parser.add_argument("--proj", metavar=("projname"), help="transform vertex array to projected system.\
+ projname: name of the (projected) Coordinate Reference System (CRS) (e.g. EPSG:32646 for UTM46N)", required=True)
+
+
 args = parser.parse_args()
 
 # set projection
 
 from pyproj import Transformer
-transformer = Transformer.from_crs("epsg:4326", args.proj[0], always_xy=True)
+transformer = Transformer.from_crs(args.inputCRS, args.proj, always_xy=True)
 
 # read Ts file
 fid = open(args.ts_file)
diff --git a/creating_geometric_models/refine_and_smooth_mesh.py b/creating_geometric_models/refine_and_smooth_mesh.py
index 635ad31..b338e91 100755
--- a/creating_geometric_models/refine_and_smooth_mesh.py
+++ b/creating_geometric_models/refine_and_smooth_mesh.py
@@ -3,38 +3,95 @@
 import argparse
 import os
 import trimesh
+import numpy as np
 
 parser = argparse.ArgumentParser(description="refine and smooth")
 parser.add_argument("input_file", help="surface in ts file format")
 parser.add_argument(
     "--N",
-    nargs=1,
-    metavar=("N"),
     type=int,
     required=True,
     help="number of initial refine steps",
 )
 parser.add_argument(
     "--P",
-    nargs=1,
-    metavar=("P"),
     type=int,
     required=True,
     help="number of and refine/smoothing steps",
 )
+
+parser.add_argument(
+    "--remesh_first",
+    nargs=1,
+    metavar=("mesh_size"),
+    type=float,
+    help="remesh first the surface with pygalmesh",
+)
+
+parser.add_argument(
+    "--fix_boundary",
+    dest="fix_boundary",
+    action="store_true",
+    help="fix boundary edge when interpolating",
+)
+
+
 args = parser.parse_args()
 
-fid = open(args.input_file)
-myFace = Face.from_ts(fid)
+myFace = Face.from_file(args.input_file)
+
+if args.remesh_first:
+    """
+    # pygalmesh installed with:
+    sudo apt install libcgal-dev libeigen3-dev
+    pip install pygalmesh
+    """
+    import pygalmesh
+    import meshio
 
-a = trimesh.Trimesh(vertices=myFace.vertex, faces=myFace.connect)
+    mesh = meshio.Mesh(points=myFace.vertex, cells=[("triangle", myFace.connect)])
+    mesh.write("trash_me.vtk")
+    mesh_size = args.remesh_first[0]
+    print(f"remeshing with pygalmesh aiming for mesh size {mesh_size}")
+    mesh = pygalmesh.remesh_surface(
+        "trash_me.vtk",
+        max_edge_size_at_feature_edges=mesh_size,
+        min_facet_angle=25,
+        max_radius_surface_delaunay_ball=mesh_size,
+        max_facet_distance=mesh_size,
+        verbose=False,
+    )
+    os.remove("trash_me.vtk")
+    print("done remeshing")
+    a = trimesh.Trimesh(vertices=mesh.points, faces=mesh.cells[0].data)
+else:
+    a = trimesh.Trimesh(vertices=myFace.vertex, faces=myFace.connect)
 
-for i in range(args.N[0]):
+for i in range(args.N):
     a = a.subdivide()
-for i in range(args.P[0]):
+for i in range(args.P):
     a = a.subdivide()
+
+    if args.fix_boundary:
+        # Identify boundary edges using grouping
+        unique_edges = a.edges[
+            trimesh.grouping.group_rows(a.edges_sorted, require_count=1)
+        ]
+        # Extract unique vertices on boundary
+        boundary_vertices = np.unique(unique_edges)
+        a_before_smoothing = a.copy()
+
     a = trimesh.smoothing.filter_taubin(a)
 
+    if args.fix_boundary:
+        # Replace only the internal vertices
+        a.vertices[np.isin(np.arange(len(a.vertices)), boundary_vertices)] = (
+            a_before_smoothing.vertices[
+                np.isin(np.arange(len(a.vertices)), boundary_vertices)
+            ]
+        )
+
+
 myFace = Face(a.vertices, a.faces)
 basename, ext = os.path.splitext(args.input_file)
-myFace.write(f"{basename}_refined_smooth_{args.N[0]}{ext}")
+myFace.write(f"{basename}_refined_smooth_{args.N}.ts")