diff --git a/nes/__init__.py b/nes/__init__.py
index 884723bd725552252df224dc28550b1f849950ee..51ac4bcf7fd79bf7dd42ed0c07395fa722a384f7 100644
--- a/nes/__init__.py
+++ b/nes/__init__.py
@@ -2,6 +2,7 @@ __date__ = "2023-06-22"
 __version__ = "1.1.3"
 
 from .load_nes import open_netcdf, concatenate_netcdfs
+# from .load_nes import open_raster
 from .create_nes import create_nes, from_shapefile
 from .nc_projections import *
 from .methods.cell_measures import calculate_geometry_area
diff --git a/nes/create_nes.py b/nes/create_nes.py
index 7444bb5e2dfa83a5558a32de0ae9fc082abaf3a7..98f81f07a0eed533b945c88ac86db85e6e6fe752 100644
--- a/nes/create_nes.py
+++ b/nes/create_nes.py
@@ -92,32 +92,32 @@ def create_nes(comm=None, info=False, projection=None, parallel_method='Y', bala
             parallel_method = 'X'
         elif parallel_method == 'T':
             raise NotImplementedError("Parallel method T not implemented yet")
-        nessy = PointsNes(comm=comm, dataset=None, xarray=False, info=info, parallel_method=parallel_method,
+        nessy = PointsNes(comm=comm, dataset=None, info=info, parallel_method=parallel_method,
                           avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                           first_level=first_level, last_level=last_level, balanced=balanced,
                           create_nes=True, times=times, **kwargs)
     elif projection in ['regular', 'global']:
-        nessy = LatLonNes(comm=comm, dataset=None, xarray=False, info=info, parallel_method=parallel_method,
+        nessy = LatLonNes(comm=comm, dataset=None, info=info, parallel_method=parallel_method,
                           avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                           first_level=first_level, last_level=last_level, balanced=balanced,
                           create_nes=True, times=times, **kwargs)
     elif projection == 'rotated':
-        nessy = RotatedNes(comm=comm, dataset=None, xarray=False, info=info, parallel_method=parallel_method,
+        nessy = RotatedNes(comm=comm, dataset=None, info=info, parallel_method=parallel_method,
                            avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                            first_level=first_level, last_level=last_level, balanced=balanced,
                            create_nes=True, times=times, **kwargs)
     elif projection == 'rotated-nested':
-        nessy = RotatedNestedNes(comm=comm, dataset=None, xarray=False, info=info, parallel_method=parallel_method,
+        nessy = RotatedNestedNes(comm=comm, dataset=None, info=info, parallel_method=parallel_method,
                                  avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                  first_level=first_level, last_level=last_level, balanced=balanced,
                                  create_nes=True, times=times, **kwargs)
     elif projection == 'lcc':
-        nessy = LCCNes(comm=comm, dataset=None, xarray=False, info=info, parallel_method=parallel_method,
+        nessy = LCCNes(comm=comm, dataset=None, info=info, parallel_method=parallel_method,
                        avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                        first_level=first_level, last_level=last_level, balanced=balanced,
                        create_nes=True, times=times, **kwargs)
     elif projection == 'mercator':
-        nessy = MercatorNes(comm=comm, dataset=None, xarray=False, info=info, parallel_method=parallel_method,
+        nessy = MercatorNes(comm=comm, dataset=None, info=info, parallel_method=parallel_method,
                             avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                             first_level=first_level, last_level=last_level, balanced=balanced,
                             create_nes=True, times=times, **kwargs)
diff --git a/nes/load_nes.py b/nes/load_nes.py
index dd08addaae9d2a99a4723aef36e5efdb2ae5c0a8..0fb271eb2e33b150e51bb41188034d26df61f63f 100644
--- a/nes/load_nes.py
+++ b/nes/load_nes.py
@@ -12,7 +12,7 @@ DIM_VAR_NAMES = ['lat', 'latitude', 'lat_bnds', 'lon', 'longitude', 'lon_bnds',
                  'cell_area', 'crs', 'rotated_pole', 'x', 'y', 'rlat', 'rlon', 'Lambert_conformal', 'mercator']
 
 
-def open_netcdf(path, comm=None, xarray=False, info=False, parallel_method='Y', avoid_first_hours=0, avoid_last_hours=0,
+def open_netcdf(path, comm=None, info=False, parallel_method='Y', avoid_first_hours=0, avoid_last_hours=0,
                 first_level=0, last_level=None, balanced=False):
     """
     Open a netCDF file.
@@ -23,8 +23,6 @@ def open_netcdf(path, comm=None, xarray=False, info=False, parallel_method='Y',
         Path to the NetCDF file to read.
     comm : MPI.COMM
         MPI communicator to use in that netCDF. Default: MPI.COMM_WORLD.
-    xarray : bool
-        (Not working) Indicates if you want to use xarray. Default: False.
     info : bool
         Indicates if you want to print (stdout) the reading/writing steps.
     avoid_first_hours : int
@@ -54,19 +52,17 @@ def open_netcdf(path, comm=None, xarray=False, info=False, parallel_method='Y',
 
     if not os.path.exists(path):
         raise FileNotFoundError(path)
-    if xarray:
-        dataset = None
-    else:
-        dataset = Dataset(path, format="NETCDF4", mode='r', parallel=False)
-        # Parallel is not needed for reading
-        # if comm.Get_size() == 1:
-        #     dataset = Dataset(path, format="NETCDF4", mode='r', parallel=False)
-        # else:
-        #     dataset = Dataset(path, format="NETCDF4", mode='r', parallel=True, comm=comm, info=MPI.Info())
+
+    dataset = Dataset(path, format="NETCDF4", mode='r', parallel=False)
+    # Parallel is not needed for reading
+    # if comm.Get_size() == 1:
+    #     dataset = Dataset(path, format="NETCDF4", mode='r', parallel=False)
+    # else:
+    #     dataset = Dataset(path, format="NETCDF4", mode='r', parallel=True, comm=comm, info=MPI.Info())
 
     if __is_rotated(dataset):
         # Rotated grids
-        nessy = RotatedNes(comm=comm, dataset=dataset, xarray=xarray, info=info, parallel_method=parallel_method,
+        nessy = RotatedNes(comm=comm, dataset=dataset, info=info, parallel_method=parallel_method,
                            avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                            first_level=first_level, last_level=last_level, create_nes=False, balanced=balanced,)
     elif __is_points(dataset):
@@ -76,35 +72,35 @@ def open_netcdf(path, comm=None, xarray=False, info=False, parallel_method='Y',
             parallel_method = 'X'
         if __is_points_ghost(dataset):
             # Points - GHOST
-            nessy = PointsNesGHOST(comm=comm, dataset=dataset, xarray=xarray, info=info,
+            nessy = PointsNesGHOST(comm=comm, dataset=dataset, info=info,
                                    parallel_method=parallel_method,
                                    avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                    first_level=first_level, last_level=last_level, create_nes=False, balanced=balanced,)
         elif __is_points_providentia(dataset):
             # Points - Providentia
-            nessy = PointsNesProvidentia(comm=comm, dataset=dataset, xarray=xarray, info=info,
+            nessy = PointsNesProvidentia(comm=comm, dataset=dataset, info=info,
                                          parallel_method=parallel_method,
                                          avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                          first_level=first_level, last_level=last_level, create_nes=False, 
                                          balanced=balanced,)
         else:
             # Points - non-GHOST
-            nessy = PointsNes(comm=comm, dataset=dataset, xarray=xarray, info=info, parallel_method=parallel_method,
+            nessy = PointsNes(comm=comm, dataset=dataset, info=info, parallel_method=parallel_method,
                               avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                               first_level=first_level, last_level=last_level, create_nes=False, balanced=balanced,)
     elif __is_lcc(dataset):
         # Lambert conformal conic grids
-        nessy = LCCNes(comm=comm, dataset=dataset, xarray=xarray, info=info, parallel_method=parallel_method,
+        nessy = LCCNes(comm=comm, dataset=dataset, info=info, parallel_method=parallel_method,
                        avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                        first_level=first_level, last_level=last_level, create_nes=False, balanced=balanced,)
     elif __is_mercator(dataset):
         # Mercator grids
-        nessy = MercatorNes(comm=comm, dataset=dataset, xarray=xarray, info=info, parallel_method=parallel_method,
+        nessy = MercatorNes(comm=comm, dataset=dataset, info=info, parallel_method=parallel_method,
                             avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                             first_level=first_level, last_level=last_level, create_nes=False, balanced=balanced,)
     else:
         # Regular grids
-        nessy = LatLonNes(comm=comm, dataset=dataset, xarray=xarray, info=info, parallel_method=parallel_method,
+        nessy = LatLonNes(comm=comm, dataset=dataset, info=info, parallel_method=parallel_method,
                           avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                           first_level=first_level, last_level=last_level, create_nes=False, balanced=balanced,)
 
@@ -326,3 +322,34 @@ def concatenate_netcdfs(nessy_list, comm=None, info=False, parallel_method='Y',
             nessy_first.concatenate(aux_nessy)
         
     return nessy_first
+
+
+def open_raster(path, comm=None, info=False):
+    """
+    Open a Raster (TIF) file
+
+    Parameters
+    ----------
+    path : str
+        Path to the Raster (TIFF) file to read.
+    comm : MPI.COMM
+        MPI communicator to use in that netCDF. Default: MPI.COMM_WORLD.
+    info : bool
+        Indicates if you want to print (stdout) the reading/writing steps.
+
+    Returns
+    -------
+    RasterNes
+        Raster NES object
+    """
+    if comm is None:
+        comm = MPI.COMM_WORLD
+    else:
+        comm = comm
+
+    if not os.path.exists(path):
+        raise FileNotFoundError(path)
+
+    nessy = RasterNes(path=path, comm=comm, info=info)
+
+    return nessy
diff --git a/nes/methods/horizontal_interpolation.py b/nes/methods/horizontal_interpolation.py
index 3887bd905b1dc47b8fac119a31fb9a331484b56e..aacd498fc8bdd48e2aed2911a78e687a8d1dcbb7 100644
--- a/nes/methods/horizontal_interpolation.py
+++ b/nes/methods/horizontal_interpolation.py
@@ -84,7 +84,7 @@ def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='Neares
 
     # Remove original file information
     final_dst.__ini_path = None
-    final_dst.netcdf = None
+    final_dst.dataset = None
     final_dst.dataset = None
 
     # Return final_dst
diff --git a/nes/methods/vertical_interpolation.py b/nes/methods/vertical_interpolation.py
index 7260fdbe191484129649b0feff7e3efcc4d95728..3100cf93106efe913a696dfd15c41376b11ee083 100644
--- a/nes/methods/vertical_interpolation.py
+++ b/nes/methods/vertical_interpolation.py
@@ -328,6 +328,5 @@ def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear',
     # Remove original file information
     self.__ini_path = None
     self.dataset = None
-    self.netcdf = None
 
     return self
diff --git a/nes/nc_projections/__init__.py b/nes/nc_projections/__init__.py
index d4c4b9f8656b43a1646228e46629ade913236006..fa530f9f0e1985c67066d46225d4f06f96545697 100644
--- a/nes/nc_projections/__init__.py
+++ b/nes/nc_projections/__init__.py
@@ -7,3 +7,4 @@ from .points_nes_ghost import PointsNesGHOST
 from .points_nes_providentia import PointsNesProvidentia
 from .lcc_nes import LCCNes
 from .mercator_nes import MercatorNes
+# from .raster_nes import RasterNes
diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py
index 5c2fcd0350e1c05a4af589ad571891270d3e4143..f8640e14492e279517412b28f76af50816c90b2b 100644
--- a/nes/nc_projections/default_nes.py
+++ b/nes/nc_projections/default_nes.py
@@ -6,7 +6,6 @@ import warnings
 import numpy as np
 import pandas as pd
 from datetime import timedelta
-from xarray import open_dataset
 from netCDF4 import Dataset, num2date, date2num, stringtochar
 from mpi4py import MPI
 from cfunits import Units
@@ -37,17 +36,13 @@ class Nes(object):
         Size of the communicator.
     info : bool
         Indicates if you want to print reading/writing info.
-    is_xarray : bool
-        (Not working) Indicates if you want to use xarray as default.
     __ini_path : str
         Path to the original file to read when open_netcdf is called.
     hours_start : int
         Number of hours to avoid from the first original values.
     hours_end : int
         Number of hours to avoid from the last original values.
-    dataset : xr.Dataset
-        (not working) xArray Dataset.
-    netcdf : Dataset
+    dataset : Dataset
         netcdf4-python Dataset.
     variables : dict
         Variables information.
@@ -101,7 +96,7 @@ class Nes(object):
         Dictionary with the projection information.
 
     """
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='Y',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
                  balanced=False, times=None, **kwargs):
         """
@@ -115,10 +110,8 @@ class Nes(object):
             Path to the NetCDF to initialize the object.
         info: bool
             Indicates if you want to get reading/writing info.
-        dataset: Dataset
+        dataset: Dataset or None
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default over Y axis
             accepted values: ['X', 'Y', 'T'].
@@ -150,7 +143,6 @@ class Nes(object):
 
         # General info
         self.info = info
-        self.is_xarray = xarray
         self.__ini_path = path
         self.shapefile = None
 
@@ -180,7 +172,6 @@ class Nes(object):
         # NetCDF object
         if create_nes:
 
-            self.netcdf = None
             self.dataset = None
 
             # Set string length
@@ -218,21 +209,10 @@ class Nes(object):
             self.global_attrs = self.__get_global_attributes(create_nes)
 
         else:
-
             if dataset is not None:
-                if self.is_xarray:
-                    self.dataset = dataset
-                    self.netcdf = None
-                else:
-                    self.dataset = None
-                    self.netcdf = dataset
+                self.dataset = dataset
             elif self.__ini_path is not None:
-                if self.is_xarray:
-                    self.dataset = self.__open_dataset()
-                    self.netcdf = None
-                else:
-                    self.dataset = None
-                    self.netcdf = self.__open_netcdf4()
+                self.open()
 
             # Get string length
             self.strlen = self._get_strlen()
@@ -296,7 +276,7 @@ class Nes(object):
         self.first_level = None
 
     @staticmethod
-    def new(comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y', 
+    def new(comm=None, path=None, info=False, dataset=None, parallel_method='Y',
             avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False, 
             balanced=False, times=None, **kwargs):
         """
@@ -312,8 +292,6 @@ class Nes(object):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default over Y axis
             accepted values: ['X', 'Y', 'T'].
@@ -334,7 +312,7 @@ class Nes(object):
             List of times to substitute the current ones while creation.
         """
 
-        new = Nes(comm=comm, path=path, info=info, dataset=dataset, xarray=xarray, parallel_method=parallel_method,
+        new = Nes(comm=comm, path=path, info=info, dataset=dataset, parallel_method=parallel_method,
                   avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours, first_level=first_level,
                   last_level=last_level, create_nes=create_nes, balanced=balanced, times=times, **kwargs)
 
@@ -350,8 +328,8 @@ class Nes(object):
             Max length of the string data
         """  
         
-        if 'strlen' in self.netcdf.dimensions:
-            strlen = self.netcdf.dimensions['strlen'].size
+        if 'strlen' in self.dataset.dimensions:
+            strlen = self.dataset.dimensions['strlen'].size
         else:
             return None
         
@@ -421,7 +399,7 @@ class Nes(object):
         """
 
         d = self.__dict__
-        state = {k: d[k] for k in d if k not in ['comm', 'variables', 'netcdf', 'cell_measures']}
+        state = {k: d[k] for k in d if k not in ['comm', 'variables', 'dataset', 'cell_measures']}
 
         return state
 
@@ -485,7 +463,7 @@ class Nes(object):
         """
 
         nessy = deepcopy(self)
-        nessy.netcdf = None
+        nessy.dataset = None
         if copy_vars:
             nessy.set_communicator(self.comm)
             nessy.variables = deepcopy(self.variables)
@@ -755,16 +733,12 @@ class Nes(object):
         if isinstance(var_list, str):
             var_list = [var_list]
 
-        if self.is_xarray:
-            self.dataset = self.dataset.drop_vars(var_list)
-            self.variables = self._get_lazy_variables()
-        else:
-            if self.variables is not None:
-                for var_name in var_list:
-                    if var_name in self.variables:
-                        if 'data' in self.variables[var_name].keys():
-                            del self.variables[var_name]['data']
-                        del self.variables[var_name]
+        if self.variables is not None:
+            for var_name in var_list:
+                if var_name in self.variables:
+                    if 'data' in self.variables[var_name].keys():
+                        del self.variables[var_name]['data']
+                    del self.variables[var_name]
         gc.collect()
 
         return None
@@ -1260,7 +1234,7 @@ class Nes(object):
             # Load variables if they have not been loaded previously
             if self.variables[var_name]['data'] is None:
                 self.load(var_name)
-            
+
             # Get original file shape
             nessy_shape = self.variables[var_name]['data'].shape
 
@@ -1272,7 +1246,7 @@ class Nes(object):
             for curr_time in self.time:
                 # Get previous time given a set of hours
                 prev_time = curr_time - timedelta(hours=(hours-1))
-                
+
                 # Get time indices
                 curr_time_id = self.find_time_id(curr_time)
                 prev_time_id = self.find_time_id(prev_time)
@@ -1284,7 +1258,7 @@ class Nes(object):
                     aux_nessy.variables[var_name]['data'][curr_time_id, :, :, :] = self.variables[var_name]['data'][
                         prev_time_id:curr_time_id, :, :, :].mean(axis=0, keepdims=True)
                 # Fill with nan if previous time is not available
-                else: 
+                else:
                     if self.info:
                         msg = f'Mean between {prev_time} and {curr_time} cannot be calculated '
                         msg += f'because data for {prev_time} is not available.'
@@ -1553,37 +1527,10 @@ class Nes(object):
         Open the NetCDF.
         """
 
-        if self.is_xarray:
-            self.dataset = self.__open_dataset()
-            self.netcdf = None
-        else:
-            self.dataset = None
-            self.netcdf = self.__open_netcdf4()
+        self.dataset = self.__open_netcdf4()
 
         return None
 
-    def __open_dataset(self):
-        """
-        Open the NetCDF with xarray.
-
-        Returns
-        -------
-        dataset : xr.Dataset
-            Open dataset.
-        """
-
-        if self.master:
-            warnings.filterwarnings('ignore')  # Disabling warnings while reading MONARCH original file
-            dataset = open_dataset(self.__ini_path, decode_coords='all')
-            warnings.filterwarnings('default')  # Re-activating warnings
-        else:
-            dataset = None
-
-        dataset = self.comm.bcast(dataset, root=0)
-        self.dataset = dataset
-
-        return dataset
-
     def __open_netcdf4(self, mode='r'):
         """
         Open the NetCDF with netcdf4-python.
@@ -1604,7 +1551,7 @@ class Nes(object):
         else:
             netcdf = Dataset(self.__ini_path, format="NETCDF4", mode=mode, parallel=True, comm=self.comm,
                              info=MPI.Info())
-        self.netcdf = netcdf
+        self.dataset = netcdf
 
         return netcdf
 
@@ -1616,9 +1563,9 @@ class Nes(object):
             if self.master:
                 self.serial_nc.close()
             self.serial_nc = None
-        if (hasattr(self, 'netcdf')) and (self.netcdf is not None):
-            self.netcdf.close()
-            self.netcdf = None
+        if (hasattr(self, 'netcdf')) and (self.dataset is not None):
+            self.dataset.close()
+            self.dataset = None
 
         return None
 
@@ -1756,23 +1703,20 @@ class Nes(object):
             List of times (datetime.datetime) of the NetCDF data.
         """
 
-        if self.is_xarray:
-            time = self.variables['time']
+        if self.master:
+            nc_var = self.dataset.variables['time']
+            time_data, units, calendar = self.__parse_time(nc_var)
+            # Extracting time resolution depending on the units
+            self._time_resolution = self.__get_time_resolution_from_units(units)
+            # Checking if it is a climatology dataset
+            if hasattr(nc_var, 'climatology'):
+                self._climatology = True
+                self._climatology_var_name = nc_var.climatology
+            time = num2date(time_data, units, calendar=calendar)
+            time = [aux.replace(second=0, microsecond=0) for aux in time]
         else:
-            if self.master:
-                nc_var = self.netcdf.variables['time']
-                time_data, units, calendar = self.__parse_time(nc_var)
-                # Extracting time resolution depending on the units
-                self._time_resolution = self.__get_time_resolution_from_units(units)
-                # Checking if it is a climatology dataset
-                if hasattr(nc_var, 'climatology'):
-                    self._climatology = True
-                    self._climatology_var_name = nc_var.climatology
-                time = num2date(time_data, units, calendar=calendar)
-                time = [aux.replace(second=0, microsecond=0) for aux in time]
-            else:
-                time = None
-            time = self.comm.bcast(time, root=0)
+            time = None
+        time = self.comm.bcast(time, root=0)
         self.free_vars('time')
 
         return time
@@ -1794,12 +1738,12 @@ class Nes(object):
 
         if self.master:
             if not create_nes:
-                if 'time_bnds' in self.netcdf.variables.keys() or self._climatology:
-                    time = self.netcdf.variables['time']
+                if 'time_bnds' in self.dataset.variables.keys() or self._climatology:
+                    time = self.dataset.variables['time']
                     if self._climatology:
-                        nc_var = self.netcdf.variables[self._climatology_var_name]
+                        nc_var = self.dataset.variables[self._climatology_var_name]
                     else:
-                        nc_var = self.netcdf.variables['time_bnds']
+                        nc_var = self.dataset.variables['time_bnds']
                     time_bnds = num2date(nc_var[:], self.__parse_time_unit(time.units),
                                          calendar=time.calendar).tolist()
                 else:
@@ -1832,28 +1776,24 @@ class Nes(object):
             Longitude bounds of the NetCDF data.
         """
 
-        if self.is_xarray:
-            lat_bnds = self.variables['lat_bnds']
-            lon_bnds = self.variables['lon_bnds']
-        else:
-            if self.master:
-                if not create_nes:
-                    if 'lat_bnds' in self.netcdf.variables.keys():
-                        lat_bnds = {'data': self._unmask_array(self.netcdf.variables['lat_bnds'][:])}
-                    else:
-                        lat_bnds = None
-                    if 'lon_bnds' in self.netcdf.variables.keys():
-                        lon_bnds = {'data': self._unmask_array(self.netcdf.variables['lon_bnds'][:])}
-                    else:
-                        lon_bnds = None
+        if self.master:
+            if not create_nes:
+                if 'lat_bnds' in self.dataset.variables.keys():
+                    lat_bnds = {'data': self._unmask_array(self.dataset.variables['lat_bnds'][:])}
                 else:
                     lat_bnds = None
+                if 'lon_bnds' in self.dataset.variables.keys():
+                    lon_bnds = {'data': self._unmask_array(self.dataset.variables['lon_bnds'][:])}
+                else:
                     lon_bnds = None
             else:
                 lat_bnds = None
                 lon_bnds = None
-            lat_bnds = self.comm.bcast(lat_bnds, root=0)
-            lon_bnds = self.comm.bcast(lon_bnds, root=0)
+        else:
+            lat_bnds = None
+            lon_bnds = None
+        lat_bnds = self.comm.bcast(lat_bnds, root=0)
+        lon_bnds = self.comm.bcast(lon_bnds, root=0)
 
         self.free_vars(['lat_bnds', 'lon_bnds'])
 
@@ -1877,9 +1817,9 @@ class Nes(object):
         c_measures = {}
         if self.master:
             if not create_nes:
-                if 'cell_area' in self.netcdf.variables.keys():
+                if 'cell_area' in self.dataset.variables.keys():
                     c_measures['cell_area'] = {}
-                    c_measures['cell_area']['data'] = self._unmask_array(self.netcdf.variables['cell_area'][:])
+                    c_measures['cell_area']['data'] = self._unmask_array(self.dataset.variables['cell_area'][:])
         c_measures = self.comm.bcast(c_measures, root=0)
 
         self.free_vars(['cell_area'])
@@ -1908,12 +1848,10 @@ class Nes(object):
 
         try:
             dimension_name = set(possible_names).intersection(set(self.variables.keys())).pop()
-            if self.is_xarray:
-                nc_var = self.dataset[dimension_name]
-            else:
-                nc_var = self.variables[dimension_name].copy()
-                nc_var['data'] = self.netcdf.variables[dimension_name][:]
-                if hasattr(nc_var, 'units'):
+
+            nc_var = self.variables[dimension_name].copy()
+            nc_var['data'] = self.dataset.variables[dimension_name][:]
+            if hasattr(nc_var, 'units'):
                     if nc_var['units'] in ['unitless', '-']:
                         nc_var['units'] = ''
             self.free_vars(dimension_name)
@@ -2034,32 +1972,29 @@ class Nes(object):
              ...}
         """
 
-        if self.is_xarray:
-            variables = self.dataset.variables
+        if self.master:
+            variables = {}
+            # Initialise data
+            for var_name, var_info in self.dataset.variables.items():
+                variables[var_name] = {}
+                variables[var_name]['data'] = None
+                variables[var_name]['dimensions'] = var_info.dimensions
+                variables[var_name]['dtype'] = var_info.dtype
+                if variables[var_name]['dtype'] in [str, np.object]:
+                    if self.strlen is None:
+                        self.set_strlen()
+                    variables[var_name]['dtype'] = str
+
+                # Avoid some attributes
+                for attrname in var_info.ncattrs():
+                    if attrname not in ['missing_value', '_FillValue']:
+                        value = getattr(var_info, attrname)
+                        if value in ['unitless', '-']:
+                            value = ''
+                        variables[var_name][attrname] = value
         else:
-            if self.master:
-                variables = {}
-                # Initialise data
-                for var_name, var_info in self.netcdf.variables.items():
-                    variables[var_name] = {}
-                    variables[var_name]['data'] = None
-                    variables[var_name]['dimensions'] = var_info.dimensions
-                    variables[var_name]['dtype'] = var_info.dtype
-                    if variables[var_name]['dtype'] in [str, np.object]:
-                        if self.strlen is None:
-                            self.set_strlen()
-                        variables[var_name]['dtype'] = str
-
-                    # Avoid some attributes
-                    for attrname in var_info.ncattrs():
-                        if attrname not in ['missing_value', '_FillValue']:
-                            value = getattr(var_info, attrname)
-                            if value in ['unitless', '-']:
-                                value = ''
-                            variables[var_name][attrname] = value
-            else:
-                variables = None
-            variables = self.comm.bcast(variables, root=0)
+            variables = None
+        variables = self.comm.bcast(variables, root=0)
 
         return variables
 
@@ -2078,7 +2013,7 @@ class Nes(object):
             Portion of the variable data corresponding to the rank.
         """
 
-        nc_var = self.netcdf.variables[var_name]
+        nc_var = self.dataset.variables[var_name]
         var_dims = nc_var.dimensions
 
         # Read data in 4 dimensions
@@ -2153,10 +2088,10 @@ class Nes(object):
             List (or single string) of the variables to be loaded.
         """
         
-        if (self.__ini_path is None) and (self.dataset is None) and (self.netcdf is None):
+        if (self.__ini_path is None) and (self.dataset is None):
             raise RuntimeError('Only data from existing files can be loaded.')
 
-        if self.netcdf is None:
+        if self.dataset is None:
             self.__open_dataset()
             close = True
         else:
@@ -2250,7 +2185,6 @@ class Nes(object):
 
         if isinstance(aux_nessy, str):
             aux_nessy = self.new(path=aux_nessy, comm=self.comm, parallel_method=self.parallel_method,
-                                 xarray=self.is_xarray,
                                  avoid_first_hours=self.hours_start, avoid_last_hours=self.hours_end,
                                  first_level=self.first_level, last_level=self.last_level)
             new = True
@@ -2288,12 +2222,10 @@ class Nes(object):
         """
 
         gl_attrs = {}
-        if self.is_xarray:
-            gl_attrs = self.dataset.attrs
-        else:
-            if not create_nes:
-                for attrname in self.netcdf.ncattrs():
-                    gl_attrs[attrname] = getattr(self.netcdf, attrname)
+
+        if not create_nes:
+            for attrname in self.dataset.ncattrs():
+                gl_attrs[attrname] = getattr(self.dataset, attrname)
 
         return gl_attrs
 
@@ -2865,7 +2797,7 @@ class Nes(object):
                         if att_value is not None:
                             if self.info:
                                 print("Rank {0:03d}: Filling {1})".format(self.rank, var_name))
-                            var = self.netcdf.variables[var_name]
+                            var = self.dataset.variables[var_name]
                             if isinstance(att_value, int) and att_value == 0:
                                 var[i_time,
                                     self.write_axis_limits['z_min']:self.write_axis_limits['z_max'],
@@ -2976,7 +2908,7 @@ class Nes(object):
         netcdf.setncattr('Conventions', 'CF-1.7')
 
         if keep_open:
-            self.netcdf = netcdf
+            self.dataset = netcdf
         else:
             netcdf.close()
 
@@ -2985,7 +2917,6 @@ class Nes(object):
     def __to_netcdf_cams_ra(self, path):
         return to_netcdf_cams_ra(self, path)
 
-
     def to_netcdf(self, path, compression_level=0, serial=False, info=False, chunking=False, type='NES',
                   keep_open=False):
         """
@@ -3013,56 +2944,54 @@ class Nes(object):
         self.info = info
         self.serial_nc = None
         self.zip_lvl = compression_level
-        if self.is_xarray:
-            raise NotImplementedError("Writing with xarray not implemented")
-        else:
-            # if serial:
-            if serial and self.size > 1:
-                try:
-                    data = self._gather_data(self.variables)
-                except KeyError:
-                    data = self.__gather_data_py_object(self.variables)
-                try:
-                    c_measures = self._gather_data(self.cell_measures)
-                except KeyError:
-                    c_measures = self.__gather_data_py_object(self.cell_measures)
-                if self.master:
-                    new_nc = self.copy(copy_vars=False)
-                    new_nc.set_communicator(MPI.COMM_SELF)
-                    new_nc.variables = data
-                    new_nc.cell_measures = c_measures
-                    if type in ['NES', 'DEFAULT']:
-                        new_nc.__to_netcdf_py(path, keep_open=keep_open)
-                    elif type == 'CAMS_RA':
-                        new_nc.__to_netcdf_cams_ra(path)
-                    elif type == 'MONARCH':
-                        to_netcdf_monarch(new_nc, path, chunking=chunking, keep_open=keep_open)
-                    elif type == 'CMAQ':
-                        to_netcdf_cmaq(new_nc, path, keep_open=keep_open)
-                    elif type == 'WRF_CHEM':
-                        to_netcdf_wrf_chem(new_nc, path, keep_open=keep_open)
-                    else:
-                        msg = "Unknown NetCDF type '{0}'. ".format(nc_type)
-                        msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)"
-                        raise ValueError(msg)
-                    self.serial_nc = new_nc
-                else:
-                    self.serial_nc = True
-            else:
+
+        # if serial:
+        if serial and self.size > 1:
+            try:
+                data = self._gather_data(self.variables)
+            except KeyError:
+                data = self.__gather_data_py_object(self.variables)
+            try:
+                c_measures = self._gather_data(self.cell_measures)
+            except KeyError:
+                c_measures = self.__gather_data_py_object(self.cell_measures)
+            if self.master:
+                new_nc = self.copy(copy_vars=False)
+                new_nc.set_communicator(MPI.COMM_SELF)
+                new_nc.variables = data
+                new_nc.cell_measures = c_measures
                 if type in ['NES', 'DEFAULT']:
-                    self.__to_netcdf_py(path, chunking=chunking, keep_open=keep_open)
-                elif nc_type == 'CAMS_RA':
-                    self.__to_netcdf_cams_ra(path)
-                elif nc_type == 'MONARCH':
-                    to_netcdf_monarch(self, path, chunking=chunking, keep_open=keep_open)
-                elif nc_type == 'CMAQ':
-                    to_netcdf_cmaq(self, path, keep_open=keep_open)
-                elif nc_type == 'WRF_CHEM':
-                    to_netcdf_wrf_chem(self, path, keep_open=keep_open)
+                    new_nc.__to_netcdf_py(path, keep_open=keep_open)
+                elif type == 'CAMS_RA':
+                    new_nc.__to_netcdf_cams_ra(path)
+                elif type == 'MONARCH':
+                    to_netcdf_monarch(new_nc, path, chunking=chunking, keep_open=keep_open)
+                elif type == 'CMAQ':
+                    to_netcdf_cmaq(new_nc, path, keep_open=keep_open)
+                elif type == 'WRF_CHEM':
+                    to_netcdf_wrf_chem(new_nc, path, keep_open=keep_open)
                 else:
                     msg = "Unknown NetCDF type '{0}'. ".format(nc_type)
                     msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)"
                     raise ValueError(msg)
+                self.serial_nc = new_nc
+            else:
+                self.serial_nc = True
+        else:
+            if type in ['NES', 'DEFAULT']:
+                self.__to_netcdf_py(path, chunking=chunking, keep_open=keep_open)
+            elif nc_type == 'CAMS_RA':
+                self.__to_netcdf_cams_ra(path)
+            elif nc_type == 'MONARCH':
+                to_netcdf_monarch(self, path, chunking=chunking, keep_open=keep_open)
+            elif nc_type == 'CMAQ':
+                to_netcdf_cmaq(self, path, keep_open=keep_open)
+            elif nc_type == 'WRF_CHEM':
+                to_netcdf_wrf_chem(self, path, keep_open=keep_open)
+            else:
+                msg = "Unknown NetCDF type '{0}'. ".format(nc_type)
+                msg += "Use CAMS_RA, MONARCH or NES (or DEFAULT)"
+                raise ValueError(msg)
 
         self.info = old_info
 
diff --git a/nes/nc_projections/latlon_nes.py b/nes/nc_projections/latlon_nes.py
index 01516d4c1c10301aeb94d403ef3e1f1f2262e415..e42eb610f265a07a236f15426c9224d994f74761 100644
--- a/nes/nc_projections/latlon_nes.py
+++ b/nes/nc_projections/latlon_nes.py
@@ -20,7 +20,7 @@ class LatLonNes(Nes):
         Tuple with the name of the dimensions of the Longitude values.
         ('lon',) for a regular latitude-longitude projection.
     """
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='Y',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
                  balanced=False, times=None, **kwargs):
         """
@@ -36,8 +36,6 @@ class LatLonNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -59,7 +57,7 @@ class LatLonNes(Nes):
         """
         
         super(LatLonNes, self).__init__(comm=comm, path=path, info=info, dataset=dataset,
-                                        xarray=xarray, parallel_method=parallel_method, balanced=balanced,
+                                        parallel_method=parallel_method, balanced=balanced,
                                         avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                         first_level=first_level, last_level=last_level, create_nes=create_nes,
                                         times=times, **kwargs)
@@ -79,7 +77,7 @@ class LatLonNes(Nes):
         self.free_vars('crs')
 
     @staticmethod
-    def new(comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y', 
+    def new(comm=None, path=None, info=False, dataset=None, parallel_method='Y',
             avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False, 
             balanced=False, times=None, **kwargs):
         """
@@ -95,8 +93,6 @@ class LatLonNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -117,7 +113,7 @@ class LatLonNes(Nes):
             List of times to substitute the current ones while creation.
         """
 
-        new = LatLonNes(comm=comm, path=path, info=info, dataset=dataset, xarray=xarray, 
+        new = LatLonNes(comm=comm, path=path, info=info, dataset=dataset,
                         parallel_method=parallel_method, avoid_first_hours=avoid_first_hours, 
                         avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level,
                         create_nes=create_nes, balanced=balanced, times=times, **kwargs)
diff --git a/nes/nc_projections/lcc_nes.py b/nes/nc_projections/lcc_nes.py
index 8225c2d6bdefc7fbbe889858d4f7af7a9e85ee44..c225d5b769956068cee4c84d8e33bdc9cf24cd4b 100644
--- a/nes/nc_projections/lcc_nes.py
+++ b/nes/nc_projections/lcc_nes.py
@@ -35,7 +35,7 @@ class LCCNes(Nes):
         Tuple with the name of the dimensions of the Longitude values.
         ('y', 'x') for a LCC projection.
     """
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='Y',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
                  balanced=False, times=None, **kwargs):
         """
@@ -51,8 +51,6 @@ class LCCNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -74,7 +72,7 @@ class LCCNes(Nes):
         """
 
         super(LCCNes, self).__init__(comm=comm, path=path, info=info, dataset=dataset,
-                                     xarray=xarray, parallel_method=parallel_method, balanced=balanced,
+                                     parallel_method=parallel_method, balanced=balanced,
                                      avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                      first_level=first_level, last_level=last_level, create_nes=create_nes,
                                      times=times, **kwargs)
@@ -102,7 +100,7 @@ class LCCNes(Nes):
         self.free_vars('crs')
 
     @staticmethod
-    def new(comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y', 
+    def new(comm=None, path=None, info=False, dataset=None, parallel_method='Y',
             avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False, 
             balanced=False, times=None, **kwargs):
         """
@@ -118,8 +116,6 @@ class LCCNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -140,7 +136,7 @@ class LCCNes(Nes):
             List of times to substitute the current ones while creation.
         """
 
-        new = LCCNes(comm=comm, path=path, info=info, dataset=dataset, xarray=xarray,
+        new = LCCNes(comm=comm, path=path, info=info, dataset=dataset,
                      parallel_method=parallel_method, avoid_first_hours=avoid_first_hours,
                      avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level,
                      create_nes=create_nes, balanced=balanced, times=times, **kwargs)
diff --git a/nes/nc_projections/mercator_nes.py b/nes/nc_projections/mercator_nes.py
index cc02a302b3fe8ad58c72ff94ca640585cdd70c0b..3cd46182389d258bd853518e9e67b08f707a0f95 100644
--- a/nes/nc_projections/mercator_nes.py
+++ b/nes/nc_projections/mercator_nes.py
@@ -35,7 +35,7 @@ class MercatorNes(Nes):
         Tuple with the name of the dimensions of the Longitude values.
         ('y', 'x') for a Mercator projection.
     """
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='Y',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
                  balanced=False, times=None, **kwargs):
         """
@@ -51,8 +51,6 @@ class MercatorNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -75,7 +73,7 @@ class MercatorNes(Nes):
         """
 
         super(MercatorNes, self).__init__(comm=comm, path=path, info=info, dataset=dataset,
-                                          xarray=xarray, parallel_method=parallel_method, balanced=balanced,
+                                          parallel_method=parallel_method, balanced=balanced,
                                           avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                           first_level=first_level, last_level=last_level, create_nes=create_nes,
                                           times=times, **kwargs)
@@ -103,7 +101,7 @@ class MercatorNes(Nes):
         self.free_vars('crs')
 
     @staticmethod
-    def new(comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y', 
+    def new(comm=None, path=None, info=False, dataset=None, parallel_method='Y',
             avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
             balanced=False, times=None, **kwargs):
         """
@@ -119,8 +117,6 @@ class MercatorNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -141,7 +137,7 @@ class MercatorNes(Nes):
             List of times to substitute the current ones while creation.
         """
 
-        new = MercatorNes(comm=comm, path=path, info=info, dataset=dataset, xarray=xarray, 
+        new = MercatorNes(comm=comm, path=path, info=info, dataset=dataset,
                           parallel_method=parallel_method, avoid_first_hours=avoid_first_hours,
                           avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level,
                           create_nes=create_nes, balanced=balanced, times=times, **kwargs)
diff --git a/nes/nc_projections/points_nes.py b/nes/nc_projections/points_nes.py
index ee541842d81dbd9a081ade61323c10bc503d2885..9df46e2e86e9881dbf0ebee03a60c9720289e81d 100644
--- a/nes/nc_projections/points_nes.py
+++ b/nes/nc_projections/points_nes.py
@@ -28,7 +28,7 @@ class PointsNes(Nes):
         Tuple with the name of the dimensions of the station values.
         ('station',) for a points grid.
     """
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='X',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='X',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False, 
                  balanced=False, times=None, **kwargs):
         """
@@ -44,8 +44,6 @@ class PointsNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset, None
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'X'.
             accepted values: ['X', 'T'].
@@ -69,7 +67,7 @@ class PointsNes(Nes):
         """
 
         super(PointsNes, self).__init__(comm=comm, path=path, info=info, dataset=dataset,
-                                        xarray=xarray, parallel_method=parallel_method,
+                                        parallel_method=parallel_method,
                                         avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                         first_level=first_level, last_level=last_level, create_nes=create_nes,
                                         times=times, **kwargs)
@@ -93,7 +91,7 @@ class PointsNes(Nes):
         self._lon_dim = ('station',)
 
     @staticmethod
-    def new(comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='X', 
+    def new(comm=None, path=None, info=False, dataset=None, parallel_method='X',
             avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None,
             create_nes=False, balanced=False, times=None, **kwargs):
         """
@@ -109,8 +107,6 @@ class PointsNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'X'.
             accepted values: ['X', 'T'].
@@ -131,7 +127,7 @@ class PointsNes(Nes):
             List of times to substitute the current ones while creation.
         """
 
-        new = PointsNes(comm=comm, path=path, info=info, dataset=dataset, xarray=xarray,
+        new = PointsNes(comm=comm, path=path, info=info, dataset=dataset,
                         parallel_method=parallel_method, avoid_first_hours=avoid_first_hours,
                         avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level,
                         create_nes=create_nes, balanced=balanced, times=times, **kwargs)
@@ -314,7 +310,7 @@ class PointsNes(Nes):
             Portion of the variable data corresponding to the rank.
         """
 
-        nc_var = self.netcdf.variables[var_name]
+        nc_var = self.dataset.variables[var_name]
         var_dims = nc_var.dimensions
 
         # Read data in 1 or 2 dimensions
diff --git a/nes/nc_projections/points_nes_ghost.py b/nes/nc_projections/points_nes_ghost.py
index 09c5a443d1daa6fd7dc1d9d7ef0bf3d9f5733341..b528edf530cfa3dce984942e2480a874933086e8 100644
--- a/nes/nc_projections/points_nes_ghost.py
+++ b/nes/nc_projections/points_nes_ghost.py
@@ -25,7 +25,7 @@ class PointsNesGHOST(PointsNes):
         Data flags (given by data provider) dictionary with the portion of 'data' corresponding to the rank values.
     """
 
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='X',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='X',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
                  balanced=False, times=None, **kwargs):
         """
@@ -41,8 +41,6 @@ class PointsNesGHOST(PointsNes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'X'.
             Accepted values: ['X'].
@@ -64,7 +62,7 @@ class PointsNesGHOST(PointsNes):
         """
 
         super(PointsNesGHOST, self).__init__(comm=comm, path=path, info=info, dataset=dataset,
-                                             xarray=xarray, parallel_method=parallel_method,
+                                             parallel_method=parallel_method,
                                              avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                              first_level=first_level, last_level=last_level, create_nes=create_nes,
                                              times=times, **kwargs)
@@ -78,7 +76,7 @@ class PointsNesGHOST(PointsNes):
         self.qa = self._get_coordinate_values(self._qa, 'X')
 
     @staticmethod
-    def new(comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='X', 
+    def new(comm=None, path=None, info=False, dataset=None, parallel_method='X',
             avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False, 
             balanced=False, times=None, **kwargs):
         """
@@ -94,8 +92,6 @@ class PointsNesGHOST(PointsNes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'X'.
             Accepted values: ['X'].
@@ -116,7 +112,7 @@ class PointsNesGHOST(PointsNes):
             List of times to substitute the current ones while creation.
         """
 
-        new = PointsNesGHOST(comm=comm, path=path, info=info, dataset=dataset, xarray=xarray, 
+        new = PointsNesGHOST(comm=comm, path=path, info=info, dataset=dataset,
                              parallel_method=parallel_method, avoid_first_hours=avoid_first_hours,
                              avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level,
                              create_nes=create_nes, balanced=balanced, times=times, **kwargs)
@@ -283,7 +279,7 @@ class PointsNesGHOST(PointsNes):
             Portion of the variable data corresponding to the rank.
         """
 
-        nc_var = self.netcdf.variables[var_name]
+        nc_var = self.dataset.variables[var_name]
         var_dims = nc_var.dimensions
 
         # Read data in 1 or 2 dimensions
diff --git a/nes/nc_projections/points_nes_providentia.py b/nes/nc_projections/points_nes_providentia.py
index 533e4a86bc069b5af48cec87cfa8370bfd9f6e04..4a72bc0e0749c183d904882f23021113fd1e5fa6 100644
--- a/nes/nc_projections/points_nes_providentia.py
+++ b/nes/nc_projections/points_nes_providentia.py
@@ -34,9 +34,10 @@ class PointsNesProvidentia(PointsNes):
     grid_edge_lat : dict 
         Grid edge latitudes dictionary with the portion of 'data' corresponding to the rank values.
     """
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='X',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='X',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
-                 balanced=False, times=None, model_centre_lon=None, model_centre_lat=None, grid_edge_lon=None, grid_edge_lat=None,
+                 balanced=False, times=None, model_centre_lon=None, model_centre_lat=None, grid_edge_lon=None,
+                 grid_edge_lat=None,
                  **kwargs):
         """
         Initialize the PointsNesProvidentia class
@@ -51,8 +52,6 @@ class PointsNesProvidentia(PointsNes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'X'.
             Accepted values: ['X'].
@@ -82,7 +81,7 @@ class PointsNesProvidentia(PointsNes):
         """
 
         super(PointsNesProvidentia, self).__init__(comm=comm, path=path, info=info, dataset=dataset,
-                                                   xarray=xarray, parallel_method=parallel_method,
+                                                   parallel_method=parallel_method,
                                                    avoid_first_hours=avoid_first_hours, 
                                                    avoid_last_hours=avoid_last_hours,
                                                    first_level=first_level, last_level=last_level, 
@@ -108,7 +107,7 @@ class PointsNesProvidentia(PointsNes):
         self.grid_edge_lat = self._get_coordinate_values(self._grid_edge_lat, '') 
 
     @staticmethod
-    def new(comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='X', 
+    def new(comm=None, path=None, info=False, dataset=None, parallel_method='X',
             avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None,
             create_nes=False, balanced=False, times=None,
             model_centre_lon=None, model_centre_lat=None, grid_edge_lon=None, grid_edge_lat=None,
@@ -126,8 +125,6 @@ class PointsNesProvidentia(PointsNes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'X'.
             Accepted values: ['X'].
@@ -156,7 +153,7 @@ class PointsNesProvidentia(PointsNes):
             Grid edge latitudes dictionary with the portion of 'data' corresponding to the rank values.
         """
 
-        new = PointsNesProvidentia(comm=comm, path=path, info=info, dataset=dataset, xarray=xarray,
+        new = PointsNesProvidentia(comm=comm, path=path, info=info, dataset=dataset,
                                    parallel_method=parallel_method, avoid_first_hours=avoid_first_hours,
                                    avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level,
                                    create_nes=create_nes, balanced=balanced, times=times,
@@ -319,7 +316,7 @@ class PointsNesProvidentia(PointsNes):
             Portion of the variable data corresponding to the rank.
         """
 
-        nc_var = self.netcdf.variables[var_name]
+        nc_var = self.dataset.variables[var_name]
         var_dims = nc_var.dimensions
 
         # Read data in 1, 2 or 3 dimensions
diff --git a/nes/nc_projections/raster_nes.py b/nes/nc_projections/raster_nes.py
new file mode 100644
index 0000000000000000000000000000000000000000..51d5e7d317811d8d71870a71b5569f6eff6a3778
--- /dev/null
+++ b/nes/nc_projections/raster_nes.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+
+import numpy as np
+from pyproj import Proj
+from .default_nes import Nes
+
+
+class RasterNes(Nes):
+
+    def __init__(self, comm=None, path=None, info=False, **kwargs):
+        super().__init__(comm=comm, path=path, info=info, dataset=None, parallel_method='Y',
+                         balanced=False, avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None,
+                         create_nes=False, times=None, **kwargs)
diff --git a/nes/nc_projections/rotated_nes.py b/nes/nc_projections/rotated_nes.py
index eccbe03320eb11485b293d6eb167e830dc14e288..b6f86c1501d67910bc647a27a51a8c34b0d8b425 100644
--- a/nes/nc_projections/rotated_nes.py
+++ b/nes/nc_projections/rotated_nes.py
@@ -36,7 +36,7 @@ class RotatedNes(Nes):
         Tuple with the name of the dimensions of the Longitude values.
         ('rlat', 'rlon') for a rotated projection.
     """
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='Y',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
                  balanced=False, times=None, **kwargs):
         """
@@ -52,8 +52,6 @@ class RotatedNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -76,7 +74,7 @@ class RotatedNes(Nes):
 
         super(RotatedNes, self).__init__(comm=comm, path=path,
                                          info=info, dataset=dataset, balanced=balanced,
-                                         xarray=xarray, parallel_method=parallel_method,
+                                         parallel_method=parallel_method,
                                          avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                          first_level=first_level, last_level=last_level, create_nes=create_nes,
                                          times=times, **kwargs)
@@ -102,7 +100,7 @@ class RotatedNes(Nes):
         self._lon_dim = ('rlat', 'rlon')
 
     @staticmethod
-    def new(comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y', 
+    def new(comm=None, path=None, info=False, dataset=None, parallel_method='Y',
             avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None,
             create_nes=False, balanced=False, times=None, **kwargs):
         """
@@ -118,8 +116,6 @@ class RotatedNes(Nes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -136,7 +132,7 @@ class RotatedNes(Nes):
             List of times to substitute the current ones while creation.
         """
 
-        new = RotatedNes(comm=comm, path=path, info=info, dataset=dataset, xarray=xarray,
+        new = RotatedNes(comm=comm, path=path, info=info, dataset=dataset,
                          parallel_method=parallel_method, avoid_first_hours=avoid_first_hours,
                          avoid_last_hours=avoid_last_hours, first_level=first_level, last_level=last_level,
                          create_nes=create_nes, balanced=balanced, times=times, **kwargs)
diff --git a/nes/nc_projections/rotated_nested_nes.py b/nes/nc_projections/rotated_nested_nes.py
index e56f42726718c7ffce7a9b794154c542a14697b0..cafc607677293000c5ef4e7e1cd2d602318cabdd 100644
--- a/nes/nc_projections/rotated_nested_nes.py
+++ b/nes/nc_projections/rotated_nested_nes.py
@@ -4,9 +4,10 @@ import numpy as np
 from netCDF4 import Dataset
 from .rotated_nes import RotatedNes
 
+
 class RotatedNestedNes(RotatedNes):
 
-    def __init__(self, comm=None, path=None, info=False, dataset=None, xarray=False, parallel_method='Y',
+    def __init__(self, comm=None, path=None, info=False, dataset=None, parallel_method='Y',
                  avoid_first_hours=0, avoid_last_hours=0, first_level=0, last_level=None, create_nes=False,
                  balanced=False, times=None, **kwargs):
         """
@@ -22,8 +23,6 @@ class RotatedNestedNes(RotatedNes):
             Indicates if you want to get reading/writing info.
         dataset: Dataset
             NetCDF4-python Dataset to initialize the class.
-        xarray: bool:
-            (Not working) Indicates if you want to use xarray as default.
         parallel_method : str
             Indicates the parallelization method that you want. Default: 'Y'.
             Accepted values: ['X', 'Y', 'T'].
@@ -46,7 +45,7 @@ class RotatedNestedNes(RotatedNes):
 
         super(RotatedNestedNes, self).__init__(comm=comm, path=path,
                                                info=info, dataset=dataset, balanced=balanced,
-                                               xarray=xarray, parallel_method=parallel_method,
+                                               parallel_method=parallel_method,
                                                avoid_first_hours=avoid_first_hours, avoid_last_hours=avoid_last_hours,
                                                first_level=first_level, last_level=last_level, create_nes=create_nes,
                                                times=times, **kwargs)
diff --git a/nes/nes_formats/cmaq_format.py b/nes/nes_formats/cmaq_format.py
index f715b765562da6a544fc1075138f01a1e8964785..e8fd3990ee7d31dddb81a2c97025dbd76f362d19 100644
--- a/nes/nes_formats/cmaq_format.py
+++ b/nes/nes_formats/cmaq_format.py
@@ -60,7 +60,7 @@ def to_netcdf_cmaq(self, path, chunking=False, keep_open=False):
 
     # Close NetCDF
     if keep_open:
-        self.netcdf = netcdf
+        self.dataset = netcdf
     else:
         netcdf.close()
 
diff --git a/nes/nes_formats/monarch_format.py b/nes/nes_formats/monarch_format.py
index a8f5e24753395200911efa47f3b3a9e2c139227f..34a22d10af47ac09313623379d42fb1bbc5778a1 100644
--- a/nes/nes_formats/monarch_format.py
+++ b/nes/nes_formats/monarch_format.py
@@ -73,7 +73,7 @@ def to_netcdf_monarch(self, path, chunking=False, keep_open=False):
     netcdf.setncattr('Conventions', 'CF-1.7')
 
     if keep_open:
-        self.netcdf = netcdf
+        self.dataset = netcdf
     else:
         netcdf.close()
 
diff --git a/nes/nes_formats/wrf_chem_format.py b/nes/nes_formats/wrf_chem_format.py
index f37a74d4618be1c83889b5dc1c3e0bd7fad19a77..d2a71cabe0d4f919de8cd5dde76c5f9af82e9418 100644
--- a/nes/nes_formats/wrf_chem_format.py
+++ b/nes/nes_formats/wrf_chem_format.py
@@ -68,7 +68,7 @@ def to_netcdf_wrf_chem(self, path, chunking=False, keep_open=False):
 
     # Close NetCDF
     if keep_open:
-        self.netcdf = netcdf
+        self.dataset = netcdf
     else:
         netcdf.close()
 
diff --git a/requirements.txt b/requirements.txt
index 1d674e35b979cf32c8bb57d8ffac3a9b048dcdcd..7bdad0d8fa2b2e14b2c548c652ba5ef1d817d62f 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -12,8 +12,8 @@ scipy>=1.7.3
 filelock>=3.9.0
 eccodes-python~=0.9.5
 cfunits>=3.3.5
-xarray>=0.20.2
 mpi4py>=3.1.4
+# rasterio>=1.1.3
 sphinx>=7.2.6
 sphinx-rtd-theme==2.0.0
 psutil>=5.9.6
\ No newline at end of file