Simplify conversion from NetCDF to multiband GeoTIFF

SEPTREFHUM.py

@@ -108,19 +108,14 @@ for timestep in timesteps:
     timestep_dates.append(file_date)
     # Create NetCDF result file
     subprocess.run('cdo -seltimestep,%s/%s %sresult.nc %sresult_%s.nc' % (timestep_index, timestep_index, tmpfile_path, tmpfile_path, file_date), shell=True)
-    # Split all variables into separate GeoTIFF files, using GDAL
-    subprocess.run("gdal_translate -of GTiff -sds %sresult_%s.nc %stmp.tif" % (tmpfile_path, file_date, tmpfile_path) , shell=True)
-    # The variables are: 1 = time_bnds, 2= lon, 3= lat, 4= WHS, 5 = WARNING_STATUS
+    # Convert to GeoTIFF
     # We only need WHS and WARNING_STATUS
-    # Merge the WARNING_STATUS and WHS GeoTIFF files into one file with two bands. The WARNING_STATUS should always be
-    # Band #1
-    subprocess.run("gdal_merge.py -separate -o %sresult_%s.tif %stmp_5.tif %stmp_4.tif " % (outfile_path, file_date,  tmpfile_path, tmpfile_path) , shell=True)
-    subprocess.run("rm %stmp_*.tif" % tmpfile_path, shell=True)
+    # Merge the WARNING_STATUS and WHS GeoTIFF files into one file with two bands. 
+    # The WARNING_STATUS should always be band #1
+    subprocess.run(f'gdal_merge.py -separate -o {outfile_path}result_{file_date}.tif NETCDF:"{tmpfile_path}result_{file_date}.nc":WARNING_STATUS NETCDF:"{tmpfile_path}result_{file_date}.nc":WHS', shell=True)
+    
     timestep_index = timestep_index + 1
 
-
-
-
 # Generate mapfile
 env = Environment(loader=FileSystemLoader('.'))
 template = env.get_template("mapfile/template.j2")
@@ -132,4 +127,4 @@ with open("mapfile/SEPTREFHUM.map", 'w') as f:
 
 
 # Remove all temporary/intermediary files
-subprocess.run("rm %s*" % tmpfile_path, shell=True)
+subprocess.run(f"rm {tmpfile_path}*", shell=True)
\ No newline at end of file