I am wanting to automate the preparation of raster layers from Sentinel-2 to be used in an unsupervised classification. So far, I have been able to write the following script, which will successfully prepare the images in such a way that each band will cover the desired extent even when the extent is covered by multiple satellite image tiles:

bandset = {
    "band02": Band02,
    "band03": Band03,
    "band04": Band04,
    "band05": Band05,
    "band06": Band06,
    "band07": Band07,
    "band08": Band08,
    "band8A": Band8A,
    "band11": Band11,
    "band12": Band12

for bandidx, bandimages in bandset.items():

    # Generate output file
    outputFile = os.path.join(Output_folder, bandidx+".tif")

    # Merge band images together 
    algorithmOutput = processing.runalg("gdalogr:merge", Band_images, False, False, 5, None)
    ret_merged = algorithmOutput['OUTPUT']

    # Reproject the layer to ensure that all layers will have consistent projection and resolution
    algorithmOutput = processing.runalg("gdalogr:warpreproject", ret_merged, "", Target_CRS, "", Target_resolution, 0, 5, 4, 75, 6, 1, False, 0, False, "", None)
    ret_reprojected = algorithmOutput['OUTPUT']

    # Produce buffer around region of interest
    algorithmOutput=processing.runalg("qgis:fixeddistancebuffer", Region_of_Interest, buffer_size_m, 5, False, None, progress=None)

    # Clip the merged band images to the extent of the desired area
    algorithmOutput=processing.runalg("gdalogr:cliprasterbymasklayer", ret_reprojected, ret_buffered, '0.0',  False, True, True, 5, 4, 75, 6, 1, False, 0, False, "", outputFile)        

Then I tried performing the unsupervised classification, in the Semi-Automatic Classification Plugin in QGIS 3.0. I created a band set with the raster layers that are generated by this script and tried running a clustering of the band set using ISODATA when I received the error message "Error [46] : Error reading raster. Possibly bands are not aligned". I thought this was strange because accessing the metadata using the code block pasted below showed that all raster layers have identical size, projection, origin, and pixel size.

for bandidx, bandfiles in bandset.items():
    filename = os.path.abspath(bandfiles[0])
    dataset = gdal.Open(filename, gdal.GA_ReadOnly)

    geotransform = dataset.GetGeoTransform()
    projection = dataset.GetProjection()
    srs = osr.SpatialReference(wkt=projection)

    print("Driver: %s, %s" % (dataset.GetDriver().ShortName, dataset.GetDriver().LongName))
    print("Size: %d x %d x %d" % (dataset.RasterXSize, dataset.RasterYSize, dataset.RasterCount))
    print("Projection: %s" % (srs.GetAttrValue('projcs')))
    print("Origin: (%d, %d)" % (geotransform[0], geotransform[3]))
    print("Pixel Size: %d x %d" % (geotransform[1], geotransform[5]))

However, when I tried aligning the raster layers using the "Raster > Align Rasters" tool before making the band set and then performing the unsupervised classification, there was no error. Thus, aligning the rasters solved my problem completely!

So my question is this: How do I call the Align Rasters tool from a python script, or how would I write a script that accomplishes the same task as the Align Rasters tool?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.