diff --git a/autotest/pymod/gdaltest.py b/autotest/pymod/gdaltest.py
index cd0389320fba..db37bd2eeab8 100755
--- a/autotest/pymod/gdaltest.py
+++ b/autotest/pymod/gdaltest.py
@@ -1912,7 +1912,10 @@ def runexternal(
     encoding="latin1",
 ):
     # pylint: disable=unused-argument
-    command = shlex.split(cmd)
+    if sys.platform == "win32":
+        command = cmd
+    else:
+        command = shlex.split(cmd)
     if strin is None:
         p = subprocess.Popen(command, stdout=subprocess.PIPE)
     else:
@@ -1949,7 +1952,10 @@ def _read_in_thread(f, q):
 
 def runexternal_out_and_err(cmd, check_memleak=True, encoding="ascii"):
     # pylint: disable=unused-argument
-    command = shlex.split(cmd)
+    if sys.platform == "win32":
+        command = cmd
+    else:
+        command = shlex.split(cmd)
     p = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
 
     if p.stdout is not None:
diff --git a/autotest/pyscripts/test_gdal_edit.py b/autotest/pyscripts/test_gdal_edit.py
index 035f118a68fb..c0b716ea0176 100755
--- a/autotest/pyscripts/test_gdal_edit.py
+++ b/autotest/pyscripts/test_gdal_edit.py
@@ -173,7 +173,9 @@ def test_gdal_edit_py_3(script_path):
     shutil.copy(test_py_scripts.get_data_path("gcore") + "byte.tif", filename)
 
     try:
-        test_py_scripts.run_py_script(script_path, "gdal_edit", filename + " -a_srs ''")
+        test_py_scripts.run_py_script(
+            script_path, "gdal_edit", filename + " -a_srs None"
+        )
 
         ds = gdal.Open(filename)
         wkt = ds.GetProjectionRef()
diff --git a/autotest/pyscripts/test_gdal_merge.py b/autotest/pyscripts/test_gdal_merge.py
index 94408f9178f5..96b07f70d346 100755
--- a/autotest/pyscripts/test_gdal_merge.py
+++ b/autotest/pyscripts/test_gdal_merge.py
@@ -47,21 +47,57 @@ def script_path():
     return test_py_scripts.get_py_script("gdal_merge")
 
 
+@pytest.fixture(scope="module")
+def sample_tifs(tmp_path_factory):
+    tmpdir = tmp_path_factory.mktemp("tmp")
+
+    drv = gdal.GetDriverByName("GTiff")
+    srs = osr.SpatialReference()
+    srs.SetWellKnownGeogCS("WGS84")
+    wkt = srs.ExportToWkt()
+
+    sample1_tif = str(tmpdir / "in1.tif")
+    with drv.Create(sample1_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
+        ds.GetRasterBand(1).Fill(0)
+
+    sample2_tif = str(tmpdir / "in2.tif")
+    with drv.Create(sample2_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([3, 0.1, 0, 49, 0, -0.1])
+        ds.GetRasterBand(1).Fill(63)
+
+    sample3_tif = str(tmpdir / "in3.tif")
+    with drv.Create(sample3_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([2, 0.1, 0, 48, 0, -0.1])
+        ds.GetRasterBand(1).Fill(127)
+
+    sample4_tif = str(tmpdir / "in4.tif")
+    with drv.Create(sample4_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([3, 0.1, 0, 48, 0, -0.1])
+        ds.GetRasterBand(1).Fill(255)
+
+    yield (sample1_tif, sample2_tif, sample3_tif, sample4_tif)
+
+
 ###############################################################################
 # Basic test
 
 
-def test_gdal_merge_1(script_path):
+def test_gdal_merge_1(script_path, tmp_path):
+
+    output_tif = str(tmp_path / "test_gdal_merge_1.tif")
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_merge",
-        "-o tmp/test_gdal_merge_1.tif "
-        + test_py_scripts.get_data_path("gcore")
-        + "byte.tif",
+        f"-o {output_tif} " + test_py_scripts.get_data_path("gcore") + "byte.tif",
     )
 
-    ds = gdal.Open("tmp/test_gdal_merge_1.tif")
+    ds = gdal.Open(output_tif)
     assert ds.GetRasterBand(1).Checksum() == 4672
     ds = None
 
@@ -70,44 +106,17 @@ def test_gdal_merge_1(script_path):
 # Merge 4 tiles
 
 
-def test_gdal_merge_2(script_path):
-
-    drv = gdal.GetDriverByName("GTiff")
-    srs = osr.SpatialReference()
-    srs.SetWellKnownGeogCS("WGS84")
-    wkt = srs.ExportToWkt()
-
-    ds = drv.Create("tmp/in1.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
-    ds.GetRasterBand(1).Fill(0)
-    ds = None
-
-    ds = drv.Create("tmp/in2.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([3, 0.1, 0, 49, 0, -0.1])
-    ds.GetRasterBand(1).Fill(63)
-    ds = None
+def test_gdal_merge_2(script_path, tmp_path, sample_tifs):
 
-    ds = drv.Create("tmp/in3.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([2, 0.1, 0, 48, 0, -0.1])
-    ds.GetRasterBand(1).Fill(127)
-    ds = None
-
-    ds = drv.Create("tmp/in4.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([3, 0.1, 0, 48, 0, -0.1])
-    ds.GetRasterBand(1).Fill(255)
-    ds = None
+    output_tif = str(tmp_path / "test_gdal_merge_2.tif")
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_merge",
-        "-q -o tmp/test_gdal_merge_2.tif tmp/in1.tif tmp/in2.tif tmp/in3.tif tmp/in4.tif",
+        f"-q -o {output_tif} {' '.join(sample_tifs)}",
     )
 
-    ds = gdal.Open("tmp/test_gdal_merge_2.tif")
+    ds = gdal.Open(output_tif)
     assert ds.GetProjectionRef().find("WGS 84") != -1, "Expected WGS 84\nGot : %s" % (
         ds.GetProjectionRef()
     )
@@ -133,15 +142,17 @@ def test_gdal_merge_2(script_path):
 # Test -separate and -v options
 
 
-def test_gdal_merge_3(script_path):
+def test_gdal_merge_3(script_path, tmp_path, sample_tifs):
+
+    output_tif = str(tmp_path / "test_gdal_merge_3.tif")
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_merge",
-        "-separate -v -o tmp/test_gdal_merge_3.tif tmp/in1.tif tmp/in2.tif tmp/in3.tif tmp/in4.tif",
+        f"-separate -v -o {output_tif} {' '.join(sample_tifs)}",
     )
 
-    ds = gdal.Open("tmp/test_gdal_merge_3.tif")
+    ds = gdal.Open(output_tif)
     assert ds.GetProjectionRef().find("WGS 84") != -1, "Expected WGS 84\nGot : %s" % (
         ds.GetProjectionRef()
     )
@@ -167,15 +178,17 @@ def test_gdal_merge_3(script_path):
 # Test -init option
 
 
-def test_gdal_merge_4(script_path):
+def test_gdal_merge_4(script_path, tmp_path, sample_tifs):
+
+    output_tif = str(tmp_path / "test_gdal_merge_4.tif")
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_merge",
-        "-init 255 -o tmp/test_gdal_merge_4.tif tmp/in2.tif tmp/in3.tif",
+        f"-init 255 -o {output_tif} {sample_tifs[1]} {sample_tifs[2]}",
     )
 
-    ds = gdal.Open("tmp/test_gdal_merge_4.tif")
+    ds = gdal.Open(output_tif)
 
     assert ds.GetRasterBand(1).Checksum() == 4725, "Wrong checksum"
 
@@ -184,7 +197,7 @@ def test_gdal_merge_4(script_path):
 # Test merging with alpha band (#3669)
 
 
-def test_gdal_merge_5(script_path):
+def test_gdal_merge_5(script_path, tmp_path):
     gdal_array = pytest.importorskip("osgeo.gdal_array")
     try:
         gdal_array.BandRasterIONumPy
@@ -196,71 +209,46 @@ def test_gdal_merge_5(script_path):
     srs.SetWellKnownGeogCS("WGS84")
     wkt = srs.ExportToWkt()
 
-    ds = drv.Create("tmp/in5.tif", 10, 10, 4)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
-    ds.GetRasterBand(1).Fill(255)
-    ds = None
+    input1_tif = str(tmp_path / "in5.tif")
+    with drv.Create(input1_tif, 10, 10, 4) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
+        ds.GetRasterBand(1).Fill(255)
 
-    ds = drv.Create("tmp/in6.tif", 10, 10, 4)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
-    ds.GetRasterBand(2).Fill(255)
-    ds.GetRasterBand(4).Fill(255)
-    cs = ds.GetRasterBand(4).Checksum()
-    ds = None
+    input2_tif = str(tmp_path / "in6.tif")
+    with drv.Create(input2_tif, 10, 10, 4) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
+        ds.GetRasterBand(2).Fill(255)
+        ds.GetRasterBand(4).Fill(255)
+        cs = ds.GetRasterBand(4).Checksum()
+
+    output_tif = str(tmp_path / "test_gdal_merge_5.tif")
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_merge",
-        " -o tmp/test_gdal_merge_5.tif tmp/in5.tif tmp/in6.tif",
+        f" -o {output_tif} {input1_tif} {input2_tif}",
     )
 
-    ds = gdal.Open("tmp/test_gdal_merge_5.tif")
+    with gdal.Open(output_tif) as ds:
 
-    assert ds.GetRasterBand(1).Checksum() == 0, "Wrong checksum"
-    assert ds.GetRasterBand(2).Checksum() == cs, "Wrong checksum"
-    assert ds.GetRasterBand(3).Checksum() == 0, "Wrong checksum"
-    assert ds.GetRasterBand(4).Checksum() == cs, "Wrong checksum"
-    ds = None
+        assert ds.GetRasterBand(1).Checksum() == 0, "Wrong checksum"
+        assert ds.GetRasterBand(2).Checksum() == cs, "Wrong checksum"
+        assert ds.GetRasterBand(3).Checksum() == 0, "Wrong checksum"
+        assert ds.GetRasterBand(4).Checksum() == cs, "Wrong checksum"
 
-    os.unlink("tmp/test_gdal_merge_5.tif")
+    os.unlink(output_tif)
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_merge",
-        " -o tmp/test_gdal_merge_5.tif tmp/in6.tif tmp/in5.tif",
+        f" -o {output_tif} {input1_tif} {input2_tif}",
     )
 
-    ds = gdal.Open("tmp/test_gdal_merge_5.tif")
+    with gdal.Open(output_tif) as ds:
 
-    assert ds.GetRasterBand(1).Checksum() == 0, "Wrong checksum"
-    assert ds.GetRasterBand(2).Checksum() == cs, "Wrong checksum"
-    assert ds.GetRasterBand(3).Checksum() == 0, "Wrong checksum"
-    assert ds.GetRasterBand(4).Checksum() == cs, "Wrong checksum"
-
-
-###############################################################################
-# Cleanup
-
-
-def test_gdal_merge_cleanup():
-
-    lst = [
-        "tmp/test_gdal_merge_1.tif",
-        "tmp/test_gdal_merge_2.tif",
-        "tmp/test_gdal_merge_3.tif",
-        "tmp/test_gdal_merge_4.tif",
-        "tmp/test_gdal_merge_5.tif",
-        "tmp/in1.tif",
-        "tmp/in2.tif",
-        "tmp/in3.tif",
-        "tmp/in4.tif",
-        "tmp/in5.tif",
-        "tmp/in6.tif",
-    ]
-    for filename in lst:
-        try:
-            os.remove(filename)
-        except OSError:
-            pass
+        assert ds.GetRasterBand(1).Checksum() == 0, "Wrong checksum"
+        assert ds.GetRasterBand(2).Checksum() == cs, "Wrong checksum"
+        assert ds.GetRasterBand(3).Checksum() == 0, "Wrong checksum"
+        assert ds.GetRasterBand(4).Checksum() == cs, "Wrong checksum"
diff --git a/autotest/pyscripts/test_gdal_pansharpen.py b/autotest/pyscripts/test_gdal_pansharpen.py
index 15fee660c75a..14596a0c60e5 100755
--- a/autotest/pyscripts/test_gdal_pansharpen.py
+++ b/autotest/pyscripts/test_gdal_pansharpen.py
@@ -46,38 +46,48 @@ def script_path():
     return test_py_scripts.get_py_script("gdal_pansharpen")
 
 
+@pytest.fixture(scope="module")
+def small_world_pan_tif(tmp_path_factory):
+
+    small_world_pan_tif = str(tmp_path_factory.mktemp("tmp") / "small_world_pan.tif")
+
+    with gdal.Open(
+        test_py_scripts.get_data_path("gdrivers") + "small_world.tif"
+    ) as src_ds:
+        src_data = src_ds.GetRasterBand(1).ReadRaster()
+        gt = src_ds.GetGeoTransform()
+        wkt = src_ds.GetProjectionRef()
+
+    with gdal.GetDriverByName("GTiff").Create(small_world_pan_tif, 800, 400) as pan_ds:
+        gt = [gt[i] for i in range(len(gt))]
+        gt[1] *= 0.5
+        gt[5] *= 0.5
+        pan_ds.SetGeoTransform(gt)
+        pan_ds.SetProjection(wkt)
+        pan_ds.GetRasterBand(1).WriteRaster(0, 0, 800, 400, src_data, 400, 200)
+
+    return small_world_pan_tif
+
+
 ###############################################################################
 # Simple test
 
 
-def test_gdal_pansharpen_1(script_path):
+def test_gdal_pansharpen_1(script_path, tmp_path, small_world_pan_tif):
 
-    src_ds = gdal.Open(test_py_scripts.get_data_path("gdrivers") + "small_world.tif")
-    src_data = src_ds.GetRasterBand(1).ReadRaster()
-    gt = src_ds.GetGeoTransform()
-    wkt = src_ds.GetProjectionRef()
-    src_ds = None
-    pan_ds = gdal.GetDriverByName("GTiff").Create("tmp/small_world_pan.tif", 800, 400)
-    gt = [gt[i] for i in range(len(gt))]
-    gt[1] *= 0.5
-    gt[5] *= 0.5
-    pan_ds.SetGeoTransform(gt)
-    pan_ds.SetProjection(wkt)
-    pan_ds.GetRasterBand(1).WriteRaster(0, 0, 800, 400, src_data, 400, 200)
-    pan_ds = None
+    out_tif = str(tmp_path / "out.tif")
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_pansharpen",
-        " tmp/small_world_pan.tif "
+        f" {small_world_pan_tif} "
         + test_py_scripts.get_data_path("gdrivers")
-        + "small_world.tif tmp/out.tif",
+        + "small_world.tif "
+        + out_tif,
     )
 
-    ds = gdal.Open("tmp/out.tif")
-    cs = [ds.GetRasterBand(i + 1).Checksum() for i in range(ds.RasterCount)]
-    ds = None
-    gdal.GetDriverByName("GTiff").Delete("tmp/out.tif")
+    with gdal.Open(out_tif) as ds:
+        cs = [ds.GetRasterBand(i + 1).Checksum() for i in range(ds.RasterCount)]
 
     assert cs in ([4735, 10000, 9742], [4731, 9991, 9734])  # s390x or graviton2
 
@@ -86,32 +96,26 @@ def test_gdal_pansharpen_1(script_path):
 # Full options
 
 
-def test_gdal_pansharpen_2(script_path):
+def test_gdal_pansharpen_2(script_path, tmp_path, small_world_pan_tif):
+
+    out_vrt = str(tmp_path / "out.vrt")
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_pansharpen",
-        " -q -b 3 -b 1 -bitdepth 8 -threads ALL_CPUS -spat_adjust union -w 0.33333333333333333 -w 0.33333333333333333 -w 0.33333333333333333 -of VRT -r cubic tmp/small_world_pan.tif "
+        " -q -b 3 -b 1 -bitdepth 8 -threads ALL_CPUS -spat_adjust union -w 0.33333333333333333 -w 0.33333333333333333 -w 0.33333333333333333 -of VRT -r cubic "
+        + small_world_pan_tif
+        + " "
         + test_py_scripts.get_data_path("gdrivers")
         + "small_world.tif,band=1 "
         + test_py_scripts.get_data_path("gdrivers")
         + "small_world.tif,band=2 "
         + test_py_scripts.get_data_path("gdrivers")
-        + "small_world.tif,band=3 tmp/out.vrt",
+        + "small_world.tif,band=3 "
+        + out_vrt,
     )
 
-    ds = gdal.Open("tmp/out.vrt")
-    cs = [ds.GetRasterBand(i + 1).Checksum() for i in range(ds.RasterCount)]
-    ds = None
-    gdal.GetDriverByName("VRT").Delete("tmp/out.vrt")
+    with gdal.Open(out_vrt) as ds:
+        cs = [ds.GetRasterBand(i + 1).Checksum() for i in range(ds.RasterCount)]
 
     assert cs in ([9742, 4735], [9734, 4731])  # s390x or graviton2
-
-
-###############################################################################
-# Cleanup
-
-
-def test_gdal_pansharpen_cleanup():
-
-    gdal.GetDriverByName("GTiff").Delete("tmp/small_world_pan.tif")
diff --git a/autotest/pyscripts/test_gdal_retile.py b/autotest/pyscripts/test_gdal_retile.py
index 58a917280ef7..55508fb22073 100755
--- a/autotest/pyscripts/test_gdal_retile.py
+++ b/autotest/pyscripts/test_gdal_retile.py
@@ -29,7 +29,6 @@
 ###############################################################################
 
 import os
-import shutil
 
 import pytest
 import test_py_scripts
@@ -51,60 +50,53 @@ def script_path():
 # Test gdal_retile.py
 
 
-def test_gdal_retile_1(script_path):
+def test_gdal_retile_1(script_path, tmp_path):
 
-    try:
-        os.mkdir("tmp/outretile")
-    except OSError:
-        pass
+    out_dir = tmp_path / "outretile2"
+    out_dir.mkdir()
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_retile",
-        "-v -levels 2 -r bilinear -targetDir tmp/outretile "
+        f"-v -levels 2 -r bilinear -targetDir {out_dir} "
         + test_py_scripts.get_data_path("gcore")
         + "byte.tif",
     )
 
-    ds = gdal.Open("tmp/outretile/byte_1_1.tif")
-    assert ds.GetRasterBand(1).Checksum() == 4672
-    ds = None
+    with gdal.Open(f"{out_dir}/byte_1_1.tif") as ds:
+        assert ds.GetRasterBand(1).Checksum() == 4672
 
-    ds = gdal.Open("tmp/outretile/1/byte_1_1.tif")
-    assert ds.RasterXSize == 10
-    # if ds.GetRasterBand(1).Checksum() != 1152:
-    #    print(ds.GetRasterBand(1).Checksum())
-    #    return 'fail'
-    ds = None
+    with gdal.Open(f"{out_dir}/1/byte_1_1.tif") as ds:
+        assert ds.RasterXSize == 10
+        # if ds.GetRasterBand(1).Checksum() != 1152:
+        #    print(ds.GetRasterBand(1).Checksum())
+        #    return 'fail'
 
