QGIS 3

.gitattributes

@@ -0,0 +1 @@
+*.zip binary

add_to_path.py

@@ -0,0 +1,6 @@
+import os
+import sys
+
+# Required for imports to work in QGIS
+if os.path.dirname(__file__) not in sys.path:
+    sys.path.append(os.path.dirname(__file__))

qgis_fmask/interfaces/fmask_sentinel2Stacked.py

@@ -28,6 +28,7 @@
 
 from fmask import config
 from fmask import fmask
+from fmask.cmdline import sentinel2Stacked
 
 
 def checkAnglesFile(inputAnglesFile, toafile):
@@ -73,6 +74,7 @@ def mainRoutine(cmdargs):
     """
     Main routine that calls fmask
     """
+    topMeta = sentinel2Stacked.readTopLevelMeta(cmdargs)
     anglesfile = checkAnglesFile(cmdargs.anglesfile, cmdargs.toa)
     anglesInfo = config.AnglesFileInfo(
         anglesfile, 3, anglesfile, 2, anglesfile, 1, anglesfile, 0
@@ -87,13 +89,15 @@ def mainRoutine(cmdargs):
     fmaskConfig.setKeepIntermediates(cmdargs.keepintermediates)
     fmaskConfig.setVerbose(cmdargs.verbose)
     fmaskConfig.setTempDir(cmdargs.tempdir)
-    fmaskConfig.setTOARefScaling(10000.0)
+    fmaskConfig.setTOARefScaling(topMeta.scaleVal)
+    offsetDict = sentinel2Stacked.makeRefOffsetDict(topMeta)
+    fmaskConfig.setTOARefOffsetDict(offsetDict)
     fmaskConfig.setMinCloudSize(cmdargs.mincloudsize)
-    fmaskConfig.setEqn17CloudProbThresh(
-        cmdargs.cloudprobthreshold / 100
-    )  # Note conversion from percentage
+    # Note conversion from percentage
+    fmaskConfig.setEqn17CloudProbThresh(float(cmdargs.cloudprobthreshold) / 100)
     fmaskConfig.setEqn20NirSnowThresh(cmdargs.nirsnowthreshold)
     fmaskConfig.setEqn20GreenSnowThresh(cmdargs.greensnowthreshold)
+    fmaskConfig.setSen2displacementTest(cmdargs.parallaxtest)
 
     # Work out a suitable buffer size, in pixels, dependent on the resolution of the input TOA image
     toaImgInfo = fileinfo.ImageInfo(cmdargs.toa)

qgis_fmask/stacks/landsat_stack.py

@@ -1,7 +1,11 @@
 import os
 import glob
 
-from .buildvrt import buildvrt
+import sys
+
+if os.path.dirname(__file__) not in sys.path:
+    sys.path.append(os.path.dirname(__file__))
+import buildvrt as bv
 
 
 bandfile_patterns = {
@@ -16,7 +20,7 @@ def create_landsat_stack(productdir, outfile, landsatkey, imagename):
     infiles = sorted(glob.glob(pattern))
     if not infiles:
         raise RuntimeError("No files found for pattern '{}'.".format(pattern))
-    buildvrt(infiles, outfile, separate=True)
+    bv.buildvrt(infiles, outfile, separate=True)
 
 
 def create_landsat_stacks(productdir, outfile_template, landsatkey):

qgis_fmask/stacks/sentinel_stack.py

@@ -1,24 +1,31 @@
 import os
 import fnmatch
 
-from .buildvrt import buildvrt
+import sys
+
+if os.path.dirname(__file__) not in sys.path:
+    sys.path.append(os.path.dirname(__file__))
+import buildvrt as bv
+
+from fmask.cmdline import sentinel2Stacked
 
 
 def recursive_glob(rootdir, pattern="*"):
     matches = []
-    for root, _, filenames in os.walk(rootdir):
+    for root, _, filenames in os.walk(os.path.join(rootdir, "IMG_DATA")):
         for filename in fnmatch.filter(filenames, pattern):
             matches.append(os.path.join(root, filename))
     return matches
 
 
-def create_sentinel_stack(granuledir, outfile):
+def create_sentinel_stack(safedir, outfile):
+    granuledir = sentinel2Stacked.findGranuleDir(safedir)
     infiles = []
     for fnpattern in ["*_B0[1-8].jp2", "*_B8A.jp2", "*_B09.jp2", "*_B1[0-2].jp2"]:
         bandfiles = recursive_glob(rootdir=granuledir, pattern=fnpattern)
         if not bandfiles:
             raise RuntimeError("No files found for pattern '{}'.".format(fnpattern))
         infiles += bandfiles
-    buildvrt(
+    bv.buildvrt(
         infiles, outfile, resolution="user", separate=True, extra=["-tr", "20", "20"]
     )

scripts/qgis_fmask_burnMask.py

@@ -1,56 +1,129 @@
-# Definition of inputs and outputs
-# ==================================
-##Fmask=group
-##Burn cloud mask=name
-##ParameterRaster|dataFile|An image file to burn the cloud mask into|False
-##ParameterRaster|maskFile|Could mask from FMask|False
-##*ParameterBoolean|maskNull|Mask FMask null pixels|True
-##*ParameterBoolean|maskCloud|Mask FMask cloud pixels|True
-##*ParameterBoolean|maskShadow|Mask FMask shadow pixels|True
-##*ParameterBoolean|maskSnow|Mask FMask snow pixels|True
-##*ParameterBoolean|maskWater|Mask FMask water pixels|False
-##*ParameterBoolean|maskLand|Mask FMask land pixels|False
-##OutputRaster|outputFile|Maked output image
+# TODO yields Error1: PROJ...
+from qgis.processing import alg
+from processing.core.ProcessingConfig import ProcessingConfig
+from db_manager.db_plugins import data_model
+from processing.tools import system
+from qgis import processing
 
-from processing.tools import dataobjects, system
 
-# Run GDAL warp to make sure that the mask file exactly aligns with the image file
-progress.setText("Aligning mask raster to image raster...")
-dataRaster = dataobjects.getObject(dataFile)
-proj = dataRaster.crs().authid()
-resolution = dataRaster.rasterUnitsPerPixelX()
-bandCount = dataRaster.bandCount()
-extent = dataobjects.extent([dataRaster])
-warpMask = system.getTempFilename("tif")
-params = {
-    "INPUT": maskFile,
-    "DEST_SRS": proj,
-    "TR": resolution,
-    "USE_RASTER_EXTENT": True,
-    "RASTER_EXTENT": extent,
-    "EXTENT_CRS": proj,
-    "METHOD": 0,
-    "RTYPE": 0,
-    "OUTPUT": warpMask,
-}
-processing.runalg("gdalogr:warpreproject", params)
+@alg(
+    name="burncloudmask",
+    label=alg.tr("Burn cloud mask"),
+    group="fmask",
+    group_label=alg.tr("FMask"),
+)
+@alg.input(
+    type=alg.RASTER_LAYER,
+    name="dataFile",
+    label="An image file to burn the cloud mask into",
+    optional=False,
+)
+@alg.input(
+    type=alg.RASTER_LAYER,
+    name="maskFile",
+    label="Cloud mask from FMask",
+    optional=False,
+)
+@alg.input(
+    type=bool,
+    name="maskNull",
+    label="Mask FMask null pixels",
+    default=True,
+    advanced=True,
+)
+@alg.input(
+    type=bool,
+    name="maskCloud",
+    label="Mask FMask cloud pixels",
+    default=True,
+    advanced=True,
+)
+@alg.input(
+    type=bool,
+    name="maskShadow",
+    label="Mask FMask shadow pixels",
+    default=True,
+    advanced=True,
+)
+@alg.input(
+    type=bool,
+    name="maskSnow",
+    label="Mask FMask snow pixels",
+    default=True,
+    advanced=True,
+)
+@alg.input(
+    type=bool,
+    name="maskWater",
+    label="Mask FMask water pixels",
+    default=False,
+    advanced=True,
+)
+@alg.input(
+    type=bool,
+    name="maskLand",
+    label="Mask FMask land pixels",
+    default=False,
+    advanced=True,
+)
+@alg.input(type=alg.FILE_DEST, name="outputFile", label="Maked output image")
+def burncloudmask(instance, parameters, context, feedback, inputs):
+    """burncloudmask"""
 
-progress.setText("Applying mask to image...")
-# First reclassify fmask output into two classes
-flags = [maskNull, maskLand, maskCloud, maskShadow, maskSnow, maskWater]
-maskClass = [str(i) for i in range(len(flags)) if flags[i]]
-leaveClass = [str(i) for i in range(len(flags)) if not flags[i]]
-reclassString = " ".join(maskClass) + " = 0\n" + " ".join(leaveClass) + " = 1"
-reclassMask = system.getTempFilename("tif")
-params = {
-    "input": warpMask,
-    "txtrules": reclassString,
-    "output": reclassMask,
-    "GRASS_REGION_PARAMETER": extent,
-    "GRASS_REGION_CELLSIZE_PARAMETER": resolution,
-}
-processing.runalg("grass7:r.reclass", params)
+    dataRaster = instance.parameterAsRasterLayer(parameters, "dataFile", context)
+    # Run GDAL warp to make sure that the mask file exactly aligns with the image file
 
-# Then use OTB band maths on all bands
-params = {"-il": dataFile + ";" + reclassMask, "-exp": "im1*im2b1", "-out": outputFile}
-processing.runalg("otb:bandmathx", params)
+    feedback.setProgressText("Aligning mask raster to image raster...")
+    # feedback.setProgressText(str(dir(dataRaster)))
+
+    feedback.setProgressText(str(dir(parameters["dataFile"])))
+    extent = dataRaster.extent()
+    proj = dataRaster.crs().toWkt()
+    resolution = dataRaster.rasterUnitsPerPixelX()
+    bandCount = dataRaster.bandCount()
+    warpMask = system.getTempFilename("tif")
+    params = {
+        "INPUT": parameters["maskFile"],
+        "TARGET_CRS": proj,
+        "TARGET_RESOLUTION": resolution,
+        "TARGET_EXTENT": extent,
+        "TARGET_EXTENT_CRS": proj,
+        "RESAMPLING": 0,
+        "DATA_TYPE": 0,
+        "OUTPUT": warpMask,
+    }
+    processing.run("gdal:warpreproject", params, feedback=feedback, context=context)
+
+    feedback.setProgressText("Applying mask to image...")
+    # First reclassify fmask output into two classes
+    flags = [
+        parameters["maskNull"],
+        parameters["maskLand"],
+        parameters["maskCloud"],
+        parameters["maskShadow"],
+        parameters["maskSnow"],
+        parameters["maskWater"],
+    ]
+    maskClass = [str(i) for i in range(len(flags)) if flags[i]]
+    leaveClass = [str(i) for i in range(len(flags)) if not flags[i]]
+    reclassString = " ".join(maskClass) + " = 0\n" + " ".join(leaveClass) + " = 1"
+    reclassMask = system.getTempFilename("tif")
+    params = {
+        "input": warpMask,
+        "txtrules": reclassString,
+        "output": reclassMask,
+        "GRASS_REGION_PARAMETER": extent,
+        "GRASS_REGION_CELLSIZE_PARAMETER": resolution,
+    }
+    processing.run("grass7:r.reclass", params, feedback=feedback, context=context)
+
+    # Then use OTB band maths on all bands
+    params = {
+        "il": [
+            dataRaster.source(),
+            reclassMask,
+        ],
+        "exp": "im1*im2b1",
+        "out": instance.parameterAsString(parameters, "outputFile", context),
+    }
+    processing.run("otb:BandMathX", params, feedback=feedback, context=context)

scripts/qgis_fmask_reclassifyLandsatQA.py

@@ -0,0 +1,72 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Feb 25 09:15:47 2022
+
+@author: rmgu
+"""
+
+from qgis.processing import alg
+import numpy as np
+from osgeo import gdal
+from processing.tools import dataobjects
+
+_clear_land_fmask = 0
+_clear_water_fmask = 1
+_cloud_shadow_fmask = 2
+_snow_fmask = 3
+_cloud_fmask = 4
+_filled_fmask = 255
+
+
+@alg(
+    name="reclassifylandsatqa",
+    label=alg.tr("Reclassify Landsat Quality Assesment image"),
+    group="fmask",
+    group_label=alg.tr("FMask"),
+)
+@alg.input(
+    type=alg.FILE,
+    name="qaImage",
+    label="Landsat Level-2 QA_PIXEL file",
+    optional=False,
+)
+@alg.input(
+    type=alg.RASTER_LAYER_DEST, name="outputFile", label="Reclassified output image"
+)
+def reclassifyLandsatQa(instance, parameters, context, feedback, inputs):
+    """Reclassify Landsat QA_PIXEL image to Fmask classes"""
+
+    def check_bit_flag(flag_array, bitmask):
+        return np.bitwise_and(qaRaster, bitmask).astype(bool)
+
+    qaRasterFile = instance.parameterAsFile(parameters, "qaImage", context)
+    feedback.setProgressText(qaRasterFile)
+    qaRasterImage = gdal.Open(qaRasterFile)
+    qaRaster = qaRasterImage.GetRasterBand(1).ReadAsArray()
+
+    fmaskClassRaster = np.zeros(qaRaster.shape)
+
+    # See https://www.usgs.gov/media/images/landsat-collection-2-pixel-quality-assessment-bit-index
+    # for bit flags
+    fmaskClassRaster[check_bit_flag(qaRaster, 64)] = _clear_water_fmask
+    fmaskClassRaster[check_bit_flag(qaRaster, 128)] = _clear_land_fmask
+    fmaskClassRaster[check_bit_flag(qaRaster, 1)] = _filled_fmask
+    fmaskClassRaster[check_bit_flag(qaRaster, 14)] = _cloud_fmask
+    fmaskClassRaster[check_bit_flag(qaRaster, 16)] = _cloud_shadow_fmask
+    fmaskClassRaster[check_bit_flag(qaRaster, 32)] = _snow_fmask
+    drv = gdal.GetDriverByName("GTiff")
+    outputFile = instance.parameterAsFileOutput(parameters, "outputFile", context)
+    outTif = drv.Create(
+        outputFile,
+        qaRasterImage.RasterXSize,
+        qaRasterImage.RasterYSize,
+        1,
+        gdal.GDT_Int16,
+    )
+    outTif.SetProjection(qaRasterImage.GetProjection())
+    outTif.SetGeoTransform(qaRasterImage.GetGeoTransform())
+    outTif.GetRasterBand(1).WriteArray(fmaskClassRaster)
+
+    outTif = None
+    qaRasterImage = None
+    dataobjects.load(outputFile, isRaster=True)