diff --git a/pysgems/__init__.py b/pysgems/__init__.py
index 05e28c1..3489430 100644
--- a/pysgems/__init__.py
+++ b/pysgems/__init__.py
@@ -1,11 +1,4 @@
-__name__ = 'pysgems'
-__author__ = 'Robin Thibaut'
+__name__ = "pysgems"
+__author__ = "Robin Thibaut"
 
-from . import algo
-from . import base
-from . import dis
-from . import examples
-from . import io
-from . import plot
-from . import sgems
-from . import utils
+from . import algo, base, dis, examples, io, plot, sgems, utils
diff --git a/pysgems/algo/__init__.py b/pysgems/algo/__init__.py
index 539a0ca..9c073a9 100644
--- a/pysgems/algo/__init__.py
+++ b/pysgems/algo/__init__.py
@@ -1,2 +1 @@
 #  Copyright (c) 2020. Robin Thibaut, Ghent University
-
diff --git a/pysgems/algo/sgalgo.py b/pysgems/algo/sgalgo.py
index ce34179..59120a3 100644
--- a/pysgems/algo/sgalgo.py
+++ b/pysgems/algo/sgalgo.py
@@ -10,9 +10,7 @@
 
 
 class XML(Package):
-    def __init__(self,
-                 project,
-                 algo_dir=None):
+    def __init__(self, project, algo_dir=None):
 
         Package.__init__(self, project)
 