-    ds = gdal.Open("tmp/outretile/2/byte_1_1.tif")
-    assert ds.RasterXSize == 5
-    # if ds.GetRasterBand(1).Checksum() != 215:
-    #    print(ds.GetRasterBand(1).Checksum())
-    #    return 'fail'
-    ds = None
+    with gdal.Open(f"{out_dir}/2/byte_1_1.tif") as ds:
+        assert ds.RasterXSize == 5
+        # if ds.GetRasterBand(1).Checksum() != 215:
+        #    print(ds.GetRasterBand(1).Checksum())
+        #    return 'fail'
 
 
 ###############################################################################
 # Test gdal_retile.py with RGBA dataset
 
 
-def test_gdal_retile_2(script_path):
+def test_gdal_retile_2(script_path, tmp_path):
 
-    try:
-        os.mkdir("tmp/outretile2")
-    except OSError:
-        pass
+    out_dir = tmp_path / "outretile2"
+    out_dir.mkdir()
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_retile",
-        "-v -levels 2 -r bilinear -targetDir tmp/outretile2 "
+        f"-v -levels 2 -r bilinear -targetDir {out_dir} "
         + test_py_scripts.get_data_path("gcore")
         + "rgba.tif",
     )
 
-    ds = gdal.Open("tmp/outretile2/2/rgba_1_1.tif")
+    ds = gdal.Open(f"{out_dir}/2/rgba_1_1.tif")
     assert ds.GetRasterBand(1).Checksum() == 35, "wrong checksum for band 1"
     assert ds.GetRasterBand(4).Checksum() == 35, "wrong checksum for band 4"
     ds = None
@@ -114,7 +106,7 @@ def test_gdal_retile_2(script_path):
 # Test gdal_retile.py with input images of different pixel sizes
 
 
-def test_gdal_retile_3(script_path):
+def test_gdal_retile_3(script_path, tmp_path):
 
     drv = gdal.GetDriverByName("GTiff")
     srs = osr.SpatialReference()
@@ -140,35 +132,33 @@ def test_gdal_retile_3(script_path):
     #  0 N ---------------
     #      0 E           30 E
 
-    ds = drv.Create("tmp/in1.tif", 100, 100, 1)
-    px1_x = 30.0 / ds.RasterXSize
-    px1_y = 30.0 / ds.RasterYSize
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([0, px1_x, 0, 30, 0, -px1_y])
-    ds.GetRasterBand(1).Fill(0)
-    ds = None
-
-    ds = drv.Create("tmp/in2.tif", 50, 50, 1)
-    px2_x = 30.0 / ds.RasterXSize
-    px2_y = 30.0 / ds.RasterYSize
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([0, px2_x, 0, 60, 0, -px2_y])
-    ds.GetRasterBand(1).Fill(42)
-    ds = None
-
-    try:
-        os.mkdir("tmp/outretile3")
-    except OSError:
-        pass
+    in1_tif = str(tmp_path / "in1.tif")
+    with drv.Create(in1_tif, 100, 100, 1) as ds:
+        px1_x = 30.0 / ds.RasterXSize
+        px1_y = 30.0 / ds.RasterYSize
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([0, px1_x, 0, 30, 0, -px1_y])
+        ds.GetRasterBand(1).Fill(0)
+
+    in2_tif = str(tmp_path / "in2.tif")
+    with drv.Create(in2_tif, 50, 50, 1) as ds:
+        px2_x = 30.0 / ds.RasterXSize
+        px2_y = 30.0 / ds.RasterYSize
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([0, px2_x, 0, 60, 0, -px2_y])
+        ds.GetRasterBand(1).Fill(42)
+
+    out_dir = tmp_path / "outretile4"
+    out_dir.mkdir()
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_retile",
-        "-v -levels 2 -r bilinear -targetDir tmp/outretile3 tmp/in1.tif tmp/in2.tif",
+        f"-v -levels 2 -r bilinear -targetDir {out_dir} {in1_tif} {in2_tif}",
     )
 
-    ds = gdal.Open("tmp/outretile3/in1_1_1.tif")
-    assert ds.GetProjectionRef().find("WGS 84") != -1, "Expected WGS 84\nGot : %s" % (
+    ds = gdal.Open(f"{out_dir}/in1_1_1.tif")
+    assert "WGS 84" in ds.GetProjectionRef(), "Expected WGS 84\nGot : %s" % (
         ds.GetProjectionRef()
     )
 
@@ -193,42 +183,40 @@ def test_gdal_retile_3(script_path):
 # Test gdal_retile.py -overlap
 
 
-def test_gdal_retile_4(script_path):
+def test_gdal_retile_4(script_path, tmp_path):
 
-    try:
-        os.mkdir("tmp/outretile4")
-    except OSError:
-        pass
+    out_dir = tmp_path / "outretile4"
+    out_dir.mkdir()
 
     test_py_scripts.run_py_script(
         script_path,
         "gdal_retile",
-        "-v -ps 8 7 -overlap 3 -targetDir tmp/outretile4 "
+        f"-v -ps 8 7 -overlap 3 -targetDir {out_dir} "
         + test_py_scripts.get_data_path("gcore")
         + "byte.tif",
     )
 
     expected_results = [
-        ["tmp/outretile4/byte_1_1.tif", 8, 7],
-        ["tmp/outretile4/byte_1_2.tif", 8, 7],
-        ["tmp/outretile4/byte_1_3.tif", 8, 7],
-        ["tmp/outretile4/byte_1_4.tif", 5, 7],
-        ["tmp/outretile4/byte_2_1.tif", 8, 7],
-        ["tmp/outretile4/byte_2_2.tif", 8, 7],
-        ["tmp/outretile4/byte_2_3.tif", 8, 7],
-        ["tmp/outretile4/byte_2_4.tif", 5, 7],
-        ["tmp/outretile4/byte_3_1.tif", 8, 7],
-        ["tmp/outretile4/byte_3_2.tif", 8, 7],
-        ["tmp/outretile4/byte_3_3.tif", 8, 7],
-        ["tmp/outretile4/byte_3_4.tif", 5, 7],
-        ["tmp/outretile4/byte_4_1.tif", 8, 7],
-        ["tmp/outretile4/byte_4_2.tif", 8, 7],
-        ["tmp/outretile4/byte_4_3.tif", 8, 7],
-        ["tmp/outretile4/byte_4_4.tif", 5, 7],
-        ["tmp/outretile4/byte_5_1.tif", 8, 4],
-        ["tmp/outretile4/byte_5_2.tif", 8, 4],
-        ["tmp/outretile4/byte_5_3.tif", 8, 4],
-        ["tmp/outretile4/byte_5_4.tif", 5, 4],
+        [f"{out_dir}/byte_1_1.tif", 8, 7],
+        [f"{out_dir}/byte_1_2.tif", 8, 7],
+        [f"{out_dir}/byte_1_3.tif", 8, 7],
+        [f"{out_dir}/byte_1_4.tif", 5, 7],
+        [f"{out_dir}/byte_2_1.tif", 8, 7],
+        [f"{out_dir}/byte_2_2.tif", 8, 7],
+        [f"{out_dir}/byte_2_3.tif", 8, 7],
+        [f"{out_dir}/byte_2_4.tif", 5, 7],
+        [f"{out_dir}/byte_3_1.tif", 8, 7],
+        [f"{out_dir}/byte_3_2.tif", 8, 7],
+        [f"{out_dir}/byte_3_3.tif", 8, 7],
+        [f"{out_dir}/byte_3_4.tif", 5, 7],
+        [f"{out_dir}/byte_4_1.tif", 8, 7],
+        [f"{out_dir}/byte_4_2.tif", 8, 7],
+        [f"{out_dir}/byte_4_3.tif", 8, 7],
+        [f"{out_dir}/byte_4_4.tif", 5, 7],
+        [f"{out_dir}/byte_5_1.tif", 8, 4],
+        [f"{out_dir}/byte_5_2.tif", 8, 4],
+        [f"{out_dir}/byte_5_3.tif", 8, 4],
+        [f"{out_dir}/byte_5_4.tif", 5, 4],
     ]
 
     for (filename, width, height) in expected_results:
@@ -240,32 +228,32 @@ def test_gdal_retile_4(script_path):
     test_py_scripts.run_py_script(
         script_path,
         "gdal_retile",
-        "-v -levels 1 -ps 8 8 -overlap 4 -targetDir tmp/outretile4 "
+        f"-v -levels 1 -ps 8 8 -overlap 4 -targetDir {out_dir} "
         + test_py_scripts.get_data_path("gcore")
         + "byte.tif",
     )
 
     expected_results = [
-        ["tmp/outretile4/byte_1_1.tif", 8, 8],
-        ["tmp/outretile4/byte_1_2.tif", 8, 8],
-        ["tmp/outretile4/byte_1_3.tif", 8, 8],
-        ["tmp/outretile4/byte_1_4.tif", 8, 8],
-        ["tmp/outretile4/byte_2_1.tif", 8, 8],
-        ["tmp/outretile4/byte_2_2.tif", 8, 8],
-        ["tmp/outretile4/byte_2_3.tif", 8, 8],
-        ["tmp/outretile4/byte_2_4.tif", 8, 8],
-        ["tmp/outretile4/byte_3_1.tif", 8, 8],
-        ["tmp/outretile4/byte_3_2.tif", 8, 8],
-        ["tmp/outretile4/byte_3_3.tif", 8, 8],
-        ["tmp/outretile4/byte_3_4.tif", 8, 8],
-        ["tmp/outretile4/byte_4_1.tif", 8, 8],
-        ["tmp/outretile4/byte_4_2.tif", 8, 8],
-        ["tmp/outretile4/byte_4_3.tif", 8, 8],
-        ["tmp/outretile4/byte_4_4.tif", 8, 8],
-        ["tmp/outretile4/1/byte_1_1.tif", 8, 8],
-        ["tmp/outretile4/1/byte_1_2.tif", 6, 8],
-        ["tmp/outretile4/1/byte_2_1.tif", 8, 6],
-        ["tmp/outretile4/1/byte_2_2.tif", 6, 6],
+        [f"{out_dir}/byte_1_1.tif", 8, 8],
+        [f"{out_dir}/byte_1_2.tif", 8, 8],
+        [f"{out_dir}/byte_1_3.tif", 8, 8],
+        [f"{out_dir}/byte_1_4.tif", 8, 8],
+        [f"{out_dir}/byte_2_1.tif", 8, 8],
+        [f"{out_dir}/byte_2_2.tif", 8, 8],
+        [f"{out_dir}/byte_2_3.tif", 8, 8],
+        [f"{out_dir}/byte_2_4.tif", 8, 8],
+        [f"{out_dir}/byte_3_1.tif", 8, 8],
+        [f"{out_dir}/byte_3_2.tif", 8, 8],
+        [f"{out_dir}/byte_3_3.tif", 8, 8],
+        [f"{out_dir}/byte_3_4.tif", 8, 8],
+        [f"{out_dir}/byte_4_1.tif", 8, 8],
+        [f"{out_dir}/byte_4_2.tif", 8, 8],
+        [f"{out_dir}/byte_4_3.tif", 8, 8],
+        [f"{out_dir}/byte_4_4.tif", 8, 8],
+        [f"{out_dir}/1/byte_1_1.tif", 8, 8],
+        [f"{out_dir}/1/byte_1_2.tif", 6, 8],
+        [f"{out_dir}/1/byte_2_1.tif", 8, 6],
+        [f"{out_dir}/1/byte_2_2.tif", 6, 6],
     ]
 
     for (filename, width, height) in expected_results:
@@ -279,7 +267,7 @@ def test_gdal_retile_4(script_path):
 # Test gdal_retile.py with input having a NoData value
 
 
-def test_gdal_retile_5(script_path):
+def test_gdal_retile_5(script_path, tmp_path):
 
     np = pytest.importorskip("numpy")
 
@@ -291,42 +279,39 @@ def test_gdal_retile_5(script_path):
     srs.SetWellKnownGeogCS("WGS84")
     wkt = srs.ExportToWkt()
 
-    ds = drv.Create("tmp/in5.tif", 2, 2, 1, gdal.GDT_Float32)
-    px1_x = 0.1 / ds.RasterXSize
-    px1_y = 0.1 / ds.RasterYSize
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([0, px1_x, 0, 30, 0, -px1_y])
-    raster_band = ds.GetRasterBand(1)
-    raster_band.SetNoDataValue(nodata_value)
-    raster_band.WriteArray(raster_array)
-    raster_band = None
-    ds = None
+    input_tif = str(tmp_path / "in5.tif")
 
-    try:
-        os.mkdir("tmp/outretile5")
-    except OSError:
-        pass
+    with drv.Create(input_tif, 2, 2, 1, gdal.GDT_Float32) as ds:
+        px1_x = 0.1 / ds.RasterXSize
+        px1_y = 0.1 / ds.RasterYSize
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([0, px1_x, 0, 30, 0, -px1_y])
+        raster_band = ds.GetRasterBand(1)
+        raster_band.SetNoDataValue(nodata_value)
+        raster_band.WriteArray(raster_array)
+        raster_band = None
+
+    out_dir = tmp_path / "outretile5"
+    out_dir.mkdir()
 
     test_py_scripts.run_py_script(
-        script_path, "gdal_retile", "-v -targetDir tmp/outretile5 tmp/in5.tif"
+        script_path, "gdal_retile", f"-v -targetDir {out_dir} {input_tif}"
     )
 
-    ds = gdal.Open("tmp/outretile5/in5_1_1.tif")
-    raster_band = ds.GetRasterBand(1)
+    with gdal.Open(str(out_dir / "in5_1_1.tif")) as ds:
+        raster_band = ds.GetRasterBand(1)
 
-    assert (
-        raster_band.GetNoDataValue() == nodata_value
-    ), "Wrong nodata value.\nExpected %f, Got: %f" % (
-        nodata_value,
-        raster_band.GetNoDataValue(),
-    )
-
-    min_val, max_val = raster_band.ComputeRasterMinMax()
-    assert max_val, "Wrong maximum value.\nExpected 2.0, Got: %f" % max_val
+        assert (
+            raster_band.GetNoDataValue() == nodata_value
+        ), "Wrong nodata value.\nExpected %f, Got: %f" % (
+            nodata_value,
+            raster_band.GetNoDataValue(),
+        )
 
-    assert min_val == -1.0, "Wrong minimum value.\nExpected -1.0, Got: %f" % min_val
+        min_val, max_val = raster_band.ComputeRasterMinMax()
+        assert max_val, "Wrong maximum value.\nExpected 2.0, Got: %f" % max_val
 
-    ds = None
+        assert min_val == -1.0, "Wrong minimum value.\nExpected -1.0, Got: %f" % min_val
 
 
 ###############################################################################
@@ -334,71 +319,24 @@ def test_gdal_retile_5(script_path):
 
 
 @pytest.mark.require_driver("PNG")
-def test_gdal_retile_png(script_path):
-
-    out_dirname = os.path.join("tmp", "outretile_png")
-    os.mkdir(out_dirname)
-
-    try:
-        test_py_scripts.run_py_script(
-            script_path,
-            "gdal_retile",
-            "-v -levels 2 -r bilinear -of PNG -targetDir "
-            + '"'
-            + out_dirname
-            + '"'
-            + " "
-            + test_py_scripts.get_data_path("gcore")
-            + "byte.tif",
-        )
+def test_gdal_retile_png(script_path, tmp_path):
 
-        ds = gdal.Open(out_dirname + "/byte_1_1.png")
-        assert ds.GetRasterBand(1).Checksum() == 4672
-        ds = None
-        assert os.path.exists(out_dirname + "/byte_1_1.png.aux.xml")
-    finally:
-        shutil.rmtree(out_dirname)
+    out_dir = tmp_path / "outretile_png"
+    out_dir.mkdir()
+
+    test_py_scripts.run_py_script(
+        script_path,
+        "gdal_retile",
+        "-v -levels 2 -r bilinear -of PNG -targetDir "
+        + '"'
+        + str(out_dir)
+        + '"'
+        + " "
+        + test_py_scripts.get_data_path("gcore")
+        + "byte.tif",
+    )
 
+    with gdal.Open(str(out_dir / "byte_1_1.png")) as ds:
+        assert ds.GetRasterBand(1).Checksum() == 4672
 
-###############################################################################
-# Cleanup
-
-
-def test_gdal_retile_cleanup():
-
-    lst = [
-        "tmp/outretile/1/byte_1_1.tif",
-        "tmp/outretile/2/byte_1_1.tif",
-        "tmp/outretile/byte_1_1.tif",
-        "tmp/outretile/1",
-        "tmp/outretile/2",
-        "tmp/outretile",
-        "tmp/outretile2/1/rgba_1_1.tif",
-        "tmp/outretile2/2/rgba_1_1.tif",
-        "tmp/outretile2/1",
-        "tmp/outretile2/2",
-        "tmp/outretile2/rgba_1_1.tif",
-        "tmp/outretile2",
-        "tmp/in1.tif",
-        "tmp/in2.tif",
-        "tmp/outretile3/1/in1_1_1.tif",
-        "tmp/outretile3/2/in1_1_1.tif",
-        "tmp/outretile3/1",
-        "tmp/outretile3/2",
-        "tmp/outretile3/in1_1_1.tif",
-        "tmp/outretile3",
-        "tmp/in5.tif",
-    ]
-    for filename in lst:
-        try:
-            os.remove(filename)
-        except OSError:
-            try:
-                os.rmdir(filename)
-            except OSError:
-                pass
-
-    shutil.rmtree("tmp/outretile4")
-
-    if os.path.exists("tmp/outretile5"):
-        shutil.rmtree("tmp/outretile5")
+    assert os.path.exists(out_dir / "byte_1_1.png.aux.xml")
diff --git a/autotest/utilities/test_gdal_contour.py b/autotest/utilities/test_gdal_contour.py
index 10d1d8b37130..f195214565b7 100755
--- a/autotest/utilities/test_gdal_contour.py
+++ b/autotest/utilities/test_gdal_contour.py
@@ -29,7 +29,6 @@
 # DEALINGS IN THE SOFTWARE.
 ###############################################################################
 
-import os
 import struct
 
 import gdaltest
@@ -50,24 +49,9 @@ def gdal_contour_path():
     return test_cli_utilities.get_gdal_contour_path()
 
 
-###############################################################################
-# Test with -a and -i options
-
-
-def test_gdal_contour_1(gdal_contour_path):
-
-    try:
-        os.remove("tmp/contour.shp")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/contour.dbf")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/contour.shx")
-    except OSError:
-        pass
+@pytest.fixture()
+def testdata_tif(tmp_path):
+    tif_fname = str(tmp_path / "gdal_contour.tif")
 
     drv = gdal.GetDriverByName("GTiff")
     sr = osr.SpatialReference()
@@ -77,7 +61,7 @@ def test_gdal_contour_1(gdal_contour_path):
     size = 160
     precision = 1.0 / size
 
-    ds = drv.Create("tmp/gdal_contour.tif", size, size, 1)
+    ds = drv.Create(tif_fname, size, size, 1)
     ds.SetProjection(wkt)
     ds.SetGeoTransform([1, precision, 0, 50, 0, -precision])
 
@@ -119,12 +103,23 @@ def test_gdal_contour_1(gdal_contour_path):
 
     ds = None
 
