2

I am using SAGA GIS and python modules in QGIS to compute the Topographic Position Index (TPI) for a group of Digital Elevation Models (DEMs) and I need to export each computed TPI in a separate raster file. How can I do it? See my written code below:

base_dir = 'D:/Selected DEMs'
fnames = [os.path.join(base_dir, fname) for fname in os.listdir(base_dir)]
dfs_names = (os.listdir(base_dir))

for i in range(len (fnames)):
    print (i)
    tpi = processing.run("saga:topographicpositionindextpi", {'DEM':fnames[i],'TPI':'TEMPORARY_OUTPUT','STANDARD':False,'DW_WEIGHTING':0,'DW_IDW_POWER':2,'DW_BANDWIDTH':75})
hamid mohebzadeh
  • 127
  • 5

2 Answers2

0

I found the answer. QGIS saves the computed TPI with SAGA GIS format (sdat) in a temp folder. The computed raster is in dict format as below:

dict_items([('TPI', 'C:/Users/.../.../TPI.sdat')])).

Therefore, firstly it is required to read the TPI.sdat from the dict and then export it to Geotiff using rasterio as below:

import rasterio as rio

base_dir = 'D:/Selected DEMs'
fnames = [os.path.join(base_dir, fname) for fname in os.listdir(base_dir)]
dfs_names = (os.listdir(base_dir))

for i in range(len (fnames)):
    print (i)
    tpi = processing.run("saga:topographicpositionindextpi", {'DEM':fnames[i],'TPI':'TEMPORARY_OUTPUT','STANDARD':False,'DW_WEIGHTING':0,'DW_IDW_POWER':2,'DW_BANDWIDTH':75})
    tpi.items()
    tpi['TPI']
    TPI_sdat_format = rio.open(tpi['TPI'])
    TPI_sdat_format_meta = TPI_sdat_format.profile   
    TPI_sdat_format_meta['driver'] = 'GTiff'
    TPI_data=TPI_sdat_format.read(1)
    
    with rio.open('D:/TPI_%s.tif'%i, 'w', **TPI_sdat_format_meta) as dst:
        dst.write(TPI_data, 1)
hamid mohebzadeh
  • 127
  • 5
0

Because QGIS is used, I would use the modules available in QGIS. First the file is written out as .sdat and then loaded back into QGIS as a raster layer. Then the raster can be written as a .tif file

from qgis.core import *
import os

base_dir = 'D:/Selected DEMs'
out_dir = 'D:/Precessed_DEMs'
fnames = os.listdir(base_dir)
EPSG = 'EPSG:25832' # example crs or set from layer r

for i in fnames:
    print (i)
    out_filename =  os.path.join(out_dir,os.path.splitext(i)[0] + '.sdat')
    tpi = processing.run("saga:topographicpositionindextpi",{'DEM':i,'TPI':out_filename ,'STANDARD':False,'DW_WEIGHTING':0,'DW_IDW_POWER':2,'DW_BANDWIDTH':75})

    rlayer = QgsRasterLayer(os.path.join(folder_preproc, os.path.splitext(i)[0]+".sdat"),os.path.splitext(i)[0])
    pipe = QgsRasterPipe()
    pipe.set(rlayer.dataProvider().clone())
    file_writer = QgsRasterFileWriter(os.path.join(folder_preproc, os.path.splitext(i)[0]+".tif"))
    file_writer.writeRaster(pipe, rlayer.width(),rlayer.height(), rlayer.extent(), QgsCoordinateReferenceSystem(EPSG))
Muesgen
  • 152
  • 8