diff --git a/autotest/pymod/__init__.py b/autotest/pymod/__init__.py
new file mode 100644
index 000000000000..e69de29bb2d1
diff --git a/autotest/pyscripts/test_gdal_retile.py b/autotest/pyscripts/test_gdal_retile.py
index afd948705b94..adc1cebc6246 100755
--- a/autotest/pyscripts/test_gdal_retile.py
+++ b/autotest/pyscripts/test_gdal_retile.py
@@ -31,7 +31,6 @@
 import shutil
 import os
 
-
 from osgeo import gdal
 from osgeo import osr
 import test_py_scripts
@@ -241,7 +240,66 @@ def test_gdal_retile_4():
         assert ds.RasterYSize == height, filename
         ds = None
 
-    
+
+###############################################################################
+# Test gdal_retile.py with input having a NoData value
+
+
+def test_gdal_retile_5():
+
+    try:
+        from osgeo import gdalnumeric
+        gdalnumeric.zeros
+        import numpy as np
+    except (ImportError, AttributeError):
+        pytest.skip()
+
+    nodata_value = -3.4028234663852886e+38
+    raster_array = np.array(([0.0, 2.0], [-1.0, nodata_value]))
+
+    script_path = test_py_scripts.get_py_script('gdal_retile')
+    if script_path is None:
+        pytest.skip()
+
+    drv = gdal.GetDriverByName('GTiff')
+    srs = osr.SpatialReference()
+    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
+
+    try:
+        os.mkdir('tmp/outretile5')
+    except OSError:
+        pass
+
+    test_py_scripts.run_py_script(script_path, 'gdal_retile', '-v -targetDir tmp/outretile5 tmp/in5.tif')
+
+    ds = gdal.Open('tmp/outretile5/in5_1_1.tif')
+    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 min_val == -1.0, \
+        ('Wrong minimum value.\nExpected -1.0, Got: %f' % min_val)
+
+    ds = None
+
+
 ###############################################################################
 # Cleanup
 
@@ -267,7 +325,8 @@ def test_gdal_retile_cleanup():
            'tmp/outretile3/1',
            'tmp/outretile3/2',
            'tmp/outretile3/in1_1_1.tif',
-           'tmp/outretile3']
+           'tmp/outretile3',
+           'tmp/in5.tif']
     for filename in lst:
         try:
             os.remove(filename)
@@ -279,6 +338,6 @@ def test_gdal_retile_cleanup():
 
     shutil.rmtree('tmp/outretile4')
 
-
-
+    if os.path.exists('tmp/outretile5'):
+        shutil.rmtree('tmp/outretile5')
 
diff --git a/gdal/swig/python/scripts/gdal_retile.py b/gdal/swig/python/scripts/gdal_retile.py
index 9ff131583797..d2b8091dbd99 100755
--- a/gdal/swig/python/scripts/gdal_retile.py
+++ b/gdal/swig/python/scripts/gdal_retile.py
@@ -148,6 +148,7 @@ def __init__(self, filename, inputDS):
         self.bands = fhInputTile.RasterCount
         self.band_type = fhInputTile.GetRasterBand(1).DataType
         self.projection = fhInputTile.GetProjection()
+        self.nodata = fhInputTile.GetRasterBand(1).GetNoDataValue()
 
         dec = AffineTransformDecorator(fhInputTile.GetGeoTransform())
         self.scaleX = dec.scaleX
@@ -209,6 +210,12 @@ def getDataSet(self, minx, miny, maxx, maxy):
         resultDS = self.TempDriver.Create("TEMP", resultSizeX, resultSizeY, self.bands, self.band_type, [])
         resultDS.SetGeoTransform([minx, self.scaleX, 0, maxy, 0, self.scaleY])
 
+        for bandNr in range(1, self.bands + 1):
+            t_band = resultDS.GetRasterBand(bandNr)
+            if self.nodata is not None:
+                t_band.Fill(self.nodata)
+                t_band.SetNoDataValue(self.nodata)
+
         for feature in features:
             featureName = feature.GetField(0)
             sourceDS = self.cache.get(featureName)
@@ -240,8 +247,6 @@ def getDataSet(self, minx, miny, maxx, maxy):
             tw_yoff = int((tgw_uly - maxy) / self.scaleY + 0.5)
             tw_xsize = min(resultDS.RasterXSize, int((tgw_lrx - minx) / self.scaleX + 0.5)) - tw_xoff
             tw_ysize = min(resultDS.RasterYSize, int((tgw_lry - maxy) / self.scaleY + 0.5)) - tw_yoff
-            # print(sw_xoff, sw_yoff, sw_xsize, sw_ysize, sourceDS.RasterXSize, sourceDS.RasterYSize)
-            # print(tw_xoff, tw_yoff, tw_xsize, tw_ysize, resultDS.RasterXSize, resultDS.RasterYSize)
             if tw_xsize <= 0 or tw_ysize <= 0:
                 continue
 
@@ -456,6 +461,12 @@ def createPyramidTile(levelMosaicInfo, offsetX, offsetY, width, height, tileName
             t_band.SetRasterColorTable(levelMosaicInfo.ct)
         t_band.SetRasterColorInterpretation(levelMosaicInfo.ci[band - 1])
 
+    if levelMosaicInfo.nodata is not None:
+        for band in range(1, bands + 1):
+            t_band = t_fh.GetRasterBand(band)
+            t_band.Fill(levelMosaicInfo.nodata)
+            t_band.SetNoDataValue(levelMosaicInfo.nodata)
+
     res = gdal.ReprojectImage(s_fh, t_fh, None, None, ResamplingMethod)
     if res != 0:
         print("Reprojection failed for %s, error %d" % (tileName, res))
@@ -522,8 +533,10 @@ def createTile(minfo, offsetX, offsetY, width, height, tilename, OGRDS):
         t_band = t_fh.GetRasterBand(band)
         if minfo.ct is not None:
             t_band.SetRasterColorTable(minfo.ct)
+        if minfo.nodata is not None:
+            t_band.Fill(minfo.nodata)
+            t_band.SetNoDataValue(minfo.nodata)
 
-#        data = s_band.ReadRaster( offsetX,offsetY,width,height,width,height, t_band.DataType )
         data = s_band.ReadRaster(0, 0, readX, readY, readX, readY, t_band.DataType)
         t_band.WriteRaster(0, 0, readX, readY, data, readX, readY, t_band.DataType)