+    yield tif_fname
+
+
+###############################################################################
+# Test with -a and -i options
+
+
+def test_gdal_contour_1(gdal_contour_path, testdata_tif, tmp_path):
+
+    contour_shp = str(tmp_path / "contour.shp")
+
     (_, err) = gdaltest.runexternal_out_and_err(
-        gdal_contour_path + " -a elev -i 10 tmp/gdal_contour.tif tmp/contour.shp"
+        gdal_contour_path + f" -a elev -i 10 {testdata_tif} {contour_shp}"
     )
     assert err is None or err == "", "got error/warning"
 
-    ds = ogr.Open("tmp/contour.shp")
+    ds = ogr.Open(contour_shp)
 
     expected_envelopes = [
         [1.25, 1.75, 49.25, 49.75],
@@ -134,10 +129,17 @@ def test_gdal_contour_1(gdal_contour_path):
 
     lyr = ds.ExecuteSQL("select * from contour order by elev asc")
 
-    assert lyr.GetSpatialRef().ExportToWkt() == wkt, "Did not get expected spatial ref"
+    raster_srs_wkt = gdal.Open(testdata_tif).GetSpatialRef().ExportToWkt()
+
+    assert (
+        lyr.GetSpatialRef().ExportToWkt() == raster_srs_wkt
+    ), "Did not get expected spatial ref"
 
     assert lyr.GetFeatureCount() == len(expected_envelopes)
 
+    size = 160
+    precision = 1.0 / size
+
     i = 0
     feat = lyr.GetNextFeature()
     while feat is not None:
@@ -163,31 +165,20 @@ def test_gdal_contour_1(gdal_contour_path):
 # Test with -fl option and -3d option
 
 
-def test_gdal_contour_2(gdal_contour_path):
+def test_gdal_contour_2(gdal_contour_path, testdata_tif, tmp_path):
 
-    try:
-        os.remove("tmp/contour.shp")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/contour.dbf")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/contour.shx")
-    except OSError:
-        pass
+    contour_shp = str(tmp_path / "contour.shp")
 
     # put -3d just after -fl to test #2793
-    gdaltest.runexternal(
-        gdal_contour_path
-        + " -a elev -fl 10 20 25 -3d tmp/gdal_contour.tif tmp/contour.shp"
+    _, err = gdaltest.runexternal_out_and_err(
+        gdal_contour_path + f" -a elev -fl 10 20 25 -3d {testdata_tif} {contour_shp}"
     )
+    assert not err
 
     size = 160
     precision = 1.0 / size
 
-    ds = ogr.Open("tmp/contour.shp")
+    ds = ogr.Open(contour_shp)
 
     expected_envelopes = [
         [1.25, 1.75, 49.25, 49.75],
@@ -231,27 +222,16 @@ def test_gdal_contour_2(gdal_contour_path):
 # Test on a real DEM
 
 
-def test_gdal_contour_3(gdal_contour_path):
+def test_gdal_contour_3(gdal_contour_path, tmp_path):
 
-    try:
-        os.remove("tmp/contour.shp")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/contour.dbf")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/contour.shx")
-    except OSError:
-        pass
+    contour_shp = str(tmp_path / "contour.shp")
 
     # put -3d just after -fl to test #2793
     gdaltest.runexternal(
-        gdal_contour_path + " -a elev -i 50 ../gdrivers/data/n43.tif tmp/contour.shp"
+        gdal_contour_path + f" -a elev -i 50 ../gdrivers/data/n43.tif {contour_shp}"
     )
 
-    ds = ogr.Open("tmp/contour.shp")
+    ds = ogr.Open(contour_shp)
 
     lyr = ds.ExecuteSQL("select distinct elev from contour order by elev asc")
 
@@ -273,20 +253,10 @@ def test_gdal_contour_3(gdal_contour_path):
 # Test contour orientation
 
 
-def test_gdal_contour_4(gdal_contour_path):
+def test_gdal_contour_4(gdal_contour_path, tmp_path):
 
-    try:
-        os.remove("tmp/contour_orientation.shp")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/contour_orientation.dbf")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/contour_orientation.shx")
-    except OSError:
-        pass
+    contour_orientation1_shp = str(tmp_path / "contour_orientation1.shp")
+    contour_orientation_tif = str(tmp_path / "contour_orientation.tif")
 
     drv = gdal.GetDriverByName("GTiff")
     wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]'
@@ -294,7 +264,7 @@ def test_gdal_contour_4(gdal_contour_path):
     size = 160
     precision = 1.0 / size
 
-    ds = drv.Create("tmp/gdal_contour_orientation.tif", size, size, 1)
+    ds = drv.Create(contour_orientation_tif, size, size, 1)
     ds.SetProjection(wkt)
     ds.SetGeoTransform([1, precision, 0, 50, 0, -precision])
 
@@ -335,10 +305,10 @@ def test_gdal_contour_4(gdal_contour_path):
 
     gdaltest.runexternal(
         gdal_contour_path
-        + " -a elev -i 10 tmp/gdal_contour_orientation.tif tmp/contour_orientation1.shp"
+        + f" -a elev -i 10 {contour_orientation_tif} {contour_orientation1_shp}"
     )
 