@@ -22,13 +20,14 @@ def __init__(self,
         self.res_dir = self.parent.res_dir
         self.algo_dir = algo_dir
         if self.algo_dir is None:
-            self.algo_dir = jp(self.cwd, 'algorithms')
+            self.algo_dir = jp(self.cwd, "algorithms")
 
-        self.op_file = jp(self.algo_dir, f'{uuid.uuid4().hex}.xml')  # Temporary saves a modified XML
+        # Temporary saves a modified XML
+        self.op_file = jp(self.algo_dir, f"{uuid.uuid4().hex}.xml")
         self.tree = None
         self.root = None
 
-        setattr(self.parent, 'algo', self)
+        setattr(self.parent, "algo", self)
 
     def xml_reader(self, algo_name):
         """
@@ -40,14 +39,16 @@ def xml_reader(self, algo_name):
         except FileNotFoundError:
             pass
 
-        self.tree = ET.parse(jp(self.algo_dir, '{}.xml'.format(algo_name)))
+        self.tree = ET.parse(jp(self.algo_dir, "{}.xml".format(algo_name)))
         self.root = self.tree.getroot()
 
-        name = self.root.find('algorithm').attrib['name']  # Algorithm name
+        name = self.root.find("algorithm").attrib["name"]  # Algorithm name
 
         # By default, replace the grid name by 'computation_grid', and the name by the algorithm name.
-        replace = [['Grid_Name', {'value': 'computation_grid', 'region': ''}],
-                   ['Property_Name', {'value': name}]]
+        replace = [
+            ["Grid_Name", {"value": "computation_grid", "region": ""}],
+            ["Property_Name", {"value": name}],
+        ]
 
         for r in replace:
             try:
@@ -55,8 +56,8 @@ def xml_reader(self, algo_name):
             except AttributeError:
                 pass
 
-        setattr(self.parent.algo, 'tree', self.tree)
-        setattr(self.parent.algo, 'root', self.root)
+        setattr(self.parent.algo, "tree", self.tree)
+        setattr(self.parent.algo, "root", self.root)
 
     def show_tree(self):
         """
@@ -72,19 +73,17 @@ def show_tree(self):
                     elems = list(element)
                     for e in elems:
                         c_list.append(e.tag)
-                        print('//'.join(c_list))
+                        print("//".join(c_list))
                         print(e.attrib)
                         element = list(e)
                         if len(element) == 0:
                             c_list.pop(-1)
         except TypeError:
-            print('No loaded XML file')
+            print("No loaded XML file")
 
-    def xml_update(self, path,
-                   attribute_name=None,
-                   value=None,
-                   new_attribute_dict=None,
-                   show=1):
+    def xml_update(
+        self, path, attribute_name=None, value=None, new_attribute_dict=None, show=1
+    ):
         """
         Given a path in the algorithm XML file, changes the corresponding attribute to the new attribute
         :param path: object path
@@ -95,14 +94,14 @@ def xml_update(self, path,
         """
 
         if new_attribute_dict is not None:
-            if (auto_update is True) and ('property' in new_attribute_dict):
+            if (auto_update is True) and ("property" in new_attribute_dict):
                 # If one property point set needs to be used
-                pp = new_attribute_dict['property']  # Name property
+                pp = new_attribute_dict["property"]  # Name property
                 if pp in self.parent.point_set.columns:
-                    if 'grid' in new_attribute_dict:  # ensure default grid name
-                        new_attribute_dict['grid'] = '{}_grid'.format(pp)
-                    if 'value' in new_attribute_dict:  # ensure default grid name
-                        new_attribute_dict['value'] = '{}_grid'.format(pp)
+                    if "grid" in new_attribute_dict:  # ensure default grid name
+                        new_attribute_dict["grid"] = "{}_grid".format(pp)
+                    if "value" in new_attribute_dict:  # ensure default grid name
+                        new_attribute_dict["value"] = "{}_grid".format(pp)
             self.root.find(path).attrib = new_attribute_dict
             self.tree.write(self.op_file)
 
@@ -111,12 +110,12 @@ def xml_update(self, path,
             self.tree.write(self.op_file)
 
         if show:
-            print('Updated')
+            print("Updated")
             print(self.root.find(path).tag)
             print(self.root.find(path).attrib)
 
-        setattr(self.parent.algo, 'tree', self.tree)
-        setattr(self.parent.algo, 'root', self.root)
+        setattr(self.parent.algo, "tree", self.tree)
+        setattr(self.parent.algo, "root", self.root)
 
     def auto_fill(self):
         """
@@ -138,41 +137,63 @@ def auto_fill(self):
                 trk = list(element.attrib.keys())
 
                 for i in range(len(trv)):
-                    if (trv[i] in self.parent.point_set.columns) \
-                            and ('Variable' or 'Hard_Data' in element.tag):
+                    if (trv[i] in self.parent.point_set.columns) and (
+                        "Variable" or "Hard_Data" in element.tag
+                    ):
                         if trv[i] not in self.parent.object_file_names:
                             self.parent.object_file_names.append(trv[i])
                             try:
-                                if trk[i - 1] == 'grid':  # ensure default grid name
+                                if trk[i - 1] == "grid":  # ensure default grid name
                                     print(element.attrib)
-                                    element.attrib['grid'] = '{}_grid'.format(trv[i])
-                                    self.xml_update('//'.join(c_list), 'grid', '{}_grid'.format(trv[i]))
-                                    print('>>>')
+                                    element.attrib["grid"] = "{}_grid".format(trv[i])
+                                    self.xml_update(
+                                        "//".join(c_list),
+                                        "grid",
+                                        "{}_grid".format(trv[i]),
+                                    )
+                                    print(">>>")
                                     print(element.attrib)
-                                if trk[i - 1] == 'value' and trk[i] == 'property':  # ensure default grid name
+                                # ensure default grid name
+                                if trk[i - 1] == "value" and trk[i] == "property":
                                     print(element.attrib)
-                                    element.attrib['value'] = '{}_grid'.format(trv[i])
-                                    self.xml_update('//'.join(c_list), 'value', '{}_grid'.format(trv[i]))
-                                    print('>>>')
+                                    element.attrib["value"] = "{}_grid".format(trv[i])
+                                    self.xml_update(
+                                        "//".join(c_list),
+                                        "value",
+                                        "{}_grid".format(trv[i]),
+                                    )
+                                    print(">>>")
                                     print(element.attrib)
                             except IndexError:
                                 pass
                             try:
-                                if 'Grid' in elist[-2].tag:
+                                if "Grid" in elist[-2].tag:
                                     tp = list(elist[-2].attrib.keys())
-                                    if 'grid' in tp:
-                                        print('//'.join(c_list[:-2]))
+                                    if "grid" in tp:
+                                        print("//".join(c_list[:-2]))
                                         print(elist[-2].attrib)
-                                        elist[-2].attrib['grid'] = '{}_grid'.format(trv[i])
-                                        self.xml_update(elist[-2].tag, 'grid', '{}_grid'.format(trv[i]))
-                                        print('>>>')
+                                        elist[-2].attrib["grid"] = "{}_grid".format(
+                                            trv[i]
+                                        )
+                                        self.xml_update(
+                                            elist[-2].tag,
+                                            "grid",
+                                            "{}_grid".format(trv[i]),
+                                        )
+                                        print(">>>")
                                         print(elist[-2].attrib)
-                                    if 'value' in tp:
-                                        print('//'.join(c_list[:-2]))
+                                    if "value" in tp:
+                                        print("//".join(c_list[:-2]))
                                         print(elist[-2].attrib)
-                                        elist[-2].attrib['value'] = '{}_grid'.format(trv[i])
-                                        self.xml_update(elist[-2].tag, 'value', '{}_grid'.format(trv[i]))
-                                        print('>>>')
+                                        elist[-2].attrib["value"] = "{}_grid".format(
+                                            trv[i]
+                                        )
+                                        self.xml_update(
+                                            elist[-2].tag,
+                                            "value",
+                                            "{}_grid".format(trv[i]),
+                                        )
+                                        print(">>>")
                                         print(elist[-2].attrib)
                             except IndexError:
                                 pass
@@ -189,26 +210,31 @@ def auto_fill(self):
                             if trv[i] in self.parent.point_set.columns:
                                 if trv[i] not in self.parent.object_file_names:
                                     self.parent.object_file_names.append(trv[i])
-                                    if trk[i] == 'grid':  # ensure default grid name
-                                        print('//'.join(c_list))
+                                    if trk[i] == "grid":  # ensure default grid name
+                                        print("//".join(c_list))
                                         print(e.attrib)
-                                        e.attrib['grid'] = '{}_grid'.format(trv[i])
-                                        self.xml_update('//'.join(c_list), 'grid', '{}_grid'.format(trv[i]))
-                                        print('>>>')
+                                        e.attrib["grid"] = "{}_grid".format(trv[i])
+                                        self.xml_update(
+                                            "//".join(c_list),
+                                            "grid",
+                                            "{}_grid".format(trv[i]),
+                                        )
+                                        print(">>>")
                                         print(e.attrib)
-                                    if trk[i] == 'value':  # ensure default grid name
-                                        print('//'.join(c_list))
+                                    if trk[i] == "value":  # ensure default grid name
+                                        print("//".join(c_list))
                                         print(e.attrib)
-                                        e.attrib['value'] = '{}_grid'.format(trv[i])
-                                        self.xml_update('//'.join(c_list), 'value', '{}_grid'.format(trv[i]))
-                                        print('>>>')
+                                        e.attrib["value"] = "{}_grid".format(trv[i])
+                                        self.xml_update(
+                                            "//".join(c_list),
+                                            "value",
+                                            "{}_grid".format(trv[i]),
+                                        )
+                                        print(">>>")
                                         print(e.attrib)
 
                         element = list(e)
                         if len(element) == 0:
                             c_list.pop(-1)
         except TypeError:
-            print('No loaded XML file')
-
-
-
+            print("No loaded XML file")
diff --git a/pysgems/base/__init__.py b/pysgems/base/__init__.py
index 539a0ca..9c073a9 100644
--- a/pysgems/base/__init__.py
+++ b/pysgems/base/__init__.py
@@ -1,2 +1 @@
 #  Copyright (c) 2020. Robin Thibaut, Ghent University
-
diff --git a/pysgems/base/packbase.py b/pysgems/base/packbase.py
index ee76c1f..a6d5ba5 100644
--- a/pysgems/base/packbase.py
+++ b/pysgems/base/packbase.py
@@ -8,15 +8,15 @@ class PackageInterface(object):
     @abc.abstractmethod
     def parent(self):
         raise NotImplementedError(
-            'must define get_model_dim_arrays in child '
-            'class to use this base class')
+            "must define get_model_dim_arrays in child " "class to use this base class"
+        )
 
     @parent.setter
     @abc.abstractmethod
     def parent(self, name):
         raise NotImplementedError(
-            'must define get_model_dim_arrays in child '
-            'class to use this base class')
+            "must define get_model_dim_arrays in child " "class to use this base class"
+        )
 
 
 class Package(PackageInterface):
diff --git a/pysgems/dis/__init__.py b/pysgems/dis/__init__.py
index 539a0ca..9c073a9 100644
--- a/pysgems/dis/__init__.py
+++ b/pysgems/dis/__init__.py
@@ -1,2 +1 @@
 #  Copyright (c) 2020. Robin Thibaut, Ghent University
-
diff --git a/pysgems/dis/sgdis.py b/pysgems/dis/sgdis.py
index cc08723..4f5f706 100644
--- a/pysgems/dis/sgdis.py
+++ b/pysgems/dis/sgdis.py
@@ -50,28 +50,28 @@ def get_node(r, c, h):
                     [c_sum[j], r_sum[i] - delr[i], l_sum[k] - dell[k]],
                     [c_sum[j] - delc[j], r_sum[i], l_sum[k] - dell[k]],
                     [c_sum[j], r_sum[i], l_sum[k] - dell[k]],
-
                     [c_sum[j] - delc[j], r_sum[i] - delr[i], l_sum[k]],
                     [c_sum[j], r_sum[i] - delr[i], l_sum[k]],
                     [c_sum[j] - delc[j], r_sum[i], l_sum[k]],
-                    [c_sum[j], r_sum[i], l_sum[k]]
+                    [c_sum[j], r_sum[i], l_sum[k]],
                 ]
                 yield get_node(i, j, k), np.array(b), np.mean(b, axis=0)
 
 
 class Discretize(Package):
-
-    def __init__(self,
-                 project,
-                 dx=1,
-                 dy=1,
-                 dz=0,
-                 xo=None,
-                 yo=None,
-                 zo=None,
-                 x_lim=None,
-                 y_lim=None,
-                 z_lim=None):
+    def __init__(
+        self,
+        project,
+        dx=1,
+        dy=1,
+        dz=0,
+        xo=None,
+        yo=None,
+        zo=None,
+        x_lim=None,
+        y_lim=None,
+        z_lim=None,
+    ):
         """
         Constructs the grid geometry. The user can not control directly the number of rows and columns
         but instead inputs the cell size in x and y dimensions.
@@ -93,44 +93,44 @@ def __init__(self,
         # Grid origins
         if xo is None:
             if self.parent.point_set is not None:
-                xs = self.parent.point_set.dataframe['x']
-                xo = np.min(xs) - 4*self.dx
+                xs = self.parent.point_set.dataframe["x"]
+                xo = np.min(xs) - 4 * self.dx
             else:
                 xo = 0
 
         if yo is None:
             if self.parent.point_set is not None:
-                ys = self.parent.point_set.dataframe['y']
-                yo = np.min(ys) - 4*self.dy
+                ys = self.parent.point_set.dataframe["y"]
+                yo = np.min(ys) - 4 * self.dy
             else:
                 yo = 0
 
         if zo is None:
             if self.parent.point_set is not None:
-                zs = self.parent.point_set.dataframe['z']
-                zo = np.min(zs) - 4*self.dz
+                zs = self.parent.point_set.dataframe["z"]
+                zo = np.min(zs) - 4 * self.dz
             else:
                 zo = 0
 
         # Grid limits
         if x_lim is None:
             if self.parent.point_set is not None:
-                xs = self.parent.point_set.dataframe['x']
-                x_lim = np.max(xs) + 4*self.dx
+                xs = self.parent.point_set.dataframe["x"]
+                x_lim = np.max(xs) + 4 * self.dx
             else:
                 x_lim = 1
 
         if y_lim is None:
             if self.parent.point_set is not None:
-                ys = self.parent.point_set.dataframe['y']
-                y_lim = np.max(ys) + 4*self.dy
+                ys = self.parent.point_set.dataframe["y"]
+                y_lim = np.max(ys) + 4 * self.dy
             else:
                 y_lim = 1
 
         if z_lim is None:
             if self.parent.point_set is not None:
-                zs = self.parent.point_set.dataframe['z']
-                z_lim = np.max(zs) + 4*self.dz
+                zs = self.parent.point_set.dataframe["z"]
+                z_lim = np.max(zs) + 4 * self.dz
             else:
                 x_lim = 1
 
@@ -148,9 +148,12 @@ def __init__(self,
         else:
             nlay = 1
 
-        along_r = np.ones(ncol) * self.dx  # Size of each cell along y-dimension - rows
-        along_c = np.ones(nrow) * self.dy  # Size of each cell along x-dimension - columns
-        along_l = np.ones(nlay) * self.dz  # Size of each cell along x-dimension - columns
+        # Size of each cell along y-dimension - rows
+        along_r = np.ones(ncol) * self.dx
+        # Size of each cell along x-dimension - columns
+        along_c = np.ones(nrow) * self.dy
+        # Size of each cell along x-dimension - columns
+        along_l = np.ones(nlay) * self.dz
 
         self.xo = xo
         self.yo = yo
@@ -168,7 +171,7 @@ def __init__(self,
         self.along_c = along_c
         self.along_l = along_l
 
-        setattr(self.parent, 'dis', self)
+        setattr(self.parent, "dis", self)
 
     def my_cell(self, xyz):
         """
@@ -182,18 +185,22 @@ def my_cell(self, xyz):
         # first check if point is within the grid
 
         crit = 0
-        if np.all(rn >= np.array([self.xo, self.yo, self.zo])) and \
-                np.all(rn <= np.array([self.x_lim, self.y_lim, self.z_lim])):
+        if np.all(rn >= np.array([self.xo, self.yo, self.zo])) and np.all(
+            rn <= np.array([self.x_lim, self.y_lim, self.z_lim])
+        ):
             crit = 1
 
         if crit:  # if point inside
             if self.parent.point_set.dimension == 3:  # if 3D
-                dmin = np.min([self.dx, self.dy, self.dz]) / 2  # minimum distance under which a point is in a cell
+                # minimum distance under which a point is in a cell
+                dmin = np.min([self.dx, self.dy, self.dz]) / 2
             else:  # if 2D
-                dmin = np.min([self.dx, self.dy]) / 2  # minimum distance under which a point is in a cell
+                # minimum distance under which a point is in a cell
+                dmin = np.min([self.dx, self.dy]) / 2
 
-            blocks = blocks_from_rc(self.along_c, self.along_r, self.along_l,
-                                    self.xo, self.yo, self.zo)  # cell generator
+            blocks = blocks_from_rc(
+                self.along_c, self.along_r, self.along_l, self.xo, self.yo, self.zo
+            )  # cell generator
 
             # mapping data points to cells:
             # slow but memory-effective method
@@ -203,12 +210,12 @@ def my_cell(self, xyz):
                 c = b[2]
                 dc = np.linalg.norm(rn - c)  # Euclidean distance
                 if dc <= dmin:  # If point is inside cell
-                    print('found 1 cell id in {} s'.format(time.time() - start))
+                    print("found 1 cell id in {} s".format(time.time() - start))
                     return b[0]
                 if dc < vmin:
                     vmin = dc
                     cell = b[0]
-            print('found 1 cell id in {} s'.format(time.time() - start))
+            print("found 1 cell id in {} s".format(time.time() - start))
             return cell
         else:
             return self.parent.nodata
@@ -222,7 +229,8 @@ def compute_cells(self, xyz):
 
         nodes = np.array([self.my_cell(c) for c in xyz])
 
-        np.save(self.cell_file, nodes)  # Save to nodes to avoid recomputing each time
+        # Save to nodes to avoid recomputing each time
+        np.save(self.cell_file, nodes)
 
     def write_hard_data(self, subdata, dis_file=None, cell_file=None, output_dir=None):
         """
@@ -241,26 +249,38 @@ def write_hard_data(self, subdata, dis_file=None, cell_file=None, output_dir=Non
             self.cell_file = cell_file
 
         if self.cell_file is None:
-            self.cell_file = jp(self.parent.res_dir, 'cells.npy')
+            self.cell_file = jp(self.parent.res_dir, "cells.npy")
 
         if output_dir is None:
             output_dir = self.parent.res_dir
 
         if dis_file is None:
-            dis_file = jp(os.path.dirname(self.cell_file), 'grid.dis')
-
-        xyz = subdata[['x', 'y', 'z']].to_numpy()
-
-        npar = np.array([self.dx, self.dy, self.dz,
-                         self.xo, self.yo, self.zo,
-                         self.x_lim, self.y_lim, self.z_lim,
-                         self.nrow, self.ncol, self.nlay])
+            dis_file = jp(os.path.dirname(self.cell_file), "grid.dis")
+
+        xyz = subdata[["x", "y", "z"]].to_numpy()
+
+        npar = np.array(
+            [
+                self.dx,
+                self.dy,
+                self.dz,
+                self.xo,
+                self.yo,
+                self.zo,
+                self.x_lim,
+                self.y_lim,
+                self.z_lim,
+                self.nrow,
+                self.ncol,
+                self.nlay,
+            ]
+        )
 
         if os.path.isfile(dis_file):  # Check previous grid info
             pdis = np.loadtxt(dis_file)
             # If different, recompute data points cell by deleting previous cell file
             if not np.array_equal(pdis, npar):
-                print('New grid found')
+                print("New grid found")
                 try:
                     os.remove(self.cell_file)
                 except FileNotFoundError:
@@ -269,7 +289,7 @@ def write_hard_data(self, subdata, dis_file=None, cell_file=None, output_dir=Non
                     np.savetxt(dis_file, npar)
                     self.compute_cells(xyz)
             else:
-                print('Using previous grid')
+                print("Using previous grid")
                 if not os.path.exists(self.cell_file):
                     self.compute_cells(xyz)
         else:
@@ -282,9 +302,13 @@ def write_hard_data(self, subdata, dis_file=None, cell_file=None, output_dir=Non
         h = subdata.columns.values[-1]
         hd = []
         # fixed nodes = [[node i, value i]....]
-        fixed_nodes = np.array([[data_nodes[dn], subdata[h][dn]] for dn in range(len(data_nodes))])
+        fixed_nodes = np.array(
+            [[data_nodes[dn], subdata[h][dn]] for dn in range(len(data_nodes))]
+        )
         # Deletes points where val == nodata
-        hard_data = np.delete(fixed_nodes, np.where(fixed_nodes == self.parent.nodata)[0], axis=0)
+        hard_data = np.delete(
+            fixed_nodes, np.where(fixed_nodes == self.parent.nodata)[0], axis=0
+        )
         # If data points share the same cell, compute their mean and assign the value to the cell
         for n in unique_nodes:
             where = np.where(hard_data[:, 0] == n)[0]
@@ -294,14 +318,17 @@ def write_hard_data(self, subdata, dis_file=None, cell_file=None, output_dir=Non
 
             fn = hard_data.tolist()
             # For each, feature X, saves a file X.hard
-            cell_values_name = jp(os.path.dirname(self.cell_file), '{}.hard'.format(h))
-            with open(cell_values_name, 'w') as nd:
+            cell_values_name = jp(os.path.dirname(self.cell_file), "{}.hard".format(h))
+            with open(cell_values_name, "w") as nd:
                 nd.write(repr(fn))
             try:
-                shutil.copyfile(cell_values_name,
-                                cell_values_name.replace(os.path.dirname(cell_values_name), output_dir))
+                shutil.copyfile(
+                    cell_values_name,
+                    cell_values_name.replace(
+                        os.path.dirname(cell_values_name), output_dir
+                    ),
+                )
             except shutil.SameFileError:
                 pass
 
-        setattr(self.parent, 'hard_data', hd.append(h))
-
+        setattr(self.parent, "hard_data", hd.append(h))
diff --git a/pysgems/examples/demo_kriging.py b/pysgems/examples/demo_kriging.py
index 7ec3c6f..4480457 100644
--- a/pysgems/examples/demo_kriging.py
+++ b/pysgems/examples/demo_kriging.py
@@ -13,12 +13,12 @@
 def main():
     # %% Initiate sgems pjt
     cwd = os.getcwd()  # Working directory
-    rdir = join_path(cwd, 'results', 'demo_kriging')  # Results directory
-    pjt = sg.Sgems(project_name='sgems_test', project_wd=cwd, res_dir=rdir)
+    rdir = join_path(cwd, "results", "demo_kriging")  # Results directory
+    pjt = sg.Sgems(project_name="sgems_test", project_wd=cwd, res_dir=rdir)
 
     # %% Load data point set
-    data_dir = join_path(cwd, 'datasets', 'demo_kriging')
-    dataset = 'sgems_dataset.eas'
+    data_dir = join_path(cwd, "datasets", "demo_kriging")
+    dataset = "sgems_dataset.eas"
     file_path = join_path(data_dir, dataset)
 
     ps = PointSet(project=pjt, pointset_path=file_path)
@@ -35,9 +35,9 @@ def main():
     print(pjt.point_set.columns)
 
     # %% Load your algorithm xml file in the 'algorithms' folder.
-    algo_dir = join_path(os.path.dirname(cwd), 'algorithms')
+    algo_dir = join_path(os.path.dirname(cwd), "algorithms")
     al = XML(project=pjt, algo_dir=algo_dir)
-    al.xml_reader('kriging')
+    al.xml_reader("kriging")
 
     # %% Show xml structure tree
     al.show_tree()
@@ -45,12 +45,12 @@ def main():
     # %% Modify xml below:
     # By default, the feature grid name of feature X is called 'X_grid'.
     # sgems.xml_update(path, attribute, new value)
-    al.xml_update('Hard_Data', 'grid', 'ag_grid')
-    al.xml_update('Hard_Data', 'property', 'ag')
+    al.xml_update("Hard_Data", "grid", "ag_grid")
+    al.xml_update("Hard_Data", "property", "ag")
 
     # %% Write binary datasets of needed features
     # sgems.export_01(['f1', 'f2'...'fn'])
-    ps.export_01('ag')
+    ps.export_01("ag")
 
     # %% Write python script
     pjt.write_command()
@@ -61,5 +61,5 @@ def main():
     pl.plot_2d(save=True)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/pysgems/examples/demo_sgsim.py b/pysgems/examples/demo_sgsim.py
index eba2767..6cb3a1d 100644
--- a/pysgems/examples/demo_sgsim.py
+++ b/pysgems/examples/demo_sgsim.py
@@ -13,12 +13,12 @@
 def main():
     # %% Initiate sgems pjt
     cwd = os.getcwd()  # Working directory
-    rdir = join_path(cwd, 'results', 'demo_sgsim')  # Results directory
-    pjt = sg.Sgems(project_name='sgsim_test', project_wd=cwd, res_dir=rdir)
+    rdir = join_path(cwd, "results", "demo_sgsim")  # Results directory
+    pjt = sg.Sgems(project_name="sgsim_test", project_wd=cwd, res_dir=rdir)
 
     # %% Load hard data point set
-    data_dir = join_path(cwd, 'datasets', 'demo_sgsim')
-    dataset = 'sgsim_hard_data.eas'
+    data_dir = join_path(cwd, "datasets", "demo_sgsim")
+    dataset = "sgsim_hard_data.eas"
     file_path = join_path(data_dir, dataset)
 
     hd = PointSet(project=pjt, pointset_path=file_path)
@@ -36,9 +36,9 @@ def main():
     print(pjt.point_set.columns)
 
     # %% Load your algorithm xml file in the 'algorithms' folder.
-    algo_dir = join_path(os.path.dirname(cwd), 'algorithms')
+    algo_dir = join_path(os.path.dirname(cwd), "algorithms")
     al = XML(project=pjt, algo_dir=algo_dir)
-    al.xml_reader('sgsim')
+    al.xml_reader("sgsim")
 
     # %% Show xml structure tree
     al.show_tree()
@@ -46,8 +46,8 @@ def main():
     # %% Modify xml below:
     # By default, the feature grid name of feature X is called 'X_grid'.
     # sgems.xml_update(path, attribute, new value)
-    al.xml_update('Assign_Hard_Data', 'value', '1')
-    al.xml_update('Hard_Data', new_attribute_dict={'grid': 'hd_grid', 'property': 'hd'})
+    al.xml_update("Assign_Hard_Data", "value", "1")
+    al.xml_update("Hard_Data", new_attribute_dict={"grid": "hd_grid", "property": "hd"})
 
     # %% Write python script
     pjt.write_command()
@@ -58,5 +58,5 @@ def main():
     pl.plot_2d(save=True)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/pysgems/io/__init__.py b/pysgems/io/__init__.py
index 539a0ca..9c073a9 100644
--- a/pysgems/io/__init__.py
+++ b/pysgems/io/__init__.py
@@ -1,2 +1 @@
 #  Copyright (c) 2020. Robin Thibaut, Ghent University
-
diff --git a/pysgems/io/sgio.py b/pysgems/io/sgio.py
index c250d31..93c4cf2 100644
--- a/pysgems/io/sgio.py
+++ b/pysgems/io/sgio.py
@@ -12,7 +12,7 @@
 def datread(file=None, start=0, end=None):
     # end must be set to None and NOT -1
     """Reads space separated dat file"""
-    with open(file, 'r') as fr:
+    with open(file, "r") as fr:
         lines = np.copy(fr.readlines())[start:end]
         try:
             op = np.array([list(map(float, line.split())) for line in lines])
@@ -61,76 +61,73 @@ def write_point_set(file_name, sub_dataframe, nodata=-999):
     # First, rows with no data occurrence are popped
     sub_dataframe = sub_dataframe[(sub_dataframe != nodata).all(axis=1)]
 
-    xyz = sub_dataframe[['x', 'y', 'z']].to_numpy()
+    xyz = sub_dataframe[["x", "y", "z"]].to_numpy()
 
     pp = sub_dataframe.columns[-1]  # Get name of the property
 
-    grid_name = '{}_grid'.format(pp)
+    grid_name = "{}_grid".format(pp)
 
-    ext = '.sgems'
+    ext = ".sgems"
     if ext not in file_name:
         file_name += ext
 
-    with open(file_name, 'wb') as wb:
-        wb.write(struct.pack('i', int(1.561792946e+9)))  # Magic number
-        wb.write(struct.pack('>i', 10))  # Length of 'Point_set' + 1
+    with open(file_name, "wb") as wb:
+        wb.write(struct.pack("i", int(1.561792946e9)))  # Magic number
+        wb.write(struct.pack(">i", 10))  # Length of 'Point_set' + 1
 
-    with open(file_name, 'a') as wb:
-        wb.write('Point_set')  # Type of file
+    with open(file_name, "a") as wb:
+        wb.write("Point_set")  # Type of file
 
-    with open(file_name, 'ab') as wb:
-        wb.write(struct.pack('>b', 0))  # Signed character 0 after str
+    with open(file_name, "ab") as wb:
+        wb.write(struct.pack(">b", 0))  # Signed character 0 after str
 
-    with open(file_name, 'ab') as wb:
-        wb.write(struct.pack('>i', len(grid_name) + 1))  # Length of 'grid' + 1
+    with open(file_name, "ab") as wb:
+        wb.write(struct.pack(">i", len(grid_name) + 1))  # Length of 'grid' + 1
 
-    with open(file_name, 'a') as wb:
+    with open(file_name, "a") as wb:
         wb.write(grid_name)  # Name of the grid on which points are saved
 
-    with open(file_name, 'ab') as wb:
-        wb.write(struct.pack('>b', 0))  # Signed character 0 after str
+    with open(file_name, "ab") as wb:
+        wb.write(struct.pack(">b", 0))  # Signed character 0 after str
 
-    with open(file_name, 'ab') as wb:
-        wb.write(struct.pack('>i', 100))  # version number
-        wb.write(struct.pack('>i', len(xyz)))  # n data points
-        wb.write(struct.pack('>i', 1))  # n property
+    with open(file_name, "ab") as wb:
+        wb.write(struct.pack(">i", 100))  # version number
+        wb.write(struct.pack(">i", len(xyz)))  # n data points
+        wb.write(struct.pack(">i", 1))  # n property
 
-    with open(file_name, 'ab') as wb:
-        wb.write(struct.pack('>i', len(pp) + 1))  # Length of property name + 1
+    with open(file_name, "ab") as wb:
+        wb.write(struct.pack(">i", len(pp) + 1))  # Length of property name + 1
 
-    with open(file_name, 'a') as wb:
+    with open(file_name, "a") as wb:
         wb.write(pp)  # Property name
 
-    with open(file_name, 'ab') as wb:
-        wb.write(struct.pack('>b', 0))  # Signed character 0 after str
+    with open(file_name, "ab") as wb:
+        wb.write(struct.pack(">b", 0))  # Signed character 0 after str
 
-    with open(file_name, 'ab') as wb:
+    with open(file_name, "ab") as wb:
         for c in xyz:
-            ttt = struct.pack('>fff', c[0], c[1], c[2])  # Coordinates x, y, z
+            ttt = struct.pack(">fff", c[0], c[1], c[2])  # Coordinates x, y, z
             wb.write(ttt)
 
-    with open(file_name, 'ab') as wb:
+    with open(file_name, "ab") as wb:
         for v in sub_dataframe[pp]:
-            wb.write(struct.pack('>f', v))  # Values
+            wb.write(struct.pack(">f", v))  # Values
 
 
-def export_eas(dataframe, filename='dataset'):
+def export_eas(dataframe, filename="dataset"):
     """Exports a Pandas DataFrame to geo-eas format"""
     columns = list(dataframe.columns.values)
-    name = os.path.basename(filename).split('.')[0]
+    name = os.path.basename(filename).split(".")[0]
     header = [name, str(len(columns))] + columns
     ri = dataframe.iterrows()
-    rows = [' '.join(list(map(str, r[1]))) for r in ri]
+    rows = [" ".join(list(map(str, r[1]))) for r in ri]
     lines = header + rows
-    with open(filename+'.eas', 'w') as eas:
-        eas.writelines('\n'.join(lines))
+    with open(filename + ".eas", "w") as eas:
+        eas.writelines("\n".join(lines))
 
 
 class PointSet(Package):
-
-    def __init__(self,
-                 project,
-                 pointset_path=None):
+    def __init__(self, project, pointset_path=None):
 
         Package.__init__(self, project)
 
@@ -144,23 +141,25 @@ def __init__(self,
         self.dataframe = pd.DataFrame(data=self.raw_data, columns=self.columns)
 
         try:
-            self.xyz = self.dataframe[['x', 'y', 'z']].to_numpy()
+            self.xyz = self.dataframe[["x", "y", "z"]].to_numpy()
             self.dimension = 3
         except KeyError:  # Assumes 2D dataset
-            self.dataframe.insert(2, 'z', np.zeros(self.dataframe.shape[0]))
+            self.dataframe.insert(2, "z", np.zeros(self.dataframe.shape[0]))
             self.columns = list(self.dataframe.columns.values)
-            self.xyz = self.dataframe[['x', 'y', 'z']].to_numpy()
+            self.xyz = self.dataframe[["x", "y", "z"]].to_numpy()
             self.dimension = 2
 
-        setattr(self.parent, 'point_set', self)
+        setattr(self.parent, "point_set", self)
 
     # Load sgems dataset
     def loader(self):
         """Parse dataset in GSLIB format"""
         project_info = datread(self.file_path, end=2)  # Name, n features
         project_name = project_info[0][0].lower()  # Project name (lowered)
-        n_features = int(project_info[1][0])  # Number of features len([x, y, f1, f2... fn])
-        head = datread(self.file_path, start=2, end=2 + n_features)  # Name of features (x, y, z, f1...)
+        # Number of features len([x, y, f1, f2... fn])
+        n_features = int(project_info[1][0])
+        # Name of features (x, y, z, f1...)
+        head = datread(self.file_path, start=2, end=2 + n_features)
         columns_name = [h[0].lower() for h in head]  # Column names (lowered)
         data = datread(self.file_path, start=2 + n_features)  # Raw data
         return data, project_name, columns_name
@@ -177,8 +176,11 @@ def export_01(self, features=None):
             features = self.dataframe.columns.values[3:]
 
         for pp in features:
-            subframe = self.dataframe[['x', 'y', 'z', pp]]  # Extract x, y, z, values
+            # Extract x, y, z, values
+            subframe = self.dataframe[["x", "y", "z", pp]]
             ps_name = jp(self.res_dir, pp)  # Path of binary file
             write_point_set(ps_name, subframe)  # Write binary file
-            if pp not in self.parent.object_file_names:  # Adding features name to load them within sgems
+            if (
+                pp not in self.parent.object_file_names
+            ):  # Adding features name to load them within sgems
                 self.parent.object_file_names.append(pp)
diff --git a/pysgems/plot/__init__.py b/pysgems/plot/__init__.py
index 539a0ca..9c073a9 100644
--- a/pysgems/plot/__init__.py
+++ b/pysgems/plot/__init__.py
@@ -1,2 +1 @@
 #  Copyright (c) 2020. Robin Thibaut, Ghent University
-
diff --git a/pysgems/plot/sgplots.py b/pysgems/plot/sgplots.py
index 4961099..4098654 100644
--- a/pysgems/plot/sgplots.py
+++ b/pysgems/plot/sgplots.py
@@ -9,36 +9,62 @@
 
 
 class Plots(Package):
-
     def __init__(self, project):
         Package.__init__(self, project)
         self.name = self.parent.project_name
 
     def plot_coordinates(self):
         try:
-            plt.plot(self.parent.point_set.raw_data[:, 0], self.parent.point_set.raw_data[:, 1], 'ko')
+            plt.plot(
+                self.parent.point_set.raw_data[:, 0],
+                self.parent.point_set.raw_data[:, 1],
+                "ko",
+            )
         except:
             pass
         try:
-            plt.xticks(np.cumsum(self.parent.dis.along_r) + self.parent.dis.xo - self.parent.dis.dx, labels=[])
-            plt.yticks(np.cumsum(self.parent.dis.along_c) + self.parent.dis.yo - self.parent.dis.dy, labels=[])
+            plt.xticks(
+                np.cumsum(self.parent.dis.along_r)
+                + self.parent.dis.xo
+                - self.parent.dis.dx,
+                labels=[],
+            )
+            plt.yticks(
+                np.cumsum(self.parent.dis.along_c)
+                + self.parent.dis.yo
+                - self.parent.dis.dy,
+                labels=[],
+            )
         except:
             pass
 
-        plt.grid('blue')
+        plt.grid("blue")
         plt.show()
 
     def plot_2d(self, res_file=None, save=False):
         """Rudimentary 2D plot"""
         if res_file is None:
-            res_file = jp(self.parent.res_dir, 'results.grid')
+            res_file = jp(self.parent.res_dir, "results.grid")
         matrix = datread(res_file, start=3)
         matrix = np.where(matrix == -9966699, np.nan, matrix)
         matrix = matrix.reshape((self.parent.dis.nrow, self.parent.dis.ncol))
-        extent = (self.parent.dis.xo, self.parent.dis.x_lim, self.parent.dis.yo, self.parent.dis.y_lim)
-        plt.imshow(np.flipud(matrix), cmap='coolwarm', extent=extent)
-        plt.plot(self.parent.point_set.raw_data[:, 0], self.parent.point_set.raw_data[:, 1], 'k+', markersize=1, alpha=.7)
+        extent = (
+            self.parent.dis.xo,
+            self.parent.dis.x_lim,
+            self.parent.dis.yo,
+            self.parent.dis.y_lim,
+        )
+        plt.imshow(np.flipud(matrix), cmap="coolwarm", extent=extent)
+        plt.plot(
+            self.parent.point_set.raw_data[:, 0],
+            self.parent.point_set.raw_data[:, 1],
+            "k+",
+            markersize=1,
+            alpha=0.7,
+        )
         plt.colorbar()
         if save:
-            plt.savefig(jp(self.parent.res_dir, 'results.png'), bbox_inches='tight', dpi=300)
+            plt.savefig(
+                jp(self.parent.res_dir, "results.png"), bbox_inches="tight", dpi=300
+            )
         plt.show()
diff --git a/pysgems/script_templates/script_template.py b/pysgems/script_templates/script_template.py
index 94a16ee..4818efd 100644
--- a/pysgems/script_templates/script_template.py
+++ b/pysgems/script_templates/script_template.py
@@ -5,14 +5,14 @@
 nodata = -9966699
 
 os.chdir("RES_DIR")
-sgems.execute('DeleteObjects computation_grid')
-sgems.execute('DeleteObjects PROJECT_NAME')
-sgems.execute('DeleteObjects finished')
+sgems.execute("DeleteObjects computation_grid")
+sgems.execute("DeleteObjects PROJECT_NAME")
+sgems.execute("DeleteObjects finished")
 
 for file in OBJECT_FILES:
-    sgems.execute('LoadObjectFromFile  {}::All'.format(file))
+    sgems.execute("LoadObjectFromFile  {}::All".format(file))
 
-sgems.execute('NewCartesianGrid  computation_grid::GRID')
+sgems.execute("NewCartesianGrid  computation_grid::GRID")
 
-#~sgems.execute('RunGeostatAlgorithm  ALGORITHM_NAME::/GeostatParamUtils/XML::ALGORITHM_XML')
-#~sgems.execute('SaveGeostatGrid  computation_grid::FEATURE_OUTPUT.grid::gslib::0::OUTPUT_LIST')
\ No newline at end of file
+# ~sgems.execute('RunGeostatAlgorithm  ALGORITHM_NAME::/GeostatParamUtils/XML::ALGORITHM_XML')
+# ~sgems.execute('SaveGeostatGrid  computation_grid::FEATURE_OUTPUT.grid::gslib::0::OUTPUT_LIST')
diff --git a/pysgems/sgems/__init__.py b/pysgems/sgems/__init__.py
index 539a0ca..9c073a9 100644
--- a/pysgems/sgems/__init__.py
+++ b/pysgems/sgems/__init__.py
@@ -1,2 +1 @@
 #  Copyright (c) 2020. Robin Thibaut, Ghent University
-
diff --git a/pysgems/sgems/sg.py b/pysgems/sgems/sg.py
index 536b202..c4b63a5 100644
--- a/pysgems/sgems/sg.py
+++ b/pysgems/sgems/sg.py
@@ -10,30 +10,39 @@
 
 
 class Sgems:
-
-    def __init__(self,
-                 project_name='sgems_test',
-                 project_wd='',
-                 res_dir='',
-                 script_dir='',
-                 exe_name='',
-                 nodata=-9966699,
-                 check_env=True,
-                 **kwargs):
+    def __init__(
+        self,
+        project_name="sgems_test",
+        project_wd="",
+        res_dir="",
+        script_dir="",
+        exe_name="",
+        nodata=-9966699,
+        check_env=True,
+        **kwargs
+    ):
 
         if check_env:
             # First check if sgems installation files are in the user environment variables
-            gstl_home = os.environ.get('GSTLAPPLIHOME')
+            gstl_home = os.environ.get("GSTLAPPLIHOME")
             if not gstl_home:
-                warnings.warn('GSTLAPPLIHOME environment variable does not exist')
+                warnings.warn("GSTLAPPLIHOME environment variable does not exist")
             else:
-                path = os.getenv('Path')
+                path = os.getenv("Path")
                 if gstl_home not in path:
-                    warnings.warn('Variable {} does not exist in Path environment variable'.format(gstl_home))
+                    warnings.warn(
+                        "Variable {} does not exist in Path environment variable".format(
+                            gstl_home
+                        )
+                    )
                 if not exe_name:  # If no sgems exe file name is provided,
                     # checks for sgems exe file in the GSTLAPPLIHOME path
                     for file in os.listdir(gstl_home):
-                        if file.endswith(".exe") and ('sgems' in file) and ('uninstall' not in file):
+                        if (
+                            file.endswith(".exe")
+                            and ("sgems" in file)
+                            and ("uninstall" not in file)
+                        ):
                             exe_name = file
 
         self.project_name = project_name
@@ -46,7 +55,11 @@ def __init__(self,
         # result directory generated according to project and algorithm name
         if self.res_dir is None:
             # Generate result directory if none is given
-            self.res_dir = jp(self.project_wd, 'results', '_'.join([self.project_name, uuid.uuid1().hex]))
+            self.res_dir = jp(
+                self.project_wd,
+                "results",
+                "_".join([self.project_name, uuid.uuid1().hex]),
+            )
         if not os.path.exists(self.res_dir):
             os.makedirs(self.res_dir)
 
@@ -58,88 +71,101 @@ def __init__(self,
         self.nodata = nodata
 
         self.object_file_names = []  # List of features name needed for the algorithm
-        self.command_name = ''
+        self.command_name = ""
 
         if not script_dir:
             dir_path = os.path.abspath(__file__ + "/../../")
-            self.template_file = jp(dir_path, 'script_templates/script_template.py')  # Python template file path
+            # Python template file path
+            self.template_file = jp(dir_path, "script_templates/script_template.py")
 
     def write_command(self):
         """
         Write python script that sgems will run.
         """
 
-        self.command_name = jp(self.res_dir, '{}_commands.py'.format(self.project_name))
+        self.command_name = jp(self.res_dir, "{}_commands.py".format(self.project_name))
 
-        run_algo_flag = ''  # This empty str will replace the # in front of the commands meant to execute sgems
+        # This empty str will replace the # in front of the commands meant to execute sgems
+        run_algo_flag = ""
         # within its python environment
         try:
-            name = self.algo.root.find('algorithm').attrib['name']  # Algorithm name
+            name = self.algo.root.find("algorithm").attrib["name"]  # Algorithm name
             try:
                 # When performing simulations, sgems automatically add '__realn'
                 # to the name of the nth generated property.
-                nr = int(self.algo.root.find('Nb_Realizations').attrib['value'])
-                name_op = '::'.join([name + '__real' + str(i) for i in range(nr)])
+                nr = int(self.algo.root.find("Nb_Realizations").attrib["value"])
+                name_op = "::".join([name + "__real" + str(i) for i in range(nr)])
             except AttributeError:
                 name_op = name
 
             with open(self.algo.op_file) as alx:  # Remove unwanted \n
-                algo_xml = alx.read().strip('\n')
+                algo_xml = alx.read().strip("\n")
 
         except AttributeError or FileNotFoundError:
-            name = 'None'
+            name = "None"
             name_op = name
-            algo_xml = 'None'
-            run_algo_flag = '#'  # If no algorithm loaded, then just loads the data
-
-        sgrid = [self.dis.ncol, self.dis.nrow, self.dis.nlay,
-                 self.dis.dx, self.dis.dy, self.dis.dz,
-                 self.dis.xo, self.dis.yo, self.dis.zo]  # Grid information
-        grid = joinlist('::', sgrid)  # Grid in sgems format
-
-        sgems_files = [sf + '.sgems' for sf in self.object_file_names]
+            algo_xml = "None"
+            run_algo_flag = "#"  # If no algorithm loaded, then just loads the data
+
+        sgrid = [
+            self.dis.ncol,
+            self.dis.nrow,
+            self.dis.nlay,
+            self.dis.dx,
+            self.dis.dy,
+            self.dis.dz,
+            self.dis.xo,
+            self.dis.yo,
+            self.dis.zo,
+        ]  # Grid information
+        grid = joinlist("::", sgrid)  # Grid in sgems format
+
+        sgems_files = [sf + ".sgems" for sf in self.object_file_names]
 
         # The list below is the list of flags that will be replaced in the sgems python script
         # TODO: add option to change output file name (now default 'results.grid')
-        params = [[run_algo_flag, '#~'],
-                  [self.res_dir.replace('\\', '//'), 'RES_DIR'],  # for sgems convention...
-                  [grid, 'GRID'],
-                  [self.project_name, 'PROJECT_NAME'],
-                  ['results', 'FEATURE_OUTPUT'],  # results.grid = output file
-                  [name, 'ALGORITHM_NAME'],
-                  [name_op, 'OUTPUT_LIST'],
-                  [algo_xml, 'ALGORITHM_XML'],
-                  [str(sgems_files), 'OBJECT_FILES']]
+        params = [
+            [run_algo_flag, "#~"],
+            # for sgems convention...
+            [self.res_dir.replace("\\", "//"), "RES_DIR"],
+            [grid, "GRID"],
+            [self.project_name, "PROJECT_NAME"],
+            ["results", "FEATURE_OUTPUT"],  # results.grid = output file
+            [name, "ALGORITHM_NAME"],
+            [name_op, "OUTPUT_LIST"],
+            [algo_xml, "ALGORITHM_XML"],
+            [str(sgems_files), "OBJECT_FILES"],
+        ]
 
         with open(self.template_file) as sst:
             template = sst.read()
         for i in range(len(params)):  # Replaces the parameters
             template = template.replace(params[i][1], params[i][0])
 
-        with open(self.command_name, 'w') as sstw:  # Write sgems python file
+        with open(self.command_name, "w") as sstw:  # Write sgems python file
             sstw.write(template)
 
     def script_file(self):
         """Create script file"""
-        run_script = jp(self.res_dir, 'sgems.script')
-        rscpt = open(run_script, 'w')
-        rscpt.write(' '.join(['RunScript', self.command_name]))
+        run_script = jp(self.res_dir, "sgems.script")
+        rscpt = open(run_script, "w")
+        rscpt.write(" ".join(["RunScript", self.command_name]))
         rscpt.close()
 
     def bat_file(self):
         """Create bat file"""
-        if not os.path.isfile(jp(self.res_dir, 'sgems.script')):
+        if not os.path.isfile(jp(self.res_dir, "sgems.script")):
             self.script_file()
 
-        batch = jp(self.res_dir, 'RunSgems.bat')
-        bat = open(batch, 'w')
-        bat.write(' '.join(['cd', self.res_dir, '\n']))
-        bat.write(' '.join([self.exe_name, 'sgems.script']))
+        batch = jp(self.res_dir, "RunSgems.bat")
+        bat = open(batch, "w")
+        bat.write(" ".join(["cd", self.res_dir, "\n"]))
+        bat.write(" ".join([self.exe_name, "sgems.script"]))
         bat.close()
 
     def run(self):
         """Call bat file, run sgems"""
-        batch = jp(self.res_dir, 'RunSgems.bat')
+        batch = jp(self.res_dir, "RunSgems.bat")
         if not os.path.isfile(batch):
             self.bat_file()
         start = time.time()
@@ -150,4 +176,4 @@ def run(self):
             pass
 
         subprocess.call([batch])  # Opens the BAT file
-        print('ran algorithm in {} s'.format(time.time() - start))
+        print("ran algorithm in {} s".format(time.time() - start))
diff --git a/pysgems/utils/__init__.py b/pysgems/utils/__init__.py
index 539a0ca..9c073a9 100644
--- a/pysgems/utils/__init__.py
+++ b/pysgems/utils/__init__.py
@@ -1,2 +1 @@
 #  Copyright (c) 2020. Robin Thibaut, Ghent University
-
diff --git a/setup.py b/setup.py
index 65896fc..babfe27 100644
--- a/setup.py
+++ b/setup.py
@@ -1,19 +1,21 @@
-from setuptools import setup, find_packages
+from setuptools import find_packages, setup
 
-my_pckg = find_packages(include=['pysgems'], exclude=['pysgems.deprecated', 'pysgems.develop'])
+my_pckg = find_packages(
+    include=["pysgems"], exclude=["pysgems.deprecated", "pysgems.develop"]
+)
 print(my_pckg)
 
 setup(
-    name='pysgems',
-    version='1.1.7',
+    name="pysgems",
+    version="1.1.7",
     packages=my_pckg,
     include_package_data=True,
-    url='https://github.com/robinthibaut/pysgems',
-    license='MIT',
-    author='Robin Thibaut',
-    author_email='robin.thibaut@UGent.be',
-    description='Use SGeMS (Stanford Geostatistical Modeling Software) within Python.',
-    install_requires=['numpy', 'pandas', 'scipy', 'matplotlib'],
+    url="https://github.com/robinthibaut/pysgems",
+    license="MIT",
+    author="Robin Thibaut",
+    author_email="robin.thibaut@UGent.be",
+    description="Use SGeMS (Stanford Geostatistical Modeling Software) within Python.",
+    install_requires=["numpy", "pandas", "scipy", "matplotlib"],
     classifiers=[
         "Programming Language :: Python :: 3",
         "License :: OSI Approved :: MIT License",