diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 56cb0b06..98218f42 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -45,7 +45,7 @@ jobs:
     - name: Create test env
       shell: bash -l {0}
       run: |
-        conda create -n test_env python=3.12 libgdal=3.9 -c conda-forge -c defaults -y
+        conda create -n test_env python=3.12 libgdal=3.13 -c conda-forge -c defaults -y
         conda activate test_env
         PIP_NO_BINARY=rasterio pip install .
     - name: test
diff --git a/setup.py b/setup.py
index db108502..2e0a2fd2 100644
--- a/setup.py
+++ b/setup.py
@@ -35,7 +35,7 @@
                   'sphinx>=7.1.2',
                   'scipy',
                   'pyscaffold',
-                  'gdal>=3.5.0,<3.10',
+                  'gdal>=3.5.0,<3.14',
                   'tqdm>=4.66',
                   'numpy<2'
               ],
diff --git a/src/eolab/georastertools/processing/stats.py b/src/eolab/georastertools/processing/stats.py
index f99eb52e..5b53da99 100644
--- a/src/eolab/georastertools/processing/stats.py
+++ b/src/eolab/georastertools/processing/stats.py
@@ -374,7 +374,9 @@ def plot_stats(chartfile: str, stats_per_date: Dict[datetime.datetime, gpd.GeoDa
 
             for zone in zones:
                 y = np.array(all_stats.loc[all_stats[index_name] == zone][stat_name])
-                line, = plt.plot_date(x, y, '-')
+                line, = plt.plot(x, y, '-')
+                ax = plt.gca()
+                ax.xaxis.axis_date()
                 lines.append(line)
 
             plt.title(stat_name)
@@ -383,7 +385,7 @@ def plot_stats(chartfile: str, stats_per_date: Dict[datetime.datetime, gpd.GeoDa
     plt.xlabel('date')
     plt.ylabel('values')
 
-    plt.figlegend(lines, zones, loc='lower center', ncol=2, fancybox=True, shadow=True)
+    plt.figlegend(lines[:len(zones)], zones, loc='lower center', ncol=2, fancybox=True, shadow=True)
     plt.savefig(chartfile, bbox_inches='tight')
     if display:
         plt.show()
diff --git a/src/eolab/georastertools/product/rasterproduct.py b/src/eolab/georastertools/product/rasterproduct.py
index d767dde9..d8baf8b3 100644
--- a/src/eolab/georastertools/product/rasterproduct.py
+++ b/src/eolab/georastertools/product/rasterproduct.py
@@ -460,32 +460,45 @@ def __apply_masks(self, input_vrt: Path, nb_bands: int, nb_masks: int, uuid: str
         # convert parameters defined as str to Path
         outdir = utils.to_path(self._vrt_outputdir, "/vsimem/")
         basename = utils.get_basename(self.file)
-
-        # create a tempdir for generated temporary files
-        tempdir = Path(tempfile.gettempdir())
-
-        # create a new vrt with only bands_files
-        # it must be written to disk so that we can read it and add a maskband
-        temp_image = tempdir.joinpath(f"{uuid}{basename}-temp.vrt")
-        ds = gdal.BuildVRT(temp_image.as_posix(), input_vrt.as_posix(),
-                           bandList=range(1, nb_bands + 1))
-        # free resource from GDAL
-        del ds
-
-        # Add mask band
-        _logger.debug("Adding band masks")
-        masks_index = list(range(nb_bands + 1, nb_bands + nb_masks + 1))
-        vrt_new_content = add_masks_to_vrt(temp_image, input_vrt, masks_index,
-                                           self.rastertype.maskfunc)
-        driver = gdal.GetDriverByName('VRT')
-        vrt = gdal.Open(vrt_new_content)
         masked_image = outdir.joinpath(f"{uuid}{basename}-mask.vrt")
-        copy_ds = driver.CreateCopy(masked_image.as_posix(), vrt, strict=0)
-        del copy_ds
-        # delete temp image
-        temp_image.unlink()
-        # free resource from GDAL
-        del vrt
+
+        if self._vrt_outputdir:  # output written in disk
+            ds = gdal.BuildVRT(masked_image.as_posix(), input_vrt.as_posix(),
+                               bandList=range(1, nb_bands + 1))
+            # free resource from GDAL
+            del ds
+            # Add mask band
+            _logger.debug("Adding band masks")
+            masks_index = list(range(nb_bands + 1, nb_bands + nb_masks + 1))
+            add_masks_to_vrt(masked_image, input_vrt, masks_index,
+                             self.rastertype.maskfunc)
+
+        else:  # output written in memory
+            # create a tempdir for generated temporary files
+            tempdir = Path(tempfile.gettempdir())
+
+            # create a new vrt with only bands_files
+            # it must be written to disk so that we can read it and add a maskband
+            temp_image = tempdir.joinpath(f"{uuid}{basename}-temp.vrt")
+            ds = gdal.BuildVRT(temp_image.as_posix(), input_vrt.as_posix(),
+                               bandList=range(1, nb_bands + 1))
+            # free resource from GDAL
+            del ds
+
+            # Add mask band
+            _logger.debug("Adding band masks")
+            masks_index = list(range(nb_bands + 1, nb_bands + nb_masks + 1))
+            add_masks_to_vrt(temp_image, input_vrt, masks_index,
+                             self.rastertype.maskfunc)
+            driver = gdal.GetDriverByName('VRT')
+            vrt = gdal.Open(str(temp_image))
+            masked_image = outdir.joinpath(f"{uuid}{basename}-mask.vrt")
+            copy_ds = driver.CreateCopy(masked_image.as_posix(), vrt, strict=0)
+            del copy_ds
+            # delete temp image
+            temp_image.unlink()
+            # free resource from GDAL
+            del vrt
 
         return masked_image
 
diff --git a/src/eolab/georastertools/product/vrt.py b/src/eolab/georastertools/product/vrt.py
index c24b2ad7..622bf302 100644
--- a/src/eolab/georastertools/product/vrt.py
+++ b/src/eolab/georastertools/product/vrt.py
@@ -103,7 +103,7 @@ def {funcname}(in_ar, out_ar, xoff, yoff, xsize, ysize,\
 
 
 def add_masks_to_vrt(src_vrt: Union[Path, str], maskfile: Union[Path, str], bands: List[int] = [1],
-                     funcname: str = None, funcdef: str = None) -> str:
+                     funcname: str = None, funcdef: str = None):
     """Adds a mask bands to the vrt.
 
     Args:
@@ -119,9 +119,6 @@ def add_masks_to_vrt(src_vrt: Union[Path, str], maskfile: Union[Path, str], band
         funcdef (str, optional, default=None):
             Function definition (without its signature). If None, the function funcname
             must be available in the scope
-
-    Returns:
-        (str): XML content (vrt format) with the added mask band.
     """
     svrt = src_vrt.as_posix() if isinstance(src_vrt, Path) else src_vrt
     mask = maskfile.as_posix() if isinstance(maskfile, Path) else maskfile
@@ -137,7 +134,7 @@ def add_masks_to_vrt(src_vrt: Union[Path, str], maskfile: Union[Path, str], band
             _func2vrt(vrtrasterband, funcname, funcdef)
         _file2vrt(vrtrasterband, mask, bands)
 
-        return ET.tostring(root)
+    tree.write(svrt)
 
 
 def set_band_descriptions(src_vrt: Union[Path, str], descriptions: List[str]):
diff --git a/tests/cmptools.py b/tests/cmptools.py
index 5ee763f9..d6003d58 100644
--- a/tests/cmptools.py
+++ b/tests/cmptools.py
@@ -70,6 +70,11 @@ def equals(self, el, tolerance):
             p2 = shapely.wkt.loads(t2)
             if p1.hausdorff_distance(p2) > tolerance:
                 return False
+        elif "GeoTransform" in self.name:
+            transform1 = [float(x) for x in t1.split(",")]
+            transform2 = [float(x) for x in t2.split(",")]
+            if transform1 != transform2:
+                return False
         else:
             if t2 != t1:
                 return False
diff --git a/tests/utils4test.py b/tests/utils4test.py
index 88cdd998..5d56d951 100644
--- a/tests/utils4test.py
+++ b/tests/utils4test.py
@@ -64,7 +64,7 @@ def basename(infile):
     return file.name if suffix == 0 else file.name[:-suffix]
 
 
-def cmpfiles(a : str, b : str, common : list, tolerance : float =1e-9, **kwargs) -> tuple:
+def cmpfiles(a : str, b : str, common : list, tolerance : float=1e-9, **kwargs) -> tuple:
     """
     Compare common files in two directories.