-    ds = ogr.Open("tmp/contour_orientation1.shp")
+    ds = ogr.Open(contour_orientation1_shp)
 
     expected_contours = [
         "LINESTRING ("
@@ -385,16 +355,18 @@ def test_gdal_contour_4(gdal_contour_path):
 # Test contour orientation
 
 
-def test_gdal_contour_5(gdal_contour_path):
+def test_gdal_contour_5(gdal_contour_path, tmp_path):
 
     ds = None
 
+    contour_orientation2_shp = str(tmp_path / "contour_orientation2.shp")
+
     gdaltest.runexternal(
         gdal_contour_path
-        + " -a elev -i 10 data/contour_orientation.tif tmp/contour_orientation2.shp"
+        + f" -a elev -i 10 data/contour_orientation.tif {contour_orientation2_shp}"
     )
 
-    ds = ogr.Open("tmp/contour_orientation2.shp")
+    ds = ogr.Open(contour_orientation2_shp)
 
     expected_contours = [
         "LINESTRING (0.0 1.999999,"
@@ -423,23 +395,3 @@ def test_gdal_contour_5(gdal_contour_path):
             feat = lyr.GetNextFeature()
 
     ds.Destroy()
-
-
-###############################################################################
-# Cleanup
-
-
-def test_gdal_contour_cleanup():
-
-    ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource("tmp/contour.shp")
-    ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource(
-        "tmp/contour_orientation1.shp"
-    )
-    ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource(
-        "tmp/contour_orientation2.shp"
-    )
-    try:
-        os.remove("tmp/gdal_contour.tif")
-        os.remove("tmp/gdal_contour_orientation.tif")
-    except OSError:
-        pass
diff --git a/autotest/utilities/test_gdal_grid.py b/autotest/utilities/test_gdal_grid.py
index 700832ad436a..554e95395855 100755
--- a/autotest/utilities/test_gdal_grid.py
+++ b/autotest/utilities/test_gdal_grid.py
@@ -29,7 +29,6 @@
 # DEALINGS IN THE SOFTWARE.
 ###############################################################################
 
-import os
 import struct
 import sys
 
@@ -52,87 +51,54 @@ def gdal_grid_path():
     return test_cli_utilities.get_gdal_grid_path()
 
 
-# List of output TIFF files that will be created by tests and later deleted
-# in test_gdal_grid_cleanup()
-outfiles = []
-
-###############################################################################
-@pytest.fixture(autouse=True, scope="module")
-def auto_cleanup():
-
-    yield
-
-    if os.path.exists("tmp/n43.shp"):
-        ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource("tmp/n43.shp")
-    drv = gdal.GetDriverByName("GTiff")
-    for outfile in outfiles:
-        if os.path.exists(outfile):
-            drv.Delete(outfile)
+@pytest.fixture()
+def n43_shp(tmp_path):
+    n43_shp_fname = str(tmp_path / "n43.shp")
 
+    # Create an OGR grid from the values of n43.tif
+    ds = gdal.Open("../gdrivers/data/n43.tif")
+    geotransform = ds.GetGeoTransform()
 
-###############################################################################
-#
+    shape_drv = ogr.GetDriverByName("ESRI Shapefile")
+    with shape_drv.CreateDataSource(str(tmp_path)) as shape_ds:
+        shape_lyr = shape_ds.CreateLayer("n43")
 
+        data = ds.ReadRaster(0, 0, 121, 121)
+        array_val = struct.unpack("h" * 121 * 121, data)
+        for j in range(121):
+            for i in range(121):
+                wkt = "POINT(%f %f %s)" % (
+                    geotransform[0] + (i + 0.5) * geotransform[1],
+                    geotransform[3] + (j + 0.5) * geotransform[5],
+                    array_val[j * 121 + i],
+                )
+                dst_feat = ogr.Feature(feature_def=shape_lyr.GetLayerDefn())
+                dst_feat.SetGeometry(ogr.CreateGeometryFromWkt(wkt))
+                shape_lyr.CreateFeature(dst_feat)
 
-def test_gdal_grid_1(gdal_grid_path):
+        shape_ds.ExecuteSQL("CREATE SPATIAL INDEX ON n43")
 
-    shape_drv = ogr.GetDriverByName("ESRI Shapefile")
-    outfiles.append("tmp/n43.tif")
-
-    try:
-        os.remove("tmp/n43.shp")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/n43.dbf")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/n43.shx")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/n43.qix")
-    except OSError:
-        pass
+    yield n43_shp_fname
 
-    # Create an OGR grid from the values of n43.tif
-    ds = gdal.Open("../gdrivers/data/n43.tif")
-    geotransform = ds.GetGeoTransform()
 
-    shape_drv = ogr.GetDriverByName("ESRI Shapefile")
-    shape_ds = shape_drv.CreateDataSource("tmp")
-    shape_lyr = shape_ds.CreateLayer("n43")
-
-    data = ds.ReadRaster(0, 0, 121, 121)
-    array_val = struct.unpack("h" * 121 * 121, data)
-    for j in range(121):
-        for i in range(121):
-            wkt = "POINT(%f %f %s)" % (
-                geotransform[0] + (i + 0.5) * geotransform[1],
-                geotransform[3] + (j + 0.5) * geotransform[5],
-                array_val[j * 121 + i],
-            )
-            dst_feat = ogr.Feature(feature_def=shape_lyr.GetLayerDefn())
-            dst_feat.SetGeometry(ogr.CreateGeometryFromWkt(wkt))
-            shape_lyr.CreateFeature(dst_feat)
+###############################################################################
+#
 
-    dst_feat.Destroy()
 
-    shape_ds.ExecuteSQL("CREATE SPATIAL INDEX ON n43")
+def test_gdal_grid_1(gdal_grid_path, n43_shp, tmp_path):
 
-    shape_ds.Destroy()
+    output_tif = str(tmp_path / "n43.tif")
 
     # Create a GDAL dataset from the previous generated OGR grid
     (_, err) = gdaltest.runexternal_out_and_err(
         gdal_grid_path
-        + " -txe -80.0041667 -78.9958333 -tye 42.9958333 44.0041667 -outsize 121 121 -ot Int16 -a nearest:radius1=0.0:radius2=0.0:angle=0.0 -co TILED=YES -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 tmp/n43.shp "
-        + outfiles[-1]
+        + f" -txe -80.0041667 -78.9958333 -tye 42.9958333 44.0041667 -outsize 121 121 -ot Int16 -a nearest:radius1=0.0:radius2=0.0:angle=0.0 -co TILED=YES -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 {n43_shp} {output_tif}"
     )
     assert err is None or err == "", "got error/warning"
 
     # We should get the same values as in n43.td0
-    ds2 = gdal.Open(outfiles[-1])
+    ds = gdal.Open("../gdrivers/data/n43.tif")
+    ds2 = gdal.Open(output_tif)
     assert (
         ds.GetRasterBand(1).Checksum() == ds2.GetRasterBand(1).Checksum()
     ), "bad checksum : got %d, expected %d" % (
@@ -150,30 +116,24 @@ def test_gdal_grid_1(gdal_grid_path):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_2(gdal_grid_path):
+def test_gdal_grid_2(gdal_grid_path, tmp_path):
 
     # Open reference dataset
     ds_ref = gdal.Open("../gcore/data/byte.tif")
     checksum_ref = ds_ref.GetRasterBand(1).Checksum()
     ds_ref = None
 
-    #################
-    outfiles.append("tmp/grid_near.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    grid_near_tif = str(tmp_path / "grid_near.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Grid nodes are located exactly in raster nodes.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 data/grid.vrt {grid_near_tif}"
     )
 
     # We should get the same values as in "gcore/data/byte.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(grid_near_tif)
     if ds.GetRasterBand(1).Checksum() != checksum_ref:
         print(
             "bad checksum : got %d, expected %d"
@@ -183,23 +143,26 @@ def test_gdal_grid_2(gdal_grid_path):
     assert ds.GetRasterBand(1).GetNoDataValue() == 0.0, "expected a nodata value"
     ds = None
 
-    #################
-    outfiles.append("tmp/grid_near_shift.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+
+@pytest.mark.require_driver("CSV")
+def test_gdal_grid_2bis(gdal_grid_path, tmp_path):
+
+    # Open reference dataset
+    ds_ref = gdal.Open("../gcore/data/byte.tif")
+    checksum_ref = ds_ref.GetRasterBand(1).Checksum()
+    ds_ref = None
+
+    grid_near_shift_tif = str(tmp_path / "grid_near_shift.tif")
 
     # Now the same, but shift grid nodes a bit in both horizontal and vertical
     # directions.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 data/grid.vrt {grid_near_shift_tif}"
     )
 
     # We should get the same values as in "gcore/data/byte.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(grid_near_shift_tif)
     if ds.GetRasterBand(1).Checksum() != checksum_ref:
         print(
             "bad checksum : got %d, expected %d"
@@ -211,19 +174,14 @@ def test_gdal_grid_2(gdal_grid_path):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.parametrize("use_quadtree", [True, False])
-def test_gdal_grid_3(gdal_grid_path, use_quadtree):
+def test_gdal_grid_3(gdal_grid_path, tmp_path, use_quadtree):
 
     # Open reference dataset
     ds_ref = gdal.Open("../gcore/data/byte.tif")
     checksum_ref = ds_ref.GetRasterBand(1).Checksum()
     ds_ref = None
 
-    #################
-    outfiles.append("tmp/grid_near_search3.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    grid_near_search3_tif = str(tmp_path / "grid_near_search3.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Try the search ellipse larger than the raster cell.
@@ -232,12 +190,11 @@ def test_gdal_grid_3(gdal_grid_path, use_quadtree):
     ):
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=180.0:radius2=180.0:angle=0.0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=180.0:radius2=180.0:angle=0.0:nodata=0.0 data/grid.vrt {grid_near_search3_tif}"
         )
 
         # We should get the same values as in "gcore/data/byte.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(grid_near_search3_tif)
         if ds.GetRasterBand(1).Checksum() != checksum_ref:
             print(
                 "bad checksum : got %d, expected %d"
@@ -247,21 +204,16 @@ def test_gdal_grid_3(gdal_grid_path, use_quadtree):
         ds = None
 
         #################
-        outfiles.append("tmp/grid_near_search1.tif")
-        try:
-            os.remove(outfiles[-1])
-        except OSError:
-            pass
+        grid_near_search1_tif = str(tmp_path / "grid_near_search1.tif")
 
         # Search ellipse smaller than the raster cell.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=20.0:radius2=20.0:angle=0.0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=20.0:radius2=20.0:angle=0.0:nodata=0.0 data/grid.vrt {grid_near_search1_tif}"
         )
 
         # We should get the same values as in "gcore/data/byte.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(grid_near_search1_tif)
         if ds.GetRasterBand(1).Checksum() != checksum_ref:
             print(
                 "bad checksum : got %d, expected %d"
@@ -271,21 +223,16 @@ def test_gdal_grid_3(gdal_grid_path, use_quadtree):
         ds = None
 
         #################
-        outfiles.append("tmp/grid_near_shift_search3.tif")
-        try:
-            os.remove(outfiles[-1])
-        except OSError:
-            pass
+        grid_near_shift_search3_tif = str(tmp_path / "grid_near_shift_search3.tif")
 
         # Large search ellipse and the grid shift.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=180.0:radius2=180.0:angle=0.0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=180.0:radius2=180.0:angle=0.0:nodata=0.0 data/grid.vrt {grid_near_shift_search3_tif}"
         )
 
         # We should get the same values as in "gcore/data/byte.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(grid_near_shift_search3_tif)
         if ds.GetRasterBand(1).Checksum() != checksum_ref:
             print(
                 "bad checksum : got %d, expected %d"
@@ -295,21 +242,16 @@ def test_gdal_grid_3(gdal_grid_path, use_quadtree):
         ds = None
 
         #################
-        outfiles.append("tmp/grid_near_shift_search1.tif")
-        try:
-            os.remove(outfiles[-1])
-        except OSError:
-            pass
+        grid_near_shift_search1_tif = str(tmp_path / "grid_near_shift_search1.tif")
 
         # Small search ellipse and the grid shift.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=20.0:radius2=20.0:angle=0.0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=20.0:radius2=20.0:angle=0.0:nodata=0.0 data/grid.vrt {grid_near_shift_search1_tif}"
         )
 
         # We should get the same values as in "gcore/data/byte.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(grid_near_shift_search1_tif)
         if ds.GetRasterBand(1).Checksum() != checksum_ref:
             print(
                 "bad checksum : got %d, expected %d"
@@ -324,170 +266,55 @@ def test_gdal_grid_3(gdal_grid_path, use_quadtree):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_4(gdal_grid_path):
+@pytest.mark.parametrize(
+    "algorithm,threads",
+    [("Generic", None), ("SSE", None), ("AVX", None), ("AVX", 1), ("AVX", 2)],
+)
+def test_gdal_grid_4(gdal_grid_path, algorithm, threads, tmp_path):
 
-    #################
-    # Test generic implementation (no AVX, no SSE)
-    outfiles.append("tmp/grid_invdist_generic.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    alg_flags = {
+        "Generic": "--config GDAL_USE_AVX NO --config GDAL_USE_SSE NO",
+        "SSE": "--config GDAL_USE_AVX NO",
+        "AVX": "",
+    }
+
+    output_tif = str(tmp_path / "output.tif")
+
+    flags = alg_flags[algorithm]
+
+    if threads:
+        flags += " --config GDAL_NUM_THREADS {threads}"
 
     # Create a GDAL dataset from the values of "grid.csv".
-    print("Step 1: Disabling AVX/SSE optimized versions...")
     (_, err) = gdaltest.runexternal_out_and_err(
         gdal_grid_path
-        + " --debug on --config GDAL_USE_AVX NO --config GDAL_USE_SSE NO -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" --debug on {flags} -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
+
     pos = err.find(" threads")
     if pos >= 0:
         pos_blank = err[0 : pos - 1].rfind(" ")
         if pos_blank >= 0:
             print("Step 1: %s threads used" % err[pos_blank + 1 : pos])
 
-    # We should get the same values as in "ref_data/gdal_invdist.tif"
-    ds = gdal.Open(outfiles[-1])
-    ds_ref = gdal.Open("ref_data/grid_invdist.tif")
-    maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
-    if maxdiff > 1:
-        gdaltest.compare_ds(ds, ds_ref, verbose=1)
-        pytest.fail("Image too different from the reference")
-    ds_ref = None
-    ds = None
+    if algorithm == "SSE" and "SSE" not in err:
+        pytest.skip(f"{algorithm} not used")
 
-    #################
-    # Potentially test optimized SSE implementation
-
-    outfiles.append("tmp/grid_invdist_sse.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
-
-    # Create a GDAL dataset from the values of "grid.csv".
-    print("Step 2: Trying SSE optimized version...")
-    (_, err) = gdaltest.runexternal_out_and_err(
-        gdal_grid_path
-        + " --debug on --config GDAL_USE_AVX NO -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
-    )
-    if "SSE" in err:
-        print("...SSE optimized version used")
-    else:
-        print("...SSE optimized version NOT used")
-
-    # We should get the same values as in "ref_data/gdal_invdist.tif"
-    ds = gdal.Open(outfiles[-1])
-    ds_ref = gdal.Open("ref_data/grid_invdist.tif")
-    maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
-    if maxdiff > 1:
-        gdaltest.compare_ds(ds, ds_ref, verbose=1)
-        pytest.fail("Image too different from the reference")
-    ds_ref = None
-    ds = None
-
-    #################
-    # Potentially test optimized AVX implementation
-
-    outfiles.append("tmp/grid_invdist_avx.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
 
-    # Create a GDAL dataset from the values of "grid.csv".
-    print("Step 3: Trying AVX optimized version...")
-    (_, err) = gdaltest.runexternal_out_and_err(
-        gdal_grid_path
-        + " --debug on -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
-    )
-    if "AVX" in err:
-        print("...AVX optimized version used")
-    else:
-        print("...AVX optimized version NOT used")
-
-    # We should get the same values as in "ref_data/gdal_invdist.tif"
-    ds = gdal.Open(outfiles[-1])
-    ds_ref = gdal.Open("ref_data/grid_invdist.tif")
-    maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
-    if maxdiff > 1:
-        gdaltest.compare_ds(ds, ds_ref, verbose=1)
-        pytest.fail("Image too different from the reference")
-    ds_ref = None
-    ds = None
-
-    #################
-    # Test GDAL_NUM_THREADS config option to 1
-
-    outfiles.append("tmp/grid_invdist_1thread.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
-
-    # Create a GDAL dataset from the values of "grid.csv".
-    gdaltest.runexternal(
-        gdal_grid_path
-        + " --config GDAL_NUM_THREADS 1 -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
-    )
-
-    # We should get the same values as in "ref_data/gdal_invdist.tif"
-    ds = gdal.Open(outfiles[-1])
-    ds_ref = gdal.Open("ref_data/grid_invdist.tif")
-    maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
-    if maxdiff > 1:
-        gdaltest.compare_ds(ds, ds_ref, verbose=1)
-        pytest.fail("Image too different from the reference")
-    ds_ref = None
-    ds = None
-
-    #################
-    # Test GDAL_NUM_THREADS config option to 2
-
-    outfiles.append("tmp/grid_invdist_2threads.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
-
-    # Create a GDAL dataset from the values of "grid.csv".
-    gdaltest.runexternal(
-        gdal_grid_path
-        + " --config GDAL_NUM_THREADS 2 -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:angle=0.0:max_points=0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
-    )
-
-    # We should get the same values as in "ref_data/gdal_invdist.tif"
-    ds = gdal.Open(outfiles[-1])
-    ds_ref = gdal.Open("ref_data/grid_invdist.tif")
-    maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
-    if maxdiff > 1:
-        gdaltest.compare_ds(ds, ds_ref, verbose=1)
-        pytest.fail("Image too different from the reference")
-    ds_ref = None
-    ds = None
+@pytest.mark.require_driver("CSV")
+def test_gdal_grid_4bis(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_invdist_90_90_8p.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_invdist_90_90_8p.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Circular window, shifted, test min points and NODATA setting.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdist:power=2.0:radius1=90.0:radius2=90.0:angle=0.0:max_points=0:min_points=8:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdist:power=2.0:radius1=90.0:radius2=90.0:angle=0.0:max_points=0:min_points=8:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_invdist_90_90_8p.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_invdist_90_90_8p.tif")
     maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
     if maxdiff > 1:
@@ -502,26 +329,20 @@ def test_gdal_grid_4(gdal_grid_path):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_5(gdal_grid_path):
+def test_gdal_grid_5(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_average.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_average.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # We are using all the points from input dataset to average, so
     # the result is a raster filled with the same value in each node.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_average.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_average.tif")
     maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
     ds_ref = None
@@ -530,23 +351,21 @@ def test_gdal_grid_5(gdal_grid_path):
         pytest.fail("Image too different from the reference")
     ds = None
 
-    #################
-    outfiles.append("tmp/grid_average_300_100_40.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+
+@pytest.mark.require_driver("CSV")
+def test_gdal_grid_5bis(gdal_grid_path, tmp_path):
+
+    output_tif = str(tmp_path / "grid_average_300_100_40.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Elliptical window, rotated.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average:radius1=300.0:radius2=100.0:angle=40.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average:radius1=300.0:radius2=100.0:angle=40.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_average_300_100_40.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_average_300_100_40.tif")
     maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
     ds_ref = None
@@ -558,14 +377,9 @@ def test_gdal_grid_5(gdal_grid_path):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.parametrize("use_quadtree", [True, False])
-def test_gdal_grid_6(gdal_grid_path, use_quadtree):
+def test_gdal_grid_6(gdal_grid_path, tmp_path, use_quadtree):
 
-    #################
-    outfiles.append("tmp/grid_average_190_190.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_average_190_190.tif")
 
     with gdaltest.config_option(
         "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
@@ -574,12 +388,11 @@ def test_gdal_grid_6(gdal_grid_path, use_quadtree):
         # This time using a circular window.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average:radius1=190.0:radius2=190.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average:radius1=190.0:radius2=190.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_average_190_190.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_average_190_190.tif")
         maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
         ds_ref = None
@@ -588,23 +401,26 @@ def test_gdal_grid_6(gdal_grid_path, use_quadtree):
             pytest.fail("Image too different from the reference")
         ds = None
 
-        #################
-        outfiles.append("tmp/grid_average_90_90_8p.tif")
-        try:
-            os.remove(outfiles[-1])
-        except OSError:
-            pass
+
+@pytest.mark.require_driver("CSV")
+@pytest.mark.parametrize("use_quadtree", [True, False])
+def test_gdal_grid_6bis(gdal_grid_path, tmp_path, use_quadtree):
+
+    output_tif = str(tmp_path / "grid_average_90_90_8p.tif")
+
+    with gdaltest.config_option(
+        "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
+    ):
 
         # Create a GDAL dataset from the values of "grid.csv".
         # Circular window, shifted, test min points and NODATA setting.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average:radius1=90.0:radius2=90.0:angle=0.0:min_points=8:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average:radius1=90.0:radius2=90.0:angle=0.0:min_points=8:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_average_90_90_8p.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_average_90_90_8p.tif")
         maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
         ds_ref = None
@@ -619,25 +435,19 @@ def test_gdal_grid_6(gdal_grid_path, use_quadtree):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_7(gdal_grid_path):
+def test_gdal_grid_7(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_minimum.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_minimum.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Search the whole dataset for minimum.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a minimum:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a minimum:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_minimum.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_minimum.tif")
     assert (
         ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -648,23 +458,21 @@ def test_gdal_grid_7(gdal_grid_path):
     ds_ref = None
     ds = None
 
-    #################
-    outfiles.append("tmp/grid_minimum_400_100_120.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+
+@pytest.mark.require_driver("CSV")
+def test_gdal_grid_7bis(gdal_grid_path, tmp_path):
+
+    output_tif = str(tmp_path / "grid_minimum_400_100_120.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Elliptical window, rotated.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a minimum:radius1=400.0:radius2=100.0:angle=120.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a minimum:radius1=400.0:radius2=100.0:angle=120.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_minimum_400_100_120.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_minimum_400_100_120.tif")
     assert (
         ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -678,14 +486,9 @@ def test_gdal_grid_7(gdal_grid_path):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.parametrize("use_quadtree", [True, False])
-def test_gdal_grid_8(gdal_grid_path, use_quadtree):
+def test_gdal_grid_8(gdal_grid_path, tmp_path, use_quadtree):
 
-    #################
-    outfiles.append("tmp/grid_minimum_180_180.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_minimum_180_180.tif")
 
     with gdaltest.config_option(
         "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
@@ -694,12 +497,11 @@ def test_gdal_grid_8(gdal_grid_path, use_quadtree):
         # Search ellipse larger than the raster cell.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a minimum:radius1=180.0:radius2=180.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a minimum:radius1=180.0:radius2=180.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_minimum_180_180.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_minimum_180_180.tif")
         assert (
             ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -710,23 +512,26 @@ def test_gdal_grid_8(gdal_grid_path, use_quadtree):
         ds_ref = None
         ds = None
 
-        #################
-        outfiles.append("tmp/grid_minimum_20_20.tif")
-        try:
-            os.remove(outfiles[-1])
-        except OSError:
-            pass
+
+@pytest.mark.require_driver("CSV")
+@pytest.mark.parametrize("use_quadtree", [True, False])
+def test_gdal_grid_8bis(gdal_grid_path, tmp_path, use_quadtree):
+
+    output_tif = str(tmp_path / "grid_minimum_20_20.tif")
+
+    with gdaltest.config_option(
+        "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
+    ):
 
         # Create a GDAL dataset from the values of "grid.csv".
         # Search ellipse smaller than the raster cell.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a minimum:radius1=20.0:radius2=20.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a minimum:radius1=20.0:radius2=20.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_minimum_20_20.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_minimum_20_20.tif")
         assert (
             ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -743,25 +548,19 @@ def test_gdal_grid_8(gdal_grid_path, use_quadtree):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_9(gdal_grid_path):
+def test_gdal_grid_9(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_maximum.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_maximum.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Search the whole dataset for maximum.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a maximum:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a maximum:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_maximum.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_maximum.tif")
     assert (
         ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -772,23 +571,21 @@ def test_gdal_grid_9(gdal_grid_path):
     ds_ref = None
     ds = None
 
-    #################
-    outfiles.append("tmp/grid_maximum_100_100.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+
+@pytest.mark.require_driver("CSV")
+def test_gdal_grid_9bis(gdal_grid_path, tmp_path):
+
+    output_tif = str(tmp_path / "grid_maximum_100_100.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Circular window.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a maximum:radius1=100.0:radius2=100.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a maximum:radius1=100.0:radius2=100.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_maximum_100_100.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_maximum_100_100.tif")
     assert (
         ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -802,14 +599,9 @@ def test_gdal_grid_9(gdal_grid_path):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.parametrize("use_quadtree", [True, False])
-def test_gdal_grid_10(gdal_grid_path, use_quadtree):
+def test_gdal_grid_10bis(gdal_grid_path, tmp_path, use_quadtree):
 
-    #################
-    outfiles.append("tmp/grid_maximum_180_180.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_maximum_180_180.tif")
 
     with gdaltest.config_option(
         "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
@@ -818,12 +610,11 @@ def test_gdal_grid_10(gdal_grid_path, use_quadtree):
         # Search ellipse larger than the raster cell.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a maximum:radius1=180.0:radius2=180.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a maximum:radius1=180.0:radius2=180.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_maximum_180_180.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_maximum_180_180.tif")
         assert (
             ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -834,23 +625,26 @@ def test_gdal_grid_10(gdal_grid_path, use_quadtree):
         ds_ref = None
         ds = None
 
-        #################
-        outfiles.append("tmp/grid_maximum_20_20.tif")
-        try:
-            os.remove(outfiles[-1])
-        except OSError:
-            pass
+
+@pytest.mark.require_driver("CSV")
+@pytest.mark.parametrize("use_quadtree", [True, False])
+def test_gdal_grid_10(gdal_grid_path, tmp_path, use_quadtree):
+
+    output_tif = str(tmp_path / "grid_maximum_20_20.tif")
+
+    with gdaltest.config_option(
+        "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
+    ):
 
         # Create a GDAL dataset from the values of "grid.csv".
         # Search ellipse smaller than the raster cell.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a maximum:radius1=20.0:radius2=20.0:angle=120.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a maximum:radius1=20.0:radius2=20.0:angle=120.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_maximum_20_20.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_maximum_20_20.tif")
         assert (
             ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -867,25 +661,19 @@ def test_gdal_grid_10(gdal_grid_path, use_quadtree):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_11(gdal_grid_path):
+def test_gdal_grid_11(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_range.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_range.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Search the whole dataset.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a range:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a range:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_range.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_range.tif")
     assert (
         ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -899,14 +687,9 @@ def test_gdal_grid_11(gdal_grid_path):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.parametrize("use_quadtree", [True, False])
-def test_gdal_grid_12(gdal_grid_path, use_quadtree):
+def test_gdal_grid_12(gdal_grid_path, tmp_path, use_quadtree):
 
-    #################
-    outfiles.append("tmp/grid_range_90_90_8p.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_range_90_90_8p.tif")
 
     with gdaltest.config_option(
         "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
@@ -916,12 +699,11 @@ def test_gdal_grid_12(gdal_grid_path, use_quadtree):
         # points found.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a range:radius1=90.0:radius2=90.0:angle=0.0:min_points=8:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a range:radius1=90.0:radius2=90.0:angle=0.0:min_points=8:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_range_90_90_8p.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_range_90_90_8p.tif")
         assert (
             ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -939,14 +721,9 @@ def test_gdal_grid_12(gdal_grid_path, use_quadtree):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.parametrize("use_quadtree", [True, False])
-def test_gdal_grid_13(gdal_grid_path, use_quadtree):
+def test_gdal_grid_13bis(gdal_grid_path, tmp_path, use_quadtree):
 
-    #################
-    outfiles.append("tmp/grid_count_70_70.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_count_70_70.tif")
 
     with gdaltest.config_option(
         "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
@@ -954,12 +731,11 @@ def test_gdal_grid_13(gdal_grid_path, use_quadtree):
         # Create a GDAL dataset from the values of "grid.csv".
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a count:radius1=70.0:radius2=70.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a count:radius1=70.0:radius2=70.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_count_70_70.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_count_70_70.tif")
         assert (
             ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -970,22 +746,24 @@ def test_gdal_grid_13(gdal_grid_path, use_quadtree):
         ds_ref = None
         ds = None
 
-        #################
-        outfiles.append("tmp/grid_count_300_300.tif")
-        try:
-            os.remove(outfiles[-1])
-        except OSError:
-            pass
 
+@pytest.mark.require_driver("CSV")
+@pytest.mark.parametrize("use_quadtree", [True, False])
+def test_gdal_grid_13(gdal_grid_path, tmp_path, use_quadtree):
+
+    output_tif = str(tmp_path / "grid_count_300_300.tif")
+
+    with gdaltest.config_option(
+        "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
+    ):
         # Create a GDAL dataset from the values of "grid.csv".
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a count:radius1=300.0:radius2=300.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a count:radius1=300.0:radius2=300.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_count_300_300.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_count_300_300.tif")
         assert (
             ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
@@ -1002,26 +780,20 @@ def test_gdal_grid_13(gdal_grid_path, use_quadtree):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_14(gdal_grid_path):
+def test_gdal_grid_14(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_avdist.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_avdist.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # We are using all the points from input dataset to average, so
     # the result is a raster filled with the same value in each node.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average_distance:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average_distance:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_avdist.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_avdist.tif")
     maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
     if maxdiff > 1:
@@ -1033,14 +805,9 @@ def test_gdal_grid_14(gdal_grid_path):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.parametrize("use_quadtree", [True, False])
-def test_gdal_grid_15(gdal_grid_path, use_quadtree):
+def test_gdal_grid_15(gdal_grid_path, tmp_path, use_quadtree):
 
-    #################
-    outfiles.append("tmp/grid_avdist_150_150.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_avdist_150_150.tif")
 
     with gdaltest.config_option(
         "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
@@ -1050,12 +817,11 @@ def test_gdal_grid_15(gdal_grid_path, use_quadtree):
         # the result is a raster filled with the same value in each node.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average_distance:radius1=150.0:radius2=150.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average_distance:radius1=150.0:radius2=150.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_avdist_150_150.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_avdist_150_150.tif")
         maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
         if maxdiff > 1:
@@ -1070,26 +836,20 @@ def test_gdal_grid_15(gdal_grid_path, use_quadtree):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_16(gdal_grid_path):
+def test_gdal_grid_16(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_avdistpts_150_50_-15.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_avdistpts_150_50_-15.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # We are using all the points from input dataset to average, so
     # the result is a raster filled with the same value in each node.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average_distance_pts:radius1=150.0:radius2=50.0:angle=-15.0:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average_distance_pts:radius1=150.0:radius2=50.0:angle=-15.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_avdistpts_150_50_-15.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_avdistpts_150_50_-15.tif")
     maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
     if maxdiff > 1:
@@ -1101,14 +861,9 @@ def test_gdal_grid_16(gdal_grid_path):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.parametrize("use_quadtree", [True, False])
-def test_gdal_grid_17(gdal_grid_path, use_quadtree):
+def test_gdal_grid_17(gdal_grid_path, tmp_path, use_quadtree):
 
-    #################
-    outfiles.append("tmp/grid_avdistpts_150_150.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_avdistpts_150_150.tif")
 
     with gdaltest.config_option(
         "GDAL_GRID_POINT_COUNT_THRESHOLD", "0" if use_quadtree else "1000000000"
@@ -1118,12 +873,11 @@ def test_gdal_grid_17(gdal_grid_path, use_quadtree):
         # the result is a raster filled with the same value in each node.
         gdaltest.runexternal(
             gdal_grid_path
-            + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average_distance_pts:radius1=150.0:radius2=150.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt "
-            + outfiles[-1]
+            + f" -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a average_distance_pts:radius1=150.0:radius2=150.0:angle=0.0:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
         )
 
         # We should get the same values as in "ref_data/grid_avdistpts_150_150.tif"
-        ds = gdal.Open(outfiles[-1])
+        ds = gdal.Open(output_tif)
         ds_ref = gdal.Open("ref_data/grid_avdistpts_150_150.tif")
         maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
         if maxdiff > 1:
@@ -1138,21 +892,20 @@ def test_gdal_grid_17(gdal_grid_path, use_quadtree):
 
 
 @pytest.mark.skipif(not gdal.HasTriangulation(), reason="qhull missing")
-def test_gdal_grid_18(gdal_grid_path):
+def test_gdal_grid_18(gdal_grid_path, tmp_path, n43_shp):
 
-    outfiles.append("tmp/n43_linear.tif")
+    output_tif = str(tmp_path / "n43_linear.tif")
 
     # Create a GDAL dataset from the previous generated OGR grid
     (_, err) = gdaltest.runexternal_out_and_err(
         gdal_grid_path
-        + " -txe -80.0041667 -78.9958333 -tye 42.9958333 44.0041667 -outsize 121 121 -ot Int16 -l n43 -a linear -co TILED=YES -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 tmp/n43.shp "
-        + outfiles[-1]
+        + f" -txe -80.0041667 -78.9958333 -tye 42.9958333 44.0041667 -outsize 121 121 -ot Int16 -l n43 -a linear -co TILED=YES -co BLOCKXSIZE=256 -co BLOCKYSIZE=256 {n43_shp} {output_tif}"
     )
     assert err is None or err == "", "got error/warning"
 
     # We should get the same values as in n43.tif
     ds = gdal.Open("../gdrivers/data/n43.tif")
-    ds2 = gdal.Open(outfiles[-1])
+    ds2 = gdal.Open(output_tif)
     assert (
         ds.GetRasterBand(1).Checksum() == ds2.GetRasterBand(1).Checksum()
     ), "bad checksum : got %d, expected %d" % (
@@ -1169,25 +922,20 @@ def test_gdal_grid_18(gdal_grid_path):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_19(gdal_grid_path):
+def test_gdal_grid_19(gdal_grid_path, tmp_path):
 
+    output_tif = str(tmp_path / "grid_invdistnn_generic.tif")
     #################
     # Test generic implementation (no AVX, no SSE)
-    outfiles.append("tmp/grid_invdistnn_generic.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
 
     # Create a GDAL dataset from the values of "grid.csv".
     (_, _) = gdaltest.runexternal_out_and_err(
         gdal_grid_path
-        + " -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdistnn:power=2.0:radius=1.0:max_points=12:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdistnn:power=2.0:radius=1.0:max_points=12:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/gdal_invdistnn.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_invdistnn.tif")
     maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
     if maxdiff > 0.00001:
@@ -1196,23 +944,21 @@ def test_gdal_grid_19(gdal_grid_path):
     ds_ref = None
     ds = None
 
-    #################
-    outfiles.append("tmp/grid_invdistnn_250_8minp.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+
+@pytest.mark.require_driver("CSV")
+def test_gdal_grid_19_250_8minp(gdal_grid_path, tmp_path):
+
+    output_tif = str(tmp_path / "grid_invdistnn_250_8minp.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     # Circular window, shifted, test min points and NODATA setting.
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdistnn:power=2.0:radius=250.0:max_points=12:min_points=8:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdistnn:power=2.0:radius=250.0:max_points=12:min_points=8:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/grid_invdistnn_250_8minp.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_invdistnn_250_8minp.tif")
     maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
     if maxdiff > 0.00001:
@@ -1221,23 +967,23 @@ def test_gdal_grid_19(gdal_grid_path):
     ds_ref = None
     ds = None
 
+
+@pytest.mark.require_driver("CSV")
+def test_gdal_grid_19_250_10maxp_3pow(gdal_grid_path, tmp_path):
+
+    output_tif = str(tmp_path / "grid_invdistnn_250_10maxp_3pow.tif")
+
     #################
     # Test generic implementation with max_points and radius specified
-    outfiles.append("tmp/grid_invdistnn_250_10maxp_3pow.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
 
     # Create a GDAL dataset from the values of "grid.csv".
     gdaltest.runexternal(
         gdal_grid_path
-        + " -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdistnn:power=3.0:radius=250.0:max_points=10:min_points=0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -txe 440721.0 441920.0 -tye 3751321.0 3750120.0 -outsize 20 20 -ot Float64 -l grid -a invdistnn:power=3.0:radius=250.0:max_points=10:min_points=0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
     # We should get the same values as in "ref_data/gdal_invdistnn_250_10maxp_3pow.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_invdistnn_250_10maxp_3pow.tif")
     maxdiff = gdaltest.compare_ds(ds, ds_ref, verbose=0)
     if maxdiff > 0.00001:
@@ -1253,16 +999,12 @@ def test_gdal_grid_19(gdal_grid_path):
 
 @pytest.mark.require_driver("CSV")
 @pytest.mark.require_geos
-def test_gdal_grid_clipsrc(gdal_grid_path):
+def test_gdal_grid_clipsrc(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_clipsrc.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_clipsrc.tif")
+    clip_csv = str(tmp_path / "clip.csv")
 
-    open("tmp/clip.csv", "wt").write(
+    open(clip_csv, "wt").write(
         'id,WKT\n1,"POLYGON((440750 3751340,440750 3750100,441900 3750100,441900 3751340,440750 3751340))"\n'
     )
 
@@ -1270,14 +1012,11 @@ def test_gdal_grid_clipsrc(gdal_grid_path):
     # Grid nodes are located exactly in raster nodes.
     gdaltest.runexternal_out_and_err(
         gdal_grid_path
-        + " -clipsrc tmp/clip.csv -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 data/grid.vrt "
-        + outfiles[-1]
+        + f" -clipsrc {clip_csv} -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -outsize 20 20 -ot Byte -l grid -a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 data/grid.vrt {output_tif}"
     )
 
-    os.unlink("tmp/clip.csv")
-
     # We should get the same values as in "gcore/data/byte.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     cs = ds.GetRasterBand(1).Checksum()
     assert not (cs == 0 or cs == 4672), "bad checksum"
     ds = None
@@ -1288,14 +1027,9 @@ def test_gdal_grid_clipsrc(gdal_grid_path):
 
 
 @pytest.mark.require_driver("CSV")
-def test_gdal_grid_tr(gdal_grid_path):
+def test_gdal_grid_tr(gdal_grid_path, tmp_path):
 
-    #################
-    outfiles.append("tmp/grid_count_70_70.tif")
-    try:
-        os.remove(outfiles[-1])
-    except OSError:
-        pass
+    output_tif = str(tmp_path / "grid_count_70_70.tif")
 
     # Create a GDAL dataset from the values of "grid.csv".
     gdaltest.runexternal(
@@ -1303,11 +1037,11 @@ def test_gdal_grid_tr(gdal_grid_path):
         + " -txe 440720.0 441920.0 -tye 3751320.0 3750120.0 -tr 60 60 -ot Byte -l grid -a count:radius1=70.0:radius2=70.0:angle=0.0:min_points=0:nodata=0.0 "
         + " -oo X_POSSIBLE_NAMES=field_1 -oo Y_POSSIBLE_NAMES=field_2 -oo Z_POSSIBLE_NAMES=field_3"
         + " data/grid.csv "
-        + outfiles[-1]
+        + output_tif
     )
 
     # We should get the same values as in "ref_data/grid_count_70_70.tif"
-    ds = gdal.Open(outfiles[-1])
+    ds = gdal.Open(output_tif)
     ds_ref = gdal.Open("ref_data/grid_count_70_70.tif")
     assert (
         ds.GetRasterBand(1).Checksum() == ds_ref.GetRasterBand(1).Checksum()
diff --git a/autotest/utilities/test_gdal_rasterize.py b/autotest/utilities/test_gdal_rasterize.py
index daa8d8ab4704..253c3c5fdf50 100755
--- a/autotest/utilities/test_gdal_rasterize.py
+++ b/autotest/utilities/test_gdal_rasterize.py
@@ -395,12 +395,11 @@ def test_gdal_rasterize_7(gdal_rasterize_path, sql_in_file):
     else:
         sql = '"' + sql + '"'
     cmds = (
-        """tmp/test_gdal_rasterize_7.csv
-              tmp/test_gdal_rasterize_7.tif
-              -init 0 -burn 1
-              -sql %s
-              -dialect sqlite -tr 1 1 -te -1 -1 51 51"""
-        % sql
+        "tmp/test_gdal_rasterize_7.csv "
+        + "tmp/test_gdal_rasterize_7.tif "
+        + "-init 0 -burn 1 "
+        + f"-sql {sql} "
+        + "-dialect sqlite -tr 1 1 -te -1 -1 51 51"
     )
 
     gdaltest.runexternal(gdal_rasterize_path + " " + cmds)
diff --git a/autotest/utilities/test_gdal_translate.py b/autotest/utilities/test_gdal_translate.py
index 2638b364dbbe..1da926b8e587 100755
--- a/autotest/utilities/test_gdal_translate.py
+++ b/autotest/utilities/test_gdal_translate.py
@@ -969,11 +969,7 @@ def test_gdal_translate_37(gdal_translate_path):
     rat = None
     ds = None
 
-
-# Test RAT is copied round trip back to hfa
-
-
-def test_gdal_translate_38(gdal_translate_path):
+    # Test RAT is copied round trip back to hfa
 
     gdaltest.runexternal(
         gdal_translate_path
diff --git a/autotest/utilities/test_gdaladdo.py b/autotest/utilities/test_gdaladdo.py
index 477349fed5e9..77a1719d6a07 100755
--- a/autotest/utilities/test_gdaladdo.py
+++ b/autotest/utilities/test_gdaladdo.py
@@ -120,12 +120,7 @@ def test_gdaladdo_3(gdaladdo_path):
     except OSError:
         pytest.fail("no external overview.")
 
-
-###############################################################################
-# Test -clean
-
-
-def test_gdaladdo_4(gdaladdo_path):
+    # Test -clean
 
     gdaltest.runexternal(gdaladdo_path + " -clean tmp/test_gdaladdo_3.tif")
 
diff --git a/autotest/utilities/test_gdalbuildvrt.py b/autotest/utilities/test_gdalbuildvrt.py
index aba7bef83036..27b7756ee074 100755
--- a/autotest/utilities/test_gdalbuildvrt.py
+++ b/autotest/utilities/test_gdalbuildvrt.py
@@ -29,13 +29,11 @@
 # DEALINGS IN THE SOFTWARE.
 ###############################################################################
 
-import os
-
 import gdaltest
 import pytest
 import test_cli_utilities
 
-from osgeo import gdal, ogr, osr
+from osgeo import gdal, osr
 
 pytestmark = pytest.mark.skipif(
     test_cli_utilities.get_gdalbuildvrt_path() is None,
@@ -43,16 +41,53 @@
 )
 
 
-@pytest.fixture()
+@pytest.fixture(scope="module")
 def gdalbuildvrt_path():
     return test_cli_utilities.get_gdalbuildvrt_path()
 
 
+@pytest.fixture(scope="module")
+def sample_tifs(tmp_path_factory):
+
+    tmpdir = tmp_path_factory.mktemp("tmp")
+
+    drv = gdal.GetDriverByName("GTiff")
+    srs = osr.SpatialReference()
+    srs.SetWellKnownGeogCS("WGS84")
+    wkt = srs.ExportToWkt()
+
+    sample1_tif = str(tmpdir / "gdalbuildvrt1.tif")
+    with drv.Create(sample1_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
+        ds.GetRasterBand(1).Fill(0)
+
+    sample2_tif = str(tmpdir / "gdalbuildvrt2.tif")
+    with drv.Create(sample2_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([3, 0.1, 0, 49, 0, -0.1])
+        ds.GetRasterBand(1).Fill(63)
+
+    sample3_tif = str(tmpdir / "gdalbuildvrt3.tif")
+    with drv.Create(sample3_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([2, 0.1, 0, 48, 0, -0.1])
+        ds.GetRasterBand(1).Fill(127)
+
+    sample4_tif = str(tmpdir / "gdalbuildvrt4.tif")
+    with drv.Create(sample4_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([3, 0.1, 0, 48, 0, -0.1])
+        ds.GetRasterBand(1).Fill(255)
+
+    yield (sample1_tif, sample2_tif, sample3_tif, sample4_tif)
+
+
 ###############################################################################
-def gdalbuildvrt_check():
+def gdalbuildvrt_check(mosaic_vrt):
+
+    with gdal.Open(mosaic_vrt) as ds:
 
-    ds = gdal.Open("tmp/mosaic.vrt")
-    try:
         assert (
             ds.GetProjectionRef().find("WGS 84") != -1
         ), "Expected WGS 84\nGot : %s" % (ds.GetProjectionRef())
@@ -72,112 +107,72 @@ def gdalbuildvrt_check():
         assert ds.RasterCount == 1, "Wrong raster count : %d " % (ds.RasterCount)
 
         assert ds.GetRasterBand(1).Checksum() == 3508, "Wrong checksum"
-    finally:
-        del ds
 
 
 ###############################################################################
 # Simple test
 
 
-def test_gdalbuildvrt_1(gdalbuildvrt_path):
-
-    drv = gdal.GetDriverByName("GTiff")
-    srs = osr.SpatialReference()
-    srs.SetWellKnownGeogCS("WGS84")
-    wkt = srs.ExportToWkt()
-
-    ds = drv.Create("tmp/gdalbuildvrt1.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
-    ds.GetRasterBand(1).Fill(0)
-    ds = None
-
-    ds = drv.Create("tmp/gdalbuildvrt2.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([3, 0.1, 0, 49, 0, -0.1])
-    ds.GetRasterBand(1).Fill(63)
-    ds = None
+def test_gdalbuildvrt_1(gdalbuildvrt_path, tmp_path, sample_tifs):
 
-    ds = drv.Create("tmp/gdalbuildvrt3.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([2, 0.1, 0, 48, 0, -0.1])
-    ds.GetRasterBand(1).Fill(127)
-    ds = None
-
-    ds = drv.Create("tmp/gdalbuildvrt4.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([3, 0.1, 0, 48, 0, -0.1])
-    ds.GetRasterBand(1).Fill(255)
-    ds = None
+    mosaic_vrt = str(tmp_path / "mosaic.vrt")
 
     (_, err) = gdaltest.runexternal_out_and_err(
-        gdalbuildvrt_path
-        + " tmp/mosaic.vrt tmp/gdalbuildvrt1.tif tmp/gdalbuildvrt2.tif tmp/gdalbuildvrt3.tif tmp/gdalbuildvrt4.tif"
+        gdalbuildvrt_path + f" {mosaic_vrt} {' '.join(sample_tifs)}"
     )
     assert err is None or err == "", "got error/warning"
 
-    return gdalbuildvrt_check()
+    return gdalbuildvrt_check(mosaic_vrt)
 
 
 ###############################################################################
 # Test with tile index
 
 
-def test_gdalbuildvrt_2(gdalbuildvrt_path):
+def test_gdalbuildvrt_2(gdalbuildvrt_path, tmp_path, sample_tifs):
 
     if test_cli_utilities.get_gdaltindex_path() is None:
         pytest.skip()
 
-    try:
-        os.remove("tmp/tileindex.shp")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/tileindex.dbf")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/tileindex.shx")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/mosaic.vrt")
-    except OSError:
-        pass
+    tileindex_shp = str(tmp_path / "tileindex.shp")
+    mosaic_vrt = str(tmp_path / "mosaic.vrt")
 
     gdaltest.runexternal(
         test_cli_utilities.get_gdaltindex_path()
-        + " tmp/tileindex.shp tmp/gdalbuildvrt1.tif tmp/gdalbuildvrt2.tif tmp/gdalbuildvrt3.tif tmp/gdalbuildvrt4.tif"
+        + f" {tileindex_shp} {' '.join(sample_tifs)}"
     )
 
-    gdaltest.runexternal(gdalbuildvrt_path + " tmp/mosaic.vrt tmp/tileindex.shp")
+    gdaltest.runexternal(gdalbuildvrt_path + f" {mosaic_vrt} {tileindex_shp}")
 
-    return gdalbuildvrt_check()
+    return gdalbuildvrt_check(mosaic_vrt)
 
 
 ###############################################################################
 # Test with file list
 
 
-def test_gdalbuildvrt_3(gdalbuildvrt_path):
+def test_gdalbuildvrt_3(gdalbuildvrt_path, tmp_path, sample_tifs):
 
-    open("tmp/filelist.txt", "wt").write(
-        "tmp/gdalbuildvrt1.tif\ntmp/gdalbuildvrt2.tif\ntmp/gdalbuildvrt3.tif\ntmp/gdalbuildvrt4.tif"
-    )
+    mosaic_vrt = str(tmp_path / "mosaic.vrt")
+    filelist_txt = str(tmp_path / "filelist.txt")
+
+    with open(filelist_txt, "wt") as filelist:
+        filelist.write("\n".join(sample_tifs))
 
     gdaltest.runexternal(
-        gdalbuildvrt_path + " -input_file_list tmp/filelist.txt tmp/mosaic.vrt"
+        gdalbuildvrt_path + f" -input_file_list {filelist_txt} {mosaic_vrt}"
     )
 
-    return gdalbuildvrt_check()
+    return gdalbuildvrt_check(mosaic_vrt)
 
 
 ###############################################################################
 # Try adding a raster in another projection
 
 
-def test_gdalbuildvrt_4(gdalbuildvrt_path):
+def test_gdalbuildvrt_4(gdalbuildvrt_path, tmp_path, sample_tifs):
+
+    mosaic_vrt = str(tmp_path / "mosaic.vrt")
 
     drv = gdal.GetDriverByName("GTiff")
     wkt = """GEOGCS["WGS 72",
@@ -187,57 +182,59 @@ def test_gdalbuildvrt_4(gdalbuildvrt_path):
     PRIMEM["Greenwich",0],
     UNIT["degree",0.0174532925199433]]"""
 
-    ds = drv.Create("tmp/gdalbuildvrt5.tif", 10, 10, 1)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([47, 0.1, 0, 2, 0, -0.1])
-    ds = None
+    input5_tif = str(tmp_path / "gdalbuildvrt5.tif")
+
+    with drv.Create(input5_tif, 10, 10, 1) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([47, 0.1, 0, 2, 0, -0.1])
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " tmp/mosaic.vrt tmp/gdalbuildvrt1.tif tmp/gdalbuildvrt2.tif tmp/gdalbuildvrt3.tif tmp/gdalbuildvrt4.tif tmp/gdalbuildvrt5.tif"
+        gdalbuildvrt_path + f" {mosaic_vrt} {' '.join(sample_tifs)} {input5_tif}"
     )
 
-    return gdalbuildvrt_check()
+    return gdalbuildvrt_check(mosaic_vrt)
 
 
 ###############################################################################
 # Try adding a raster with different band count
 
 
-# NOTE: fails. commented out originally in 4ef886421c99a4451f8873cb6e094d45ecc86d3f, not sure why
-@pytest.mark.skip()
-def test_gdalbuildvrt_5(gdalbuildvrt_path):
+# commented out originally in 4ef886421c99a4451f8873cb6e094d45ecc86d3f, not sure why
+def test_gdalbuildvrt_5(gdalbuildvrt_path, tmp_path, sample_tifs):
+
+    mosaic_vrt = str(tmp_path / "mosaic.vrt")
 
     drv = gdal.GetDriverByName("GTiff")
     srs = osr.SpatialReference()
     srs.SetWellKnownGeogCS("WGS84")
     wkt = srs.ExportToWkt()
 
-    ds = drv.Create("tmp/gdalbuildvrt5.tif", 10, 10, 2)
-    ds.SetProjection(wkt)
-    ds.SetGeoTransform([47, 0.1, 0, 2, 0, -0.1])
-    ds = None
+    input5_tif = str(tmp_path / "gdalbuildvrt5.tif")
+
+    with drv.Create(input5_tif, 10, 10, 2) as ds:
+        ds.SetProjection(wkt)
+        ds.SetGeoTransform([47, 0.1, 0, 2, 0, -0.1])
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " tmp/mosaic.vrt tmp/gdalbuildvrt1.tif tmp/gdalbuildvrt2.tif tmp/gdalbuildvrt3.tif tmp/gdalbuildvrt4.tif tmp/gdalbuildvrt5.tif"
+        gdalbuildvrt_path + f" {mosaic_vrt} {' '.join(sample_tifs)} {input5_tif}"
     )
 
-    return gdalbuildvrt_check()
+    return gdalbuildvrt_check(mosaic_vrt)
 
 
 ###############################################################################
 # Test -separate option
 
 
-def test_gdalbuildvrt_6(gdalbuildvrt_path):
+def test_gdalbuildvrt_6(gdalbuildvrt_path, tmp_path, sample_tifs):
+
+    output_vrt = str(tmp_path / "stacked.vrt")
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " -separate tmp/stacked.vrt tmp/gdalbuildvrt1.tif tmp/gdalbuildvrt2.tif tmp/gdalbuildvrt3.tif tmp/gdalbuildvrt4.tif"
+        gdalbuildvrt_path + f" -separate {output_vrt} {' '.join(sample_tifs)}"
     )
 
-    ds = gdal.Open("tmp/stacked.vrt")
+    ds = gdal.Open(output_vrt)
     assert ds.GetProjectionRef().find("WGS 84") != -1, "Expected WGS 84\nGot : %s" % (
         ds.GetProjectionRef()
     )
@@ -263,86 +260,86 @@ def test_gdalbuildvrt_6(gdalbuildvrt_path):
 # Test source rasters with nodata
 
 
-def test_gdalbuildvrt_7(gdalbuildvrt_path):
+def test_gdalbuildvrt_7(gdalbuildvrt_path, tmp_path):
 
-    out_ds = gdal.GetDriverByName("GTiff").Create(
-        "tmp/vrtnull1.tif", 20, 10, 3, gdal.GDT_UInt16
-    )
-    out_ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
-    srs = osr.SpatialReference()
-    srs.SetFromUserInput("EPSG:4326")
-    out_ds.SetProjection(srs.ExportToWkt())
-    out_ds.GetRasterBand(1).SetRasterColorInterpretation(gdal.GCI_RedBand)
-    out_ds.GetRasterBand(2).SetRasterColorInterpretation(gdal.GCI_GreenBand)
-    out_ds.GetRasterBand(3).SetRasterColorInterpretation(gdal.GCI_BlueBand)
-    out_ds.GetRasterBand(1).SetNoDataValue(256)
-
-    try:
-        ff = "\xff".encode("latin1")
-    except UnicodeDecodeError:
-        ff = "\xff"
-
-    out_ds.GetRasterBand(1).WriteRaster(
-        0, 0, 10, 10, ff, buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
-    )
-    out_ds.GetRasterBand(2).WriteRaster(
-        0, 0, 10, 10, "\x00", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
-    )
-    out_ds.GetRasterBand(3).WriteRaster(
-        0, 0, 10, 10, "\x00", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
-    )
-    out_ds = None
+    input1_tif = str(tmp_path / "vrtnull1.tif")
+    input2_tif = str(tmp_path / "vrtnull2.tif")
+    output_vrt = str(tmp_path / "gdalbuildvrt7.vrt")
 
-    out_ds = gdal.GetDriverByName("GTiff").Create(
-        "tmp/vrtnull2.tif", 20, 10, 3, gdal.GDT_UInt16
-    )
-    out_ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
-    srs = osr.SpatialReference()
-    srs.SetFromUserInput("EPSG:4326")
-    out_ds.SetProjection(srs.ExportToWkt())
-    out_ds.GetRasterBand(1).SetRasterColorInterpretation(gdal.GCI_RedBand)
-    out_ds.GetRasterBand(2).SetRasterColorInterpretation(gdal.GCI_GreenBand)
-    out_ds.GetRasterBand(3).SetRasterColorInterpretation(gdal.GCI_BlueBand)
-    out_ds.GetRasterBand(1).SetNoDataValue(256)
-
-    out_ds.GetRasterBand(1).WriteRaster(
-        10, 0, 10, 10, "\x00", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
-    )
-    out_ds.GetRasterBand(2).WriteRaster(
-        10, 0, 10, 10, ff, buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
-    )
-    out_ds.GetRasterBand(3).WriteRaster(
-        10, 0, 10, 10, "\x00", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
-    )
-    out_ds = None
+    with gdal.GetDriverByName("GTiff").Create(
+        input1_tif, 20, 10, 3, gdal.GDT_UInt16
+    ) as out_ds:
+        out_ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
+        srs = osr.SpatialReference()
+        srs.SetFromUserInput("EPSG:4326")
+        out_ds.SetProjection(srs.ExportToWkt())
+        out_ds.GetRasterBand(1).SetRasterColorInterpretation(gdal.GCI_RedBand)
+        out_ds.GetRasterBand(2).SetRasterColorInterpretation(gdal.GCI_GreenBand)
+        out_ds.GetRasterBand(3).SetRasterColorInterpretation(gdal.GCI_BlueBand)
+        out_ds.GetRasterBand(1).SetNoDataValue(256)
 
-    gdaltest.runexternal(
-        gdalbuildvrt_path + " tmp/gdalbuildvrt7.vrt tmp/vrtnull1.tif tmp/vrtnull2.tif"
-    )
+        try:
+            ff = "\xff".encode("latin1")
+        except UnicodeDecodeError:
+            ff = "\xff"
+
+        out_ds.GetRasterBand(1).WriteRaster(
+            0, 0, 10, 10, ff, buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
+        )
+        out_ds.GetRasterBand(2).WriteRaster(
+            0, 0, 10, 10, "\x00", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
+        )
+        out_ds.GetRasterBand(3).WriteRaster(
+            0, 0, 10, 10, "\x00", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
+        )
 
-    ds = gdal.Open("tmp/gdalbuildvrt7.vrt")
+    with gdal.GetDriverByName("GTiff").Create(
+        input2_tif, 20, 10, 3, gdal.GDT_UInt16
+    ) as out_ds:
+        out_ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
+        srs = osr.SpatialReference()
+        srs.SetFromUserInput("EPSG:4326")
+        out_ds.SetProjection(srs.ExportToWkt())
+        out_ds.GetRasterBand(1).SetRasterColorInterpretation(gdal.GCI_RedBand)
+        out_ds.GetRasterBand(2).SetRasterColorInterpretation(gdal.GCI_GreenBand)
+        out_ds.GetRasterBand(3).SetRasterColorInterpretation(gdal.GCI_BlueBand)
+        out_ds.GetRasterBand(1).SetNoDataValue(256)
+
+        out_ds.GetRasterBand(1).WriteRaster(
+            10, 0, 10, 10, "\x00", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
+        )
+        out_ds.GetRasterBand(2).WriteRaster(
+            10, 0, 10, 10, ff, buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
+        )
+        out_ds.GetRasterBand(3).WriteRaster(
+            10, 0, 10, 10, "\x00", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
+        )
 
-    assert ds.GetRasterBand(1).Checksum() == 1217, "Wrong checksum"
+    gdaltest.runexternal(gdalbuildvrt_path + f" {output_vrt} {input1_tif} {input2_tif}")
 
-    assert ds.GetRasterBand(2).Checksum() == 1218, "Wrong checksum"
+    with gdal.Open(output_vrt) as ds:
 
-    assert ds.GetRasterBand(3).Checksum() == 0, "Wrong checksum"
+        assert ds.GetRasterBand(1).Checksum() == 1217, "Wrong checksum"
 
-    ds = None
+        assert ds.GetRasterBand(2).Checksum() == 1218, "Wrong checksum"
+
+        assert ds.GetRasterBand(3).Checksum() == 0, "Wrong checksum"
 
 
 ###############################################################################
 # Test -tr option
 
 
-def test_gdalbuildvrt_8(gdalbuildvrt_path):
+def test_gdalbuildvrt_8(gdalbuildvrt_path, tmp_path, sample_tifs):
+
+    mosaic_vrt = str(tmp_path / "mosiac.vrt")
+    mosaic2_vrt = str(tmp_path / "mosaic2.vrt")
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " -tr 0.05 0.05 tmp/mosaic2.vrt tmp/gdalbuildvrt1.tif tmp/gdalbuildvrt2.tif tmp/gdalbuildvrt3.tif tmp/gdalbuildvrt4.tif"
+        gdalbuildvrt_path + f" -tr 0.05 0.05 {mosaic2_vrt} {' '.join(sample_tifs)}"
     )
 
-    ds = gdal.Open("tmp/mosaic2.vrt")
+    ds = gdal.Open(mosaic2_vrt)
 
     gt = ds.GetGeoTransform()
     expected_gt = [2, 0.05, 0, 49, 0, -0.05]
@@ -356,25 +353,25 @@ def test_gdalbuildvrt_8(gdalbuildvrt_path):
         ds.RasterXSize == 40 and ds.RasterYSize == 40
     ), "Wrong raster dimensions : %d x %d" % (ds.RasterXSize, ds.RasterYSize)
 
-    gdaltest.runexternal(
-        gdalbuildvrt_path + " -tr 0.1 0.1 tmp/mosaic.vrt tmp/mosaic2.vrt"
-    )
+    gdaltest.runexternal(gdalbuildvrt_path + f" -tr 0.1 0.1 {mosaic_vrt} {mosaic2_vrt}")
 
-    return gdalbuildvrt_check()
+    return gdalbuildvrt_check(mosaic_vrt)
 
 
 ###############################################################################
 # Test -te option
 
 
-def test_gdalbuildvrt_9(gdalbuildvrt_path):
+def test_gdalbuildvrt_9(gdalbuildvrt_path, tmp_path, sample_tifs):
+
+    mosaic_vrt = str(tmp_path / "mosiac.vrt")
+    mosaic2_vrt = str(tmp_path / "mosaic2.vrt")
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " -te 1 46 5 50 tmp/mosaic2.vrt tmp/gdalbuildvrt1.tif tmp/gdalbuildvrt2.tif tmp/gdalbuildvrt3.tif tmp/gdalbuildvrt4.tif"
+        gdalbuildvrt_path + f" -te 1 46 5 50 {mosaic2_vrt} {' '.join(sample_tifs)}"
     )
 
-    ds = gdal.Open("tmp/mosaic2.vrt")
+    ds = gdal.Open(mosaic2_vrt)
 
     gt = ds.GetGeoTransform()
     expected_gt = [1, 0.1, 0, 50, 0, -0.1]
@@ -389,60 +386,61 @@ def test_gdalbuildvrt_9(gdalbuildvrt_path):
     ), "Wrong raster dimensions : %d x %d" % (ds.RasterXSize, ds.RasterYSize)
 
     gdaltest.runexternal(
-        gdalbuildvrt_path + " -te 2 47 4 49 tmp/mosaic.vrt tmp/mosaic2.vrt"
+        gdalbuildvrt_path + f" -te 2 47 4 49 {mosaic_vrt} {mosaic2_vrt}"
     )
 
-    return gdalbuildvrt_check()
+    return gdalbuildvrt_check(mosaic_vrt)
 
 
 ###############################################################################
 # Test explicit nodata setting (#3254)
 
 
-def test_gdalbuildvrt_10(gdalbuildvrt_path):
+def test_gdalbuildvrt_10(gdalbuildvrt_path, tmp_path):
 
-    out_ds = gdal.GetDriverByName("GTiff").Create(
-        "tmp/test_gdalbuildvrt_10_1.tif",
+    input1_tif = str(tmp_path / "test_gdalbuildvrt_10_1.tif")
+    input2_tif = str(tmp_path / "test_gdalbuildvrt_10_2.tif")
+    output_vrt = str(tmp_path / "gdalbuildvrt10.vrt")
+
+    with gdal.GetDriverByName("GTiff").Create(
+        input1_tif,
         10,
         10,
         1,
         gdal.GDT_Byte,
         options=["NBITS=1", "PHOTOMETRIC=MINISWHITE"],
-    )
-    out_ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
-    srs = osr.SpatialReference()
-    srs.SetFromUserInput("EPSG:4326")
-    out_ds.SetProjection(srs.ExportToWkt())
-
-    out_ds.GetRasterBand(1).WriteRaster(
-        1, 1, 3, 3, "\x01", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
-    )
-    out_ds = None
+    ) as out_ds:
+        out_ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
+        srs = osr.SpatialReference()
+        srs.SetFromUserInput("EPSG:4326")
+        out_ds.SetProjection(srs.ExportToWkt())
+
+        out_ds.GetRasterBand(1).WriteRaster(
+            1, 1, 3, 3, "\x01", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
+        )
 
-    out_ds = gdal.GetDriverByName("GTiff").Create(
-        "tmp/test_gdalbuildvrt_10_2.tif",
+    with gdal.GetDriverByName("GTiff").Create(
+        input2_tif,
         10,
         10,
         1,
         gdal.GDT_Byte,
         options=["NBITS=1", "PHOTOMETRIC=MINISWHITE"],
-    )
-    out_ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
-    srs = osr.SpatialReference()
-    srs.SetFromUserInput("EPSG:4326")
-    out_ds.SetProjection(srs.ExportToWkt())
-
-    out_ds.GetRasterBand(1).WriteRaster(
-        6, 6, 3, 3, "\x01", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
-    )
-    out_ds = None
+    ) as out_ds:
+        out_ds.SetGeoTransform([2, 0.1, 0, 49, 0, -0.1])
+        srs = osr.SpatialReference()
+        srs.SetFromUserInput("EPSG:4326")
+        out_ds.SetProjection(srs.ExportToWkt())
+
+        out_ds.GetRasterBand(1).WriteRaster(
+            6, 6, 3, 3, "\x01", buf_type=gdal.GDT_Byte, buf_xsize=1, buf_ysize=1
+        )
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " -srcnodata 0 tmp/gdalbuildvrt10.vrt tmp/test_gdalbuildvrt_10_1.tif tmp/test_gdalbuildvrt_10_2.tif"
+        gdalbuildvrt_path + f" -srcnodata 0 {output_vrt} {input1_tif} {input2_tif}"
     )
 
-    ds = gdal.Open("tmp/gdalbuildvrt10.vrt")
+    ds = gdal.Open(output_vrt)
 
     assert ds.GetRasterBand(1).Checksum() == 18, "Wrong checksum"
 
@@ -453,56 +451,52 @@ def test_gdalbuildvrt_10(gdalbuildvrt_path):
 # Test that we can stack ungeoreference single band images with -separate (#3432)
 
 
-def test_gdalbuildvrt_11(gdalbuildvrt_path):
+def test_gdalbuildvrt_11(gdalbuildvrt_path, tmp_path):
 
-    out_ds = gdal.GetDriverByName("GTiff").Create(
-        "tmp/test_gdalbuildvrt_11_1.tif", 10, 10, 1
-    )
-    out_ds.GetRasterBand(1).Fill(255)
-    cs1 = out_ds.GetRasterBand(1).Checksum()
-    out_ds = None
+    output_vrt = str(tmp_path / "gdalbuildvrt11.vrt")
 
-    out_ds = gdal.GetDriverByName("GTiff").Create(
-        "tmp/test_gdalbuildvrt_11_2.tif", 10, 10, 1
-    )
-    out_ds.GetRasterBand(1).Fill(127)
-    cs2 = out_ds.GetRasterBand(1).Checksum()
-    out_ds = None
+    input_tif1 = str(tmp_path / "test_gdalbuildvrt_11_1.tif")
+    input_tif2 = str(tmp_path / "test_gdalbuildvrt_11_2.tif")
+
+    with gdal.GetDriverByName("GTiff").Create(input_tif1, 10, 10, 1) as out_ds:
+        out_ds.GetRasterBand(1).Fill(255)
+        cs1 = out_ds.GetRasterBand(1).Checksum()
+
+    with gdal.GetDriverByName("GTiff").Create(input_tif2, 10, 10, 1) as out_ds:
+        out_ds.GetRasterBand(1).Fill(127)
+        cs2 = out_ds.GetRasterBand(1).Checksum()
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " -separate tmp/gdalbuildvrt11.vrt tmp/test_gdalbuildvrt_11_1.tif tmp/test_gdalbuildvrt_11_2.tif"
+        gdalbuildvrt_path + f" -separate {output_vrt} {input_tif1} {input_tif2}"
     )
 
-    ds = gdal.Open("tmp/gdalbuildvrt11.vrt")
+    with gdal.Open(output_vrt) as ds:
 
-    assert ds.GetRasterBand(1).Checksum() == cs1, "Wrong checksum"
+        assert ds.GetRasterBand(1).Checksum() == cs1, "Wrong checksum"
 
-    assert ds.GetRasterBand(2).Checksum() == cs2, "Wrong checksum"
-
-    ds = None
+        assert ds.GetRasterBand(2).Checksum() == cs2, "Wrong checksum"
 
 
 ###############################################################################
 # Test -tap option
 
 
-def test_gdalbuildvrt_12(gdalbuildvrt_path):
+def test_gdalbuildvrt_12(gdalbuildvrt_path, tmp_path):
+
+    output_vrt = str(tmp_path / "gdalbuildvrt12.vrt")
 
     (_, err) = gdaltest.runexternal_out_and_err(
-        gdalbuildvrt_path + " -tap tmp/gdalbuildvrt12.vrt ../gcore/data/byte.tif",
+        gdalbuildvrt_path + f" -tap {output_vrt} ../gcore/data/byte.tif",
         check_memleak=False,
     )
-    assert (
-        err.find("-tap option cannot be used without using -tr") != -1
-    ), "expected error"
+
+    assert "-tap option cannot be used without using -tr" in err, "expected error"
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " -tr 100 50 -tap tmp/gdalbuildvrt12.vrt ../gcore/data/byte.tif"
+        gdalbuildvrt_path + f" -tr 100 50 -tap {output_vrt} ../gcore/data/byte.tif"
     )
 
-    ds = gdal.Open("tmp/gdalbuildvrt12.vrt")
+    ds = gdal.Open(output_vrt)
 
     gt = ds.GetGeoTransform()
     expected_gt = [440700.0, 100.0, 0.0, 3751350.0, 0.0, -50.0]
@@ -521,42 +515,45 @@ def test_gdalbuildvrt_12(gdalbuildvrt_path):
 # Test -a_srs
 
 
-def test_gdalbuildvrt_13(gdalbuildvrt_path):
+def test_gdalbuildvrt_13(gdalbuildvrt_path, tmp_path):
+
+    output_vrt = str(tmp_path / "gdalbuildvrt13.vrt")
 
     gdaltest.runexternal(
-        gdalbuildvrt_path
-        + " tmp/gdalbuildvrt13.vrt ../gcore/data/byte.tif -a_srs EPSG:4326"
+        gdalbuildvrt_path + f" {output_vrt} ../gcore/data/byte.tif -a_srs EPSG:4326"
     )
 
-    ds = gdal.Open("tmp/gdalbuildvrt13.vrt")
-    assert ds.GetProjectionRef().find("4326") != -1
-    ds = None
+    with gdal.Open(output_vrt) as ds:
+        assert "4326" in ds.GetProjectionRef()
 
 
 ###############################################################################
 # Test -r
 
 
-def test_gdalbuildvrt_14(gdalbuildvrt_path):
+def test_gdalbuildvrt_14(gdalbuildvrt_path, tmp_path):
     if test_cli_utilities.get_gdal_translate_path() is None:
         pytest.skip()
 
+    buildvrt_output_vrt = str(tmp_path / "test_gdalbuildvrt_14.vrt")
+
     gdaltest.runexternal(
         gdalbuildvrt_path
-        + " tmp/test_gdalbuildvrt_14.vrt ../gcore/data/byte.tif -r cubic -tr 30 30"
+        + f" {buildvrt_output_vrt} ../gcore/data/byte.tif -r cubic -tr 30 30"
     )
 
+    translate_output_vrt = str(tmp_path / "test_gdalbuildvrt_14_ref.vrt")
+
     gdaltest.runexternal(
         test_cli_utilities.get_gdal_translate_path()
-        + " -of VRT ../gcore/data/byte.tif tmp/test_gdalbuildvrt_14_ref.vrt -r cubic -outsize 40 40"
+        + f" -of VRT ../gcore/data/byte.tif {translate_output_vrt} -r cubic -outsize 40 40"
     )
 
-    ds = gdal.Open("tmp/test_gdalbuildvrt_14.vrt")
-    ds_ref = gdal.Open("tmp/test_gdalbuildvrt_14_ref.vrt")
-    cs = ds.GetRasterBand(1).Checksum()
-    cs_ref = ds_ref.GetRasterBand(1).Checksum()
-    ds = None
-    ds_ref = None
+    with gdal.Open(buildvrt_output_vrt) as ds:
+        cs = ds.GetRasterBand(1).Checksum()
+
+    with gdal.Open(translate_output_vrt) as ds_ref:
+        cs_ref = ds_ref.GetRasterBand(1).Checksum()
 
     assert cs == cs_ref
 
@@ -565,17 +562,16 @@ def test_gdalbuildvrt_14(gdalbuildvrt_path):
 # Test -b
 
 
-def test_gdalbuildvrt_15(gdalbuildvrt_path):
+def test_gdalbuildvrt_15(gdalbuildvrt_path, tmp_path):
+
+    output_vrt = str(tmp_path / "test_gdalbuildvrt_15.vrt")
 
     gdaltest.runexternal(
-        gdalbuildvrt_path + " tmp/test_gdalbuildvrt_15.vrt ../gcore/data/byte.tif -b 1"
+        gdalbuildvrt_path + f" {output_vrt}  ../gcore/data/byte.tif -b 1"
     )
 
-    ds = gdal.Open("tmp/test_gdalbuildvrt_15.vrt")
-    cs = ds.GetRasterBand(1).Checksum()
-    ds = None
-
-    assert cs == 4672
+    with gdal.Open(output_vrt) as ds:
+        assert ds.GetRasterBand(1).Checksum() == 4672
 
 
 ###############################################################################
@@ -594,45 +590,3 @@ def test_gdalbuildvrt_16(gdalbuildvrt_path):
     else:
         # We don't get the error code on Travis mingw
         assert "ERROR" in err, out
-
-
-###############################################################################
-# Cleanup
-
-
-def test_gdalbuildvrt_cleanup():
-
-    if gdalbuildvrt_path is None:
-        pytest.skip()
-
-    ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource("tmp/tileindex.shp")
-
-    gdal.GetDriverByName("VRT").Delete("tmp/mosaic.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/mosaic2.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/stacked.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/gdalbuildvrt7.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/gdalbuildvrt10.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/gdalbuildvrt11.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/gdalbuildvrt12.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/gdalbuildvrt13.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/test_gdalbuildvrt_14.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/test_gdalbuildvrt_14_ref.vrt")
-    gdal.GetDriverByName("VRT").Delete("tmp/test_gdalbuildvrt_15.vrt")
-
-    drv = gdal.GetDriverByName("GTiff")
-
-    drv.Delete("tmp/gdalbuildvrt1.tif")
-    drv.Delete("tmp/gdalbuildvrt2.tif")
-    drv.Delete("tmp/gdalbuildvrt3.tif")
-    drv.Delete("tmp/gdalbuildvrt4.tif")
-    drv.Delete("tmp/gdalbuildvrt5.tif")
-    drv.Delete("tmp/vrtnull1.tif")
-    drv.Delete("tmp/vrtnull2.tif")
-    drv.Delete("tmp/test_gdalbuildvrt_10_1.tif")
-    drv.Delete("tmp/test_gdalbuildvrt_10_2.tif")
-    drv.Delete("tmp/test_gdalbuildvrt_11_1.tif")
-    drv.Delete("tmp/test_gdalbuildvrt_11_2.tif")
-    try:
-        os.remove("tmp/filelist.txt")
-    except OSError:
-        pass
diff --git a/autotest/utilities/test_gdaldem.py b/autotest/utilities/test_gdaldem.py
index c9c38a8417dc..a1d1d1326089 100755
--- a/autotest/utilities/test_gdaldem.py
+++ b/autotest/utilities/test_gdaldem.py
@@ -100,11 +100,11 @@ def test_gdaldem_hillshade_compressed_tiled_output(gdaldem_path):
 
     ds = None
 
-    stat_uncompressed = os.stat("tmp/n43_hillshade.tif")
     stat_compressed = os.stat("tmp/n43_hillshade_compressed_tiled.tif")
+
     assert (
-        stat_uncompressed.st_size >= stat_compressed.st_size
-    ), "failure: compressed size greater than uncompressed one"
+        stat_compressed.st_size <= 15027
+    ), "compressed size greater than uncompressed one"
 
 
 ###############################################################################
@@ -299,14 +299,22 @@ def test_gdaldem_aspect(gdaldem_path):
 # Test gdaldem color relief
 
 
-def test_gdaldem_color_relief(gdaldem_path):
+@pytest.fixture()
+def n43_colorrelief_tif(gdaldem_path, tmp_path):
+    n43_colorrelief_tif = str(tmp_path / "n43_colorrelief.tif")
 
     gdaltest.runexternal(
         gdaldem_path
-        + " color-relief ../gdrivers/data/n43.tif data/color_file.txt tmp/n43_colorrelief.tif"
+        + f" color-relief ../gdrivers/data/n43.tif data/color_file.txt {n43_colorrelief_tif}"
     )
+
+    yield n43_colorrelief_tif
+
+
+def test_gdaldem_color_relief(gdaldem_path, n43_colorrelief_tif):
+
     src_ds = gdal.Open("../gdrivers/data/n43.tif")
-    ds = gdal.Open("tmp/n43_colorrelief.tif")
+    ds = gdal.Open(n43_colorrelief_tif)
     assert ds is not None
 
     assert ds.GetRasterBand(1).Checksum() == 55009, "Bad checksum"
@@ -363,7 +371,7 @@ def test_gdaldem_color_relief_cpt(gdaldem_path):
 # Test gdaldem color relief to VRT
 
 
-def test_gdaldem_color_relief_vrt(gdaldem_path):
+def test_gdaldem_color_relief_vrt(gdaldem_path, n43_colorrelief_tif):
 
     gdaltest.runexternal(
         gdaldem_path
@@ -373,7 +381,7 @@ def test_gdaldem_color_relief_vrt(gdaldem_path):
     ds = gdal.Open("tmp/n43_colorrelief.vrt")
     assert ds is not None
 
-    ds_ref = gdal.Open("tmp/n43_colorrelief.tif")
+    ds_ref = gdal.Open(n43_colorrelief_tif)
     assert gdaltest.compare_ds(ds, ds_ref, verbose=0) <= 1, "Bad checksum"
     ds_ref = None
 
@@ -393,15 +401,21 @@ def test_gdaldem_color_relief_vrt(gdaldem_path):
 # Test gdaldem color relief from a Float32 dataset
 
 
-def test_gdaldem_color_relief_from_float32(gdaldem_path):
+@pytest.fixture()
+def n43_float32_tif(tmp_path):
 
-    gdal.Translate(
-        "tmp/n43_float32.tif", "../gdrivers/data/n43.tif", options="-ot Float32"
-    )
+    n43_float32_path = str(tmp_path / "n43_float32.tif")
+
+    gdal.Translate(n43_float32_path, "../gdrivers/data/n43.tif", options="-ot Float32")
+
+    yield n43_float32_path
+
+
+def test_gdaldem_color_relief_from_float32(gdaldem_path, n43_float32_tif):
 
     gdaltest.runexternal(
         gdaldem_path
-        + " color-relief tmp/n43_float32.tif data/color_file.txt tmp/n43_colorrelief_from_float32.tif"
+        + f" color-relief {n43_float32_tif} data/color_file.txt tmp/n43_colorrelief_from_float32.tif"
     )
     ds = gdal.Open("tmp/n43_colorrelief_from_float32.tif")
     assert ds is not None
@@ -443,11 +457,11 @@ def test_gdaldem_color_relief_png(gdaldem_path):
 
 
 @pytest.mark.require_driver("PNG")
-def test_gdaldem_color_relief_from_float32_to_png(gdaldem_path):
+def test_gdaldem_color_relief_from_float32_to_png(gdaldem_path, n43_float32_tif):
 
     gdaltest.runexternal(
         gdaldem_path
-        + " color-relief -of PNG tmp/n43_float32.tif data/color_file.txt tmp/n43_colorrelief_from_float32.png"
+        + f" color-relief -of PNG {n43_float32_tif} data/color_file.txt tmp/n43_colorrelief_from_float32.png"
     )
     ds = gdal.Open("tmp/n43_colorrelief_from_float32.png")
     assert ds is not None
diff --git a/autotest/utilities/test_gdaltindex.py b/autotest/utilities/test_gdaltindex.py
index f708de634212..16690a23e641 100755
--- a/autotest/utilities/test_gdaltindex.py
+++ b/autotest/utilities/test_gdaltindex.py
@@ -37,9 +37,13 @@
 
 from osgeo import gdal, ogr, osr
 
-pytestmark = pytest.mark.skipif(
-    test_cli_utilities.get_gdaltindex_path() is None, reason="gdaltindex not available"
-)
+pytestmark = [
+    pytest.mark.skipif(
+        test_cli_utilities.get_gdaltindex_path() is None,
+        reason="gdaltindex not available",
+    ),
+    pytest.mark.random_order(disabled=True),
+]
 
 
 @pytest.fixture()
diff --git a/autotest/utilities/test_gdalwarp.py b/autotest/utilities/test_gdalwarp.py
index 0efac183e7d5..71f6c435608d 100755
--- a/autotest/utilities/test_gdalwarp.py
+++ b/autotest/utilities/test_gdalwarp.py
@@ -48,6 +48,107 @@ def gdalwarp_path():
     return test_cli_utilities.get_gdalwarp_path()
 
 
+@pytest.fixture(scope="module", autouse=True)
+def setup_and_cleanup():
+
+    yield
+
+    # We don't clean up when run in debug mode.
+    if gdal.GetConfigOption("CPL_DEBUG", "OFF") == "ON":
+        return
+
+    for i in range(37):
+        try:
+            os.remove("tmp/testgdalwarp" + str(i + 1) + ".tif")
+        except OSError:
+            pass
+        try:
+            os.remove("tmp/testgdalwarp" + str(i + 1) + ".vrt")
+        except OSError:
+            pass
+        try:
+            os.remove("tmp/testgdalwarp" + str(i + 1) + ".tif.aux.xml")
+        except OSError:
+            pass
+    try:
+        os.remove("tmp/testgdalwarp24src.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/testgdalwarp24dst.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/testgdalwarp30_1.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/testgdalwarp30_2.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/testgdalwarp30_3.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/testgdalwarp33_mask.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/testgdalwarp37.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/testgdalwarp38.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/test_gdalwarp_39.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/test_gdalwarp_40_src.tif")
+        os.remove("tmp/test_gdalwarp_40.tif")
+        os.remove("tmp/test_gdalwarp_40.vrt")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/test_gdalwarp_41_src.tif")
+        os.remove("tmp/test_gdalwarp_41.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/small_world_left.tif")
+        os.remove("tmp/small_world_right.tif")
+        os.remove("tmp/test_gdalwarp_42.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/small_world.tif")
+        os.remove("tmp/test_gdalwarp_43.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/test_gdalwarp_44.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/test_gdalwarp_45.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/test_gdalwarp_46.tif")
+    except OSError:
+        pass
+    try:
+        os.remove("tmp/cutline_4326.shp")
+        os.remove("tmp/cutline_4326.shx")
+        os.remove("tmp/cutline_4326.dbf")
+        os.remove("tmp/cutline_4326.prj")
+    except OSError:
+        pass
+
+
 ###############################################################################
 # Simple test
 
@@ -128,18 +229,28 @@ def test_gdalwarp_4(gdalwarp_path):
 # Test warping from GCPs without any explicit option
 
 
-def test_gdalwarp_5(gdalwarp_path):
+@pytest.fixture(scope="module")
+def testgdalwarp_gcp_tif(tmp_path_factory):
+
+    testgdalwarp_gcp_tif_fname = str(
+        tmp_path_factory.mktemp("tmp") / "testgdalwarp_gcp.tif"
+    )
 
     if test_cli_utilities.get_gdal_translate_path() is None:
         pytest.skip("gdal_translate missing")
 
     gdaltest.runexternal(
         test_cli_utilities.get_gdal_translate_path()
-        + " -a_srs EPSG:26711 -gcp 0 0  440720.000 3751320.000 -gcp 20 0 441920.000 3751320.000 -gcp 20 20 441920.000 3750120.000 0 -gcp 0 20 440720.000 3750120.000 ../gcore/data/byte.tif tmp/testgdalwarp_gcp.tif"
+        + f" -a_srs EPSG:26711 -gcp 0 0  440720.000 3751320.000 -gcp 20 0 441920.000 3751320.000 -gcp 20 20 441920.000 3750120.000 0 -gcp 0 20 440720.000 3750120.000 ../gcore/data/byte.tif {testgdalwarp_gcp_tif_fname}"
     )
 
+    yield testgdalwarp_gcp_tif_fname
+
+
+def test_gdalwarp_5(gdalwarp_path, testgdalwarp_gcp_tif):
+
     gdaltest.runexternal(
-        gdalwarp_path + " tmp/testgdalwarp_gcp.tif tmp/testgdalwarp5.tif"
+        gdalwarp_path + f" {testgdalwarp_gcp_tif} tmp/testgdalwarp5.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp5.tif")
@@ -160,10 +271,10 @@ def test_gdalwarp_5(gdalwarp_path):
 # Test warping from GCPs with -tps
 
 
-def test_gdalwarp_6(gdalwarp_path):
+def test_gdalwarp_6(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -tps tmp/testgdalwarp_gcp.tif tmp/testgdalwarp6.tif"
+        gdalwarp_path + f" -tps {testgdalwarp_gcp_tif} tmp/testgdalwarp6.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp6.tif")
@@ -184,10 +295,10 @@ def test_gdalwarp_6(gdalwarp_path):
 # Test -tr
 
 
-def test_gdalwarp_7(gdalwarp_path):
+def test_gdalwarp_7(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -tr 120 120 tmp/testgdalwarp_gcp.tif tmp/testgdalwarp7.tif"
+        gdalwarp_path + f" -tr 120 120 {testgdalwarp_gcp_tif} tmp/testgdalwarp7.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp7.tif")
@@ -203,10 +314,10 @@ def test_gdalwarp_7(gdalwarp_path):
 # Test -ts
 
 
-def test_gdalwarp_8(gdalwarp_path):
+def test_gdalwarp_8(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -ts 10 10 tmp/testgdalwarp_gcp.tif tmp/testgdalwarp8.tif"
+        gdalwarp_path + f" -ts 10 10 {testgdalwarp_gcp_tif} tmp/testgdalwarp8.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp8.tif")
@@ -222,11 +333,11 @@ def test_gdalwarp_8(gdalwarp_path):
 # Test -te
 
 
-def test_gdalwarp_9(gdalwarp_path):
+def test_gdalwarp_9(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
         gdalwarp_path
-        + " -te 440720.000 3750120.000 441920.000 3751320.000 tmp/testgdalwarp_gcp.tif tmp/testgdalwarp9.tif"
+        + f" -te 440720.000 3750120.000 441920.000 3751320.000 {testgdalwarp_gcp_tif} tmp/testgdalwarp9.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp9.tif")
@@ -245,10 +356,10 @@ def test_gdalwarp_9(gdalwarp_path):
 # Test -rn
 
 
-def test_gdalwarp_10(gdalwarp_path):
+def test_gdalwarp_10(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -ts 40 40 -rn tmp/testgdalwarp_gcp.tif tmp/testgdalwarp10.tif"
+        gdalwarp_path + f" -ts 40 40 -rn {testgdalwarp_gcp_tif} tmp/testgdalwarp10.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp10.tif")
@@ -263,10 +374,10 @@ def test_gdalwarp_10(gdalwarp_path):
 # Test -rb
 
 
-def test_gdalwarp_11(gdalwarp_path):
+def test_gdalwarp_11(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -ts 40 40 -rb tmp/testgdalwarp_gcp.tif tmp/testgdalwarp11.tif"
+        gdalwarp_path + f" -ts 40 40 -rb {testgdalwarp_gcp_tif} tmp/testgdalwarp11.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp11.tif")
@@ -287,10 +398,10 @@ def test_gdalwarp_11(gdalwarp_path):
 # Test -rc
 
 
-def test_gdalwarp_12(gdalwarp_path):
+def test_gdalwarp_12(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -ts 40 40 -rc tmp/testgdalwarp_gcp.tif tmp/testgdalwarp12.tif"
+        gdalwarp_path + f" -ts 40 40 -rc {testgdalwarp_gcp_tif} tmp/testgdalwarp12.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp12.tif")
@@ -312,11 +423,10 @@ def test_gdalwarp_12(gdalwarp_path):
 # Test -rcs
 
 
-def test_gdalwarp_13(gdalwarp_path):
+def test_gdalwarp_13(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path
-        + " -ts 40 40 -rcs tmp/testgdalwarp_gcp.tif tmp/testgdalwarp13.tif"
+        gdalwarp_path + f" -ts 40 40 -rcs {testgdalwarp_gcp_tif} tmp/testgdalwarp13.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp13.tif")
@@ -335,11 +445,11 @@ def test_gdalwarp_13(gdalwarp_path):
 # Test -r lanczos
 
 
-def test_gdalwarp_14(gdalwarp_path):
+def test_gdalwarp_14(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
         gdalwarp_path
-        + " -ts 40 40 -r lanczos tmp/testgdalwarp_gcp.tif tmp/testgdalwarp14.tif"
+        + f" -ts 40 40 -r lanczos {testgdalwarp_gcp_tif} tmp/testgdalwarp14.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp14.tif")
@@ -358,10 +468,10 @@ def test_gdalwarp_14(gdalwarp_path):
 # Test -of VRT which is a special case
 
 
-def test_gdalwarp_16(gdalwarp_path):
+def test_gdalwarp_16(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -of VRT tmp/testgdalwarp_gcp.tif tmp/testgdalwarp16.vrt"
+        gdalwarp_path + f" -of VRT {testgdalwarp_gcp_tif} tmp/testgdalwarp16.vrt"
     )
 
     ds = gdal.Open("tmp/testgdalwarp16.vrt")
@@ -416,10 +526,10 @@ def test_gdalwarp_18(gdalwarp_path):
 # Test -et 0 which is a special case
 
 
-def test_gdalwarp_19(gdalwarp_path):
+def test_gdalwarp_19(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -et 0 tmp/testgdalwarp_gcp.tif tmp/testgdalwarp19.tif"
+        gdalwarp_path + f" -et 0 {testgdalwarp_gcp_tif} tmp/testgdalwarp19.tif"
     )
 
     ds = gdal.Open("tmp/testgdalwarp19.tif")
@@ -434,10 +544,10 @@ def test_gdalwarp_19(gdalwarp_path):
 # Test -of VRT -et 0 which is a special case
 
 
-def test_gdalwarp_20(gdalwarp_path):
+def test_gdalwarp_20(gdalwarp_path, testgdalwarp_gcp_tif):
 
     gdaltest.runexternal(
-        gdalwarp_path + " -of VRT -et 0 tmp/testgdalwarp_gcp.tif tmp/testgdalwarp20.vrt"
+        gdalwarp_path + f" -of VRT -et 0 {testgdalwarp_gcp_tif} tmp/testgdalwarp20.vrt"
     )
 
     ds = gdal.Open("tmp/testgdalwarp20.vrt")
@@ -1518,109 +1628,3 @@ def test_gdalwarp_if_option(gdalwarp_path):
         gdalwarp_path + " -if HFA ../gcore/data/byte.tif /vsimem/out.tif"
     )
     assert err is not None
-
-
-###############################################################################
-# Cleanup
-
-
-def test_gdalwarp_cleanup():
-
-    # We don't clean up when run in debug mode.
-    if gdal.GetConfigOption("CPL_DEBUG", "OFF") == "ON":
-        return
-
-    for i in range(37):
-        try:
-            os.remove("tmp/testgdalwarp" + str(i + 1) + ".tif")
-        except OSError:
-            pass
-        try:
-            os.remove("tmp/testgdalwarp" + str(i + 1) + ".vrt")
-        except OSError:
-            pass
-        try:
-            os.remove("tmp/testgdalwarp" + str(i + 1) + ".tif.aux.xml")
-        except OSError:
-            pass
-    try:
-        os.remove("tmp/testgdalwarp_gcp.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/testgdalwarp24src.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/testgdalwarp24dst.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/testgdalwarp30_1.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/testgdalwarp30_2.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/testgdalwarp30_3.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/testgdalwarp33_mask.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/testgdalwarp37.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/testgdalwarp38.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/test_gdalwarp_39.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/test_gdalwarp_40_src.tif")
-        os.remove("tmp/test_gdalwarp_40.tif")
-        os.remove("tmp/test_gdalwarp_40.vrt")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/test_gdalwarp_41_src.tif")
-        os.remove("tmp/test_gdalwarp_41.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/small_world_left.tif")
-        os.remove("tmp/small_world_right.tif")
-        os.remove("tmp/test_gdalwarp_42.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/small_world.tif")
-        os.remove("tmp/test_gdalwarp_43.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/test_gdalwarp_44.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/test_gdalwarp_45.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/test_gdalwarp_46.tif")
-    except OSError:
-        pass
-    try:
-        os.remove("tmp/cutline_4326.shp")
-        os.remove("tmp/cutline_4326.shx")
-        os.remove("tmp/cutline_4326.dbf")
-        os.remove("tmp/cutline_4326.prj")
-    except OSError:
-        pass
diff --git a/autotest/utilities/test_gdalwarp_lib.py b/autotest/utilities/test_gdalwarp_lib.py
index 4cd8a4873c5b..568536e69931 100755
--- a/autotest/utilities/test_gdalwarp_lib.py
+++ b/autotest/utilities/test_gdalwarp_lib.py
@@ -40,6 +40,23 @@
 
 from osgeo import gdal, ogr, osr
 
+
+@pytest.fixture(scope="module", autouse=True)
+def setup_and_cleanup():
+
+    yield
+
+    # We don't clean up when run in debug mode.
+    if gdal.GetConfigOption("CPL_DEBUG", "OFF") == "ON":
+        return
+
+    for i in range(2):
+        try:
+            os.remove("tmp/testgdalwarp" + str(i + 1) + ".tif")
+        except OSError:
+            pass
+
+
 ###############################################################################
 # Simple test
 
@@ -104,7 +121,12 @@ def test_gdalwarp_lib_4():
 # Test warping from GCPs without any explicit option
 
 
-def test_gdalwarp_lib_5():
+@pytest.fixture(scope="module")
+def testgdalwarp_gcp_tif(tmp_path_factory):
+
+    testgdalwarp_gcp_tif_fname = str(
+        tmp_path_factory.mktemp("tmp") / "testgdalwarp_gcp.tif"
+    )
 
     ds = gdal.Open("../gcore/data/byte.tif")
     gcpList = [
@@ -114,8 +136,19 @@ def test_gdalwarp_lib_5():
         gdal.GCP(440720.000, 3750120.000, 0, 0, 20),
     ]
     ds1 = gdal.Translate(
-        "tmp/testgdalwarp_gcp.tif", ds, outputSRS="EPSG:26711", GCPs=gcpList
+        testgdalwarp_gcp_tif_fname, ds, outputSRS="EPSG:26711", GCPs=gcpList
     )
+    del ds1
+
+    yield testgdalwarp_gcp_tif_fname
+
+
+def test_gdalwarp_lib_5(testgdalwarp_gcp_tif):
+
+    ds = gdal.Open("../gcore/data/byte.tif")
+
+    ds1 = gdal.Open(testgdalwarp_gcp_tif)
+
     dstDS = gdal.Warp("", ds1, format="MEM", tps=True)
 
     assert dstDS.GetRasterBand(1).Checksum() == 4672, "Bad checksum"
@@ -129,9 +162,9 @@ def test_gdalwarp_lib_5():
 # Test warping from GCPs with -tps
 
 
-def test_gdalwarp_lib_6():
+def test_gdalwarp_lib_6(testgdalwarp_gcp_tif):
 
-    ds1 = gdal.Open("tmp/testgdalwarp_gcp.tif")
+    ds1 = gdal.Open(testgdalwarp_gcp_tif)
     dstDS = gdal.Warp("", ds1, format="MEM", tps=True)
 
     assert dstDS.GetRasterBand(1).Checksum() == 4672, "Bad checksum"
@@ -149,9 +182,9 @@ def test_gdalwarp_lib_6():
 # Test -tr
 
 
-def test_gdalwarp_lib_7():
+def test_gdalwarp_lib_7(testgdalwarp_gcp_tif):
 
-    ds1 = gdal.Open("tmp/testgdalwarp_gcp.tif")
+    ds1 = gdal.Open(testgdalwarp_gcp_tif)
     dstDS = gdal.Warp("", [ds1], format="MEM", xRes=120, yRes=120)
     assert dstDS is not None
 
@@ -165,9 +198,9 @@ def test_gdalwarp_lib_7():
 # Test -ts
 
 
-def test_gdalwarp_lib_8():
+def test_gdalwarp_lib_8(testgdalwarp_gcp_tif):
 
-    ds1 = gdal.Open("tmp/testgdalwarp_gcp.tif")
+    ds1 = gdal.Open(testgdalwarp_gcp_tif)
     dstDS = gdal.Warp("", [ds1], format="MEM", width=10, height=10)
     assert dstDS is not None
 
@@ -368,10 +401,10 @@ def test_gdalwarp_lib_resampling_methods(resampleAlg, resampleAlgStr):
 # Test -dstnodata
 
 
-def test_gdalwarp_lib_15():
+def test_gdalwarp_lib_15(testgdalwarp_gcp_tif):
 
     ds = gdal.Warp(
-        "", "tmp/testgdalwarp_gcp.tif", format="MEM", dstSRS="EPSG:32610", dstNodata=1
+        "", testgdalwarp_gcp_tif, format="MEM", dstSRS="EPSG:32610", dstNodata=1
     )
 
     assert ds.GetRasterBand(1).GetNoDataValue() == 1, "Bad nodata value"
@@ -385,10 +418,10 @@ def test_gdalwarp_lib_15():
 # Test -of VRT which is a special case
 
 
-def test_gdalwarp_lib_16():
+def test_gdalwarp_lib_16(testgdalwarp_gcp_tif):
 
     ds = gdal.Warp(
-        "/vsimem/test_gdalwarp_lib_16.vrt", "tmp/testgdalwarp_gcp.tif", format="VRT"
+        "/vsimem/test_gdalwarp_lib_16.vrt", testgdalwarp_gcp_tif, format="VRT"
     )
     assert ds is not None
 
@@ -402,7 +435,7 @@ def test_gdalwarp_lib_16():
     with pytest.raises(Exception):
         gdal.Warp(
             "/i_dont/exist/test_gdalwarp_lib_16.vrt",
-            "tmp/testgdalwarp_gcp.tif",
+            testgdalwarp_gcp_tif,
             format="VRT",
         )
 
@@ -425,9 +458,9 @@ def test_gdalwarp_lib_17():
 # Test -et 0 which is a special case
 
 
-def test_gdalwarp_lib_19():
+def test_gdalwarp_lib_19(testgdalwarp_gcp_tif):
 
-    ds = gdal.Warp("", "tmp/testgdalwarp_gcp.tif", format="MEM", errorThreshold=0)
+    ds = gdal.Warp("", testgdalwarp_gcp_tif, format="MEM", errorThreshold=0)
     assert ds is not None
 
     assert ds.GetRasterBand(1).Checksum() == 4672, "Bad checksum"
@@ -3481,24 +3514,3 @@ def test_gdalwarp_lib_working_data_type_with_source_dataset_of_different_types()
     minval, maxval = out_ds.GetRasterBand(1).ComputeRasterMinMax()
     assert minval == 1
     assert maxval == 1
-
-
-###############################################################################
-# Cleanup
-
-
-def test_gdalwarp_lib_cleanup():
-
-    # We don't clean up when run in debug mode.
-    if gdal.GetConfigOption("CPL_DEBUG", "OFF") == "ON":
-        return
-
-    for i in range(2):
-        try:
-            os.remove("tmp/testgdalwarp" + str(i + 1) + ".tif")
-        except OSError:
-            pass
-    try:
-        os.remove("tmp/testgdalwarp_gcp.tif")
-    except OSError:
-        pass
diff --git a/autotest/utilities/test_gnmutils.py b/autotest/utilities/test_gnmutils.py
index 4032b9d4459f..5531a3ff6aa4 100755
--- a/autotest/utilities/test_gnmutils.py
+++ b/autotest/utilities/test_gnmutils.py
@@ -46,6 +46,7 @@
         test_cli_utilities.get_gnmanalyse_path() is None,
         reason="gnmanalyse not available",
     ),
+    pytest.mark.random_order(disabled=True),
 ]
 
 
diff --git a/autotest/utilities/test_nearblack.py b/autotest/utilities/test_nearblack.py
index 9d810883c22d..c237e99896d4 100755
--- a/autotest/utilities/test_nearblack.py
+++ b/autotest/utilities/test_nearblack.py
@@ -100,12 +100,7 @@ def test_nearblack_2(nearblack_path):
 
     ds = None
 
-
-###############################################################################
-# Set existing alpha band
-
-
-def test_nearblack_3(nearblack_path):
+    # Set existing alpha band
 
     shutil.copy("tmp/nearblack2.tif", "tmp/nearblack3.tif")
     gdaltest.runexternal(
@@ -165,12 +160,7 @@ def test_nearblack_5(nearblack_path):
 
     ds = None
 
-
-###############################################################################
-# Set existing mask band
-
-
-def test_nearblack_6(nearblack_path):
+    # Set existing mask band
 
     shutil.copy("tmp/nearblack5.tif", "tmp/nearblack6.tif")
     shutil.copy("tmp/nearblack5.tif.msk", "tmp/nearblack6.tif.msk")
diff --git a/autotest/utilities/test_ogrlineref.py b/autotest/utilities/test_ogrlineref.py
index c648c620d8ed..5ca4bfb49368 100755
--- a/autotest/utilities/test_ogrlineref.py
+++ b/autotest/utilities/test_ogrlineref.py
@@ -49,38 +49,36 @@
 ]
 
 
-@pytest.fixture()
+@pytest.fixture(scope="module")
 def ogrlineref_path():
     return test_cli_utilities.get_ogrlineref_path()
 
 
-###############################################################################
-# create test
-
-
-def test_ogrlineref_1(ogrlineref_path):
+@pytest.fixture(scope="module")
+def parts_shp(ogrlineref_path, tmp_path_factory):
 
-    if os.path.exists("tmp/parts.shp"):
-        ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource("tmp/parts.shp")
+    parts_shp_fname = str(tmp_path_factory.mktemp("tmp") / "parts.shp")
 
     _, err = gdaltest.runexternal_out_and_err(
         ogrlineref_path
-        + " -create -l data/path.shp -p data/mstones.shp -pm pos -o tmp/parts.shp -s 1000"
+        + f" -create -l data/path.shp -p data/mstones.shp -pm pos -o {parts_shp_fname} -s 1000"
     )
     assert err is None or err == "", 'got error/warning: "%s"' % err
 
-    ds = ogr.Open("tmp/parts.shp")
-    assert ds is not None and ds.GetLayer(0).GetFeatureCount() == 9
+    with ogr.Open(parts_shp_fname) as ds:
+        assert ds.GetLayer(0).GetFeatureCount() == 9
+
+    yield parts_shp_fname
 
 
 ###############################################################################
 # get_pos test
 
 
-def test_ogrlineref_2(ogrlineref_path):
+def test_ogrlineref_2(ogrlineref_path, parts_shp):
 
     ret = gdaltest.runexternal(
-        ogrlineref_path + " -get_pos -r tmp/parts.shp -x -1.4345 -y 51.9497 -quiet"
+        ogrlineref_path + f" -get_pos -r {parts_shp} -x -1.4345 -y 51.9497 -quiet"
     ).strip()
 
     expected = "15977.724709"
@@ -91,10 +89,10 @@ def test_ogrlineref_2(ogrlineref_path):
 # get_coord test
 
 
-def test_ogrlineref_3(ogrlineref_path):
+def test_ogrlineref_3(ogrlineref_path, parts_shp):
 
     ret = gdaltest.runexternal(
-        ogrlineref_path + " -get_coord -r tmp/parts.shp -m 15977.724709 -quiet"
+        ogrlineref_path + f" -get_coord -r {parts_shp} -m 15977.724709 -quiet"
     ).strip()
 
     expected = "-1.435097,51.950080,0.000000"
@@ -105,14 +103,14 @@ def test_ogrlineref_3(ogrlineref_path):
 # get_subline test
 
 
-def test_ogrlineref_4(ogrlineref_path):
+def test_ogrlineref_4(ogrlineref_path, parts_shp):
 
     if os.path.exists("tmp/subline.shp"):
         ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource("tmp/subline.shp")
 
     gdaltest.runexternal(
         ogrlineref_path
-        + " -get_subline -r tmp/parts.shp -mb 13300 -me 17400 -o tmp/subline.shp"
+        + f" -get_subline -r {parts_shp} -mb 13300 -me 17400 -o tmp/subline.shp"
     )
 
     ds = ogr.Open("tmp/subline.shp")
@@ -129,24 +127,16 @@ def test_ogrlineref_4(ogrlineref_path):
 # test kml
 
 
-def test_ogrlineref_5(ogrlineref_path):
+def test_ogrlineref_5(ogrlineref_path, tmp_path):
+
+    parts_kml = str(tmp_path / "parts.kml")
 
     if os.path.exists("tmp/parts.kml"):
         ogr.GetDriverByName("KML").DeleteDataSource("tmp/parts.kml")
 
     gdaltest.runexternal_out_and_err(
         ogrlineref_path
-        + ' -create -f "KML" -l data/path.shp -p data/mstones.shp -pm pos -o tmp/parts.kml -s 222'
+        + f' -create -f "KML" -l data/path.shp -p data/mstones.shp -pm pos -o {parts_kml} -s 222'
     )
-    if os.path.exists("tmp/parts.kml"):
-        return
 
-    pytest.fail()
-
-
-def test_ogrlineref_cleanup():
-
-    if os.path.exists("tmp/parts.shp"):
-        ogr.GetDriverByName("ESRI Shapefile").DeleteDataSource("tmp/parts.shp")
-    if os.path.exists("tmp/parts.kml"):
-        ogr.GetDriverByName("KML").DeleteDataSource("tmp/parts.kml")
+    assert os.path.exists(parts_kml)
diff --git a/autotest/utilities/test_ogrtindex.py b/autotest/utilities/test_ogrtindex.py
index f557bc0d9530..eafb71f016c2 100755
--- a/autotest/utilities/test_ogrtindex.py
+++ b/autotest/utilities/test_ogrtindex.py
@@ -48,6 +48,22 @@ def ogrtindex_path():
     return test_cli_utilities.get_ogrtindex_path()
 
 
+@pytest.fixture(scope="module", autouse=True)
+def setup_and_cleanup():
+
+    yield
+
+    shape_drv = ogr.GetDriverByName("ESRI Shapefile")
+    if os.path.exists("tmp/tileindex.shp"):
+        shape_drv.DeleteDataSource("tmp/tileindex.shp")
+    shape_drv.DeleteDataSource("tmp/point1.shp")
+    shape_drv.DeleteDataSource("tmp/point2.shp")
+    if os.path.exists("tmp/point3.shp"):
+        shape_drv.DeleteDataSource("tmp/point3.shp")
+    if os.path.exists("tmp/point4.shp"):
+        shape_drv.DeleteDataSource("tmp/point4.shp")
+
+
 ###############################################################################
 # Simple test
 
@@ -246,20 +262,3 @@ def test_ogrtindex_3(ogrtindex_path):
         shape_drv.DeleteDataSource("tmp/tileindex.shp")
     if os.path.exists("tmp/tileindex.db"):
         os.unlink("tmp/tileindex.db")
-
-
-###############################################################################
-# Cleanup
-
-
-def test_ogrtindex_cleanup():
-
-    shape_drv = ogr.GetDriverByName("ESRI Shapefile")
-    if os.path.exists("tmp/tileindex.shp"):
-        shape_drv.DeleteDataSource("tmp/tileindex.shp")
-    shape_drv.DeleteDataSource("tmp/point1.shp")
-    shape_drv.DeleteDataSource("tmp/point2.shp")
-    if os.path.exists("tmp/point3.shp"):
-        shape_drv.DeleteDataSource("tmp/point3.shp")
-    if os.path.exists("tmp/point4.shp"):
-        shape_drv.DeleteDataSource("tmp/point4.shp")
diff --git a/cmake/template/pytest.ini.in b/cmake/template/pytest.ini.in
index a7ea6779511d..72c62aff93c9 100644
--- a/cmake/template/pytest.ini.in
+++ b/cmake/template/pytest.ini.in
@@ -19,6 +19,7 @@ env =
 gdal_version = @GDAL_VERSION_NO_DEV_SUFFIX@
 
 markers =
+    random_order: Indicates whether tests can be run non-sequentially
     require_curl: Skip test(s) if curl support is absent
     require_creation_option: Skip test(s) if required creation option is not available
     require_driver: Skip test(s) if driver isn't present
diff --git a/swig/python/gdal-utils/osgeo_utils/gdal_pansharpen.py b/swig/python/gdal-utils/osgeo_utils/gdal_pansharpen.py
index d27ae35d5bf0..0d8edd42c7dd 100644
--- a/swig/python/gdal-utils/osgeo_utils/gdal_pansharpen.py
+++ b/swig/python/gdal-utils/osgeo_utils/gdal_pansharpen.py
@@ -277,8 +277,13 @@ def gdal_pansharpen(
         ms_name = spectral_ds[i].GetDescription()
         if driver_name.upper() == "VRT":
             if not os.path.isabs(ms_name):
-                ms_relative = "1"
-                ms_name = os.path.relpath(ms_name, os.path.dirname(dst_filename))
+                try:
+                    ms_name = os.path.relpath(ms_name, os.path.dirname(dst_filename))
+                    ms_relative = "1"
+                except ValueError:
+                    # Thrown if generating a relative path is not possible, e.g. if
+                    # ms_name is on a different Windows drive from dst_filename
+                    pass
 
         vrt_xml += """    <SpectralBand%s>
       <SourceFilename relativeToVRT="%s">%s</SourceFilename>