Skip to main content
2 of 3
naming
Vince
  • 20.3k
  • 16
  • 48
  • 65

Get source and output paths in QGIS 3 processing algorithm

How can I find the complete path of the INPUT (r_in_name) and OUTPUT (r_out_name) rasters for the gdal datasets in the following code:

from qgis.processing import alg
from qgis.core import QgsFeature, QgsFeatureSink
import numpy as np
from osgeo import gdal, gdal_array, osr
@alg(name="quantilizearaster", label=alg.tr("Quantilize a raster)"), group="rastertools", group_label=alg.tr("Raster Tools"))
@alg.input(type=alg.RASTER_LAYER, name="INPUT", label="Input layer")
@alg.input(type=alg.RASTER_LAYER, name="OUTPUT", label="Output layer")
def quantilize(instance, parameters, context, feedback, inputs):
    """
    This algorithm takes an INPUT raster and quantilize its values in five classes; result is a quantilized OUTPUT raster.
    """

source = instance.parameterAsRasterLayer(parameters, "INPUT", context)
output_raster = instance.parameterAsRasterLayer(parameters, "OUTPUT", context)

r_in_name =   # r_in_name must be a string of the complete path to the INPUT raster, including its file name (*.tif)
r_out_name =  # r_out_name may be a temporary OUTPUT raster file
    

dataset = gdal.Open(r_in_name)    
band = dataset.GetRasterBand(1)
nodata = band.GetNoDataValue()
array = dataset.ReadAsArray()

new_array = array
nan_array = array

nan_array[array == nodata] = np.nan

percentile_80 = np.nanpercentile(nan_array, 80)
percentile_60 = np.nanpercentile(nan_array, 60)
percentile_40 = np.nanpercentile(nan_array, 40)
percentile_20 = np.nanpercentile(nan_array, 20)
percentile_0 = np.nanpercentile(nan_array, 0)

for i, v in enumerate(new_array):
    for j, element in enumerate(v):
        if element > percentile_0 and element <= percentile_20:
            new_array[i,j] = 1
        if element > percentile_20 and element <= percentile_40:
            new_array[i,j] = 2
        if element > percentile_40 and element <= percentile_60:
            new_array[i,j] = 3
        if element > percentile_60 and element <= percentile_80:
            new_array[i,j] = 4
        if element > percentile_80:
            new_array[i,j] = 5

new_array[ new_array != new_array ] = nodata

geotransform = dataset.GetGeoTransform()
wkt = dataset.GetProjection()

# Create gtif file

driver = gdal.GetDriverByName("GTiff")
output_file = r_out_name  # r_out_name may be a temporary OUTPUT raster file

dst_ds = driver.Create(output_file,
                       band.XSize,
                       band.YSize,
                       1,
                       gdal.GDT_Int16)

#writting output raster
dst_ds.GetRasterBand(1).WriteArray( new_array )
#setting nodata value
dst_ds.GetRasterBand(1).SetNoDataValue(nodata)
#setting extension of output raster
# top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
dst_ds.SetGeoTransform(geotransform)
# setting spatial reference of output raster
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#Close output raster dataset

dataset = None
dst_ds = None 

return 
MLourdes
  • 318
  • 2
  • 6