fire2a.raster.gdal_calc_norm

  1#!python
  2# fmt: off
  3"""
  4<pre><code>
  5<!-- BEGIN_ARGPARSE_DOCSTRING -->
  6usage: gdal_calc_norm.py [-h] [-i INFILE] [-o OUTFILE]
  7                         [-m {minmax,maxmin,stepup,stepdown,bipiecewiselinear,bipiecewiselinear_percent,stepdown_percent,stepup_percent}]
  8                         [-min MINIMUM] [-max MAXIMUM] [-n [NODATAVALUE]] [-f FORMAT]
  9                         [-t {Byte,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64,CInt16,CInt32,CFloat32,CFloat64}]
 10                         [-p min_x max_y max_x min_y] [-r]
 11                         [params ...]
 12
 13Raster normalization utility, wrapping on osgeo_utils.gdal_calc with a set of predefined normalization methods. Run `gdal_calc.py
 14--help` for more information.
 15
 16positional arguments:
 17  params                Float numbers according to the normalizing method. None: for minmax, maxmin; One for stepup, stepdown; Two for
 18                        bipiecewiselinear, bipiecewiselinear_percent, see also method (default: None)
 19
 20options:
 21  -h, --help            show this help message and exit
 22  -i INFILE, --infile INFILE
 23                        Input file (default: infile.tif)
 24  -o OUTFILE, --outfile OUTFILE
 25                        Output file (default: outfile.tif)
 26  -m {minmax,maxmin,stepup,stepdown,bipiecewiselinear,bipiecewiselinear_percent,stepdown_percent,stepup_percent}, --method {minmax,maxmin,stepup,stepdown,bipiecewiselinear,bipiecewiselinear_percent,stepdown_percent,stepup_percent}
 27                        Method to normalize the input, see also params (default: minmax)
 28  -min MINIMUM, --minimum MINIMUM
 29                        Minimun value for minmax/maxmin (else the value is calculated from the WHOLE input raster). (default: None)
 30  -max MAXIMUM, --maximum MAXIMUM
 31                        Maximum value for minmax/maxmin (else the value is calculated from the WHOLE input raster). (default: None)
 32  -n [NODATAVALUE], --NoDataValue [NODATAVALUE]
 33                        output nodata value (send null for default datatype specific, see `from osgeo_utils.gdal_calc import
 34                        DefaultNDVLookup`) (default: -9999)
 35  -f FORMAT, --format FORMAT
 36                        Output format (send null to infer from file, see `from osgeo_utils.auxiliary.util import GetOutputDriverFor`)
 37                        (default: GTiff)
 38  -t {Byte,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64,CInt16,CInt32,CFloat32,CFloat64}, --type {Byte,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64,CInt16,CInt32,CFloat32,CFloat64}
 39                        Output datatype (default: Float32)
 40  -p min_x max_y max_x min_y, --projwin min_x max_y max_x min_y
 41                        An optional list of 4 coordinates defining the projection window, if not provided the whole raster is
 42                        calculated (default: None)
 43  -r, --return_dataset  Return dataset (for scripting -additional keyword arguments are passed to gdal_calc.Calc) instead of return
 44                        code (default: False)
 45
 46documentation at https://fire2a.github.io/fire2a-lib/fire2a/raster/gdal_calc_norm.html
 47<!-- END_ARGPARSE_DOCSTRING -->
 48
 49Sample script usage:
 50    from fire2a.raster.gdal_calc_norm import main
 51    ds = main(["-r", ... other keyword arguments are passed to gdal_calc.Calc
 52
 53function osgeo_utils.gdal_calc.Calc(
 54    calc: Union[str, Sequence[str]],
 55    outfile: Union[str, os.PathLike, NoneType] = None,
 56    NoDataValue: Optional[numbers.Number] = None,
 57    type: Union[int, str, NoneType] = None,
 58    format: Optional[str] = None,
 59    creation_options: Optional[Sequence[str]] = None,
 60    allBands: str = '',
 61    overwrite: bool = False,
 62    hideNoData: bool = False,
 63    projectionCheck: bool = False,
 64    color_table: Union[osgeo.gdal.ColorTable,
 65    ForwardRef('ColorPalette'), str, os.PathLike, Sequence[str], NoneType] = None,
 66    extent: Optional[osgeo_utils.auxiliary.extent_util.Extent] = None,
 67    projwin: Union[Tuple, osgeo_utils.auxiliary.rectangle.GeoRectangle, NoneType] = None,
 68    user_namespace: Optional[Dict] = None,
 69    debug: bool = False,
 70    quiet: bool = False, **input_files)
 71
 72/usr/bin/gdal_calc.py
 73/usr/lib/python3/dist-packages/osgeo_utils/gdal_calc.py
 74</code></pre>
 75"""
 76import sys
 77
 78if sys.version_info < (3, 10):
 79    raise ImportError("This module requires Python 3.10 or higher")
 80
 81from pathlib import Path
 82from tempfile import NamedTemporaryFile
 83
 84from osgeo.gdal import Dataset, GA_ReadOnly, Open
 85
 86try:
 87    from osgeo_utils.auxiliary.util import GetOutputDriverFor
 88    from osgeo_utils.gdal_calc import Calc, GDALDataTypeNames
 89except ImportError:
 90    from osgeo.utils.gdal_calc import Calc, GetOutputDriverFor
 91
 92    GDALDataTypeNames = (
 93        "Byte",
 94        "UInt16",
 95        "Int16",
 96        "UInt32",
 97        "Int32",
 98        "UInt64",
 99        "Int64",
100        "Float32",
101        "Float64",
102        "CInt16",
103        "CInt32",
104        "CFloat32",
105        "CFloat64",
106    )
107
108
109def calc(
110    func,
111    outfile="outfile.tif",
112    infile="infile.tif",
113    band=1,
114    NoDataValue=None,
115    overwrite=True,
116    type="Float32",
117    format="GTiff",
118    projwin=None,
119    minimum=None,
120    maximum=None,
121    threshold=None,
122    a=None,
123    b=None,
124    r=None,
125    **kwargs,
126) -> Dataset:
127    """This is the wrapper function for the gdal_calc.Calc utility.
128
129    Although the normalization functions are symbolic strings defined in the main, passed as the `func` argument. Make sure to check them before using this function directly.
130
131    All extra keyword arguments are passed to the gdal_calc.Calc function."""
132    if isinstance(infile, Path):
133        infile = str(infile)
134    if isinstance(outfile, Path):
135        outfile = str(outfile)
136    if "method" in kwargs:
137        if kwargs["method"] in ["minmax", "maxmin", "bipiecewiselinear_percent", "stepup_percent", "stepdown_percent"]:
138            if minimum is None and maximum is None:
139                minimum, maximum = get_file_minmax(infile)
140                print(f"{minimum=}, {maximum=}")
141            if minimum is None:
142                minimum, _ = get_file_minmax(infile)
143                print(f"{minimum=}")
144            if maximum is None:
145                _, maximum = get_file_minmax(infile)
146                print(f"{maximum=}")
147        # drop before passing to Calc
148        del kwargs["method"]
149    # if not projwin:
150    #     info = locals().get("info", read_raster(infile, data=False)[1])
151    #     projwin = get_projwin(info["Transform"], info["RasterXSize"], info["RasterYSize"])
152    #     print(f"{projwin=}")
153
154    # get a the parameter values of the calc function
155    code_obj = func.__code__
156    local_var_names = code_obj.co_varnames[: code_obj.co_nlocals]
157    func_vars = []
158    for var in local_var_names:
159        func_vars += [locals().get(var)]
160    if "r" in local_var_names:
161        func_vars[-1] = (maximum - minimum) / 100
162    print(f"{dict(zip(local_var_names,func_vars))=}")
163    print(f"{local_var_names=}")
164    print(f"{func_vars=}")
165
166    dataset = Calc(
167        calc=func(*func_vars),
168        outfile=outfile,
169        A=infile,
170        A_band=band,
171        NoDataValue=NoDataValue,
172        overwrite=overwrite,
173        type=type,
174        format=format,
175        projwin=projwin,
176        **kwargs,
177    )
178    dataset.FlushCache()
179
180    return dataset
181
182
183def arg_parser(argv=None):
184    """Parse arguments list"""
185    from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
186
187    parser = ArgumentParser(
188        description="Raster normalization utility, wrapping on osgeo_utils.gdal_calc with a set of predefined normalization methods. Run `gdal_calc.py --help` for more information.",
189        formatter_class=ArgumentDefaultsHelpFormatter,
190        epilog="documentation at https://fire2a.github.io/fire2a-lib/fire2a/raster/gdal_calc_norm.html",
191    )
192    parser.add_argument(
193        "params",
194        nargs="*",
195        type=float,
196        help="Float numbers according to the normalizing method. None: for minmax, maxmin; One for stepup, stepdown; Two for bipiecewiselinear, bipiecewiselinear_percent, see also method",
197    )
198    parser.add_argument("-i", "--infile", help="Input file", type=Path, default="infile.tif")
199    parser.add_argument("-o", "--outfile", help="Output file", type=Path, default="outfile.tif")
200    parser.add_argument(
201        "-m",
202        "--method",
203        help="Method to normalize the input, see also params",
204        type=str,
205        choices=[
206            "minmax",
207            "maxmin",
208            "stepup",
209            "stepdown",
210            "bipiecewiselinear",
211            "bipiecewiselinear_percent",
212            "stepdown_percent",
213            "stepup_percent",
214        ],
215        default="minmax",
216    )
217    parser.add_argument(
218        "-min",
219        "--minimum",
220        help="Minimun value for minmax/maxmin (else the value is calculated from the WHOLE input raster).",
221        type=float,
222    )
223    parser.add_argument(
224        "-max",
225        "--maximum",
226        help="Maximum value for minmax/maxmin (else the value is calculated from the WHOLE input raster).",
227        type=float,
228    )
229    parser.add_argument(
230        "-n",
231        "--NoDataValue",
232        help="output nodata value (send null for default datatype specific, see `from osgeo_utils.gdal_calc import DefaultNDVLookup`)",
233        type=float,
234        nargs="?",
235        default=-9999,
236    )
237    parser.add_argument(
238        "-f",
239        "--format",
240        help="Output format (send null to infer from file, see `from osgeo_utils.auxiliary.util import GetOutputDriverFor`)",
241        type=str,
242        default="GTiff",
243    )
244
245    parser.add_argument(
246        "-t", "--type", help="Output datatype", type=str, default="Float32", choices=list(map(str, GDALDataTypeNames))
247    )
248    parser.add_argument(
249        "-p",
250        "--projwin",
251        nargs=4,
252        type=float,
253        metavar=("min_x", "max_y", "max_x", "min_y"),
254        help="An optional list of 4 coordinates defining the projection window, if not provided the whole raster is calculated",
255    )
256    parser.add_argument(
257        "-r",
258        "--return_dataset",
259        help="Return dataset (for scripting -additional keyword arguments are passed to gdal_calc.Calc) instead of return code",
260        action="store_true",
261    )
262    args = parser.parse_args(argv)
263    args.projwin = tuple(args.projwin) if args.projwin else None
264    if not args.infile.is_file():
265        parser.error("Input file does not exist")
266    if args.format is None:
267        args.format = GetOutputDriverFor(args.outfile)
268    return args
269
270
271def main(argv=None):
272    """
273    <pre><code>
274    minmax: (A-minimum)/(maximum - minimum)
275    maxmin: (A-maximum)/(minimum - maximum)
276    stepup: 0*(A<threshold)+1*(A>=threshold)
277    stepdown: 1*(A<threshold)+0*(A>=threshold)
278    bipiecewiselinear: (A-a)/(b-a) then "0*(A<0)+1*(A>1)"
279    bipiecewiselinear_percent: (A-a*r)/(b*r-a*r) then "0*(A<0)+1*(A>1)"
280    stepup_percent: 0*(A<threshold*r)+1*(A>=threshold*r)
281    stepdown_percent: 1*(A<threshold*r)+0*(A>=threshold*r)
282
283    r : relative delta : data.max() - data.min() / 100
284    </code></pre>
285    """
286    if argv is sys.argv:
287        argv = sys.argv[1:]
288    args = arg_parser(argv)
289
290    print(f"{args=}")
291
292    if args.method == "minmax":
293        del args.params
294        func = lambda minimum, maximum: f"(A-{minimum})/({maximum} - {minimum})"
295        ds = calc(func, **vars(args))
296    elif args.method == "maxmin":
297        del args.params
298        func = lambda minimum, maximum: f"(A-{maximum})/({minimum} - {maximum})"
299        ds = calc(func, **vars(args))
300    elif args.method == "stepup":
301        threshold = args.params[0]
302        del args.params
303        func = lambda threshold: f"0*(A<{threshold})+1*(A>={threshold})"
304        ds = calc(func, **vars(args), threshold=threshold)
305    elif args.method == "stepdown":
306        threshold = args.params[0]
307        del args.params
308        func = lambda threshold: f"1*(A<{threshold})+0*(A>={threshold})"
309        ds = calc(func, **vars(args), threshold=threshold)
310    elif args.method == "bipiecewiselinear":
311        a = args.params[0]
312        b = args.params[1]
313        del args.params
314
315        keep = args.outfile
316        args.outfile = Path(NamedTemporaryFile(suffix=".tif", delete=False).name)
317        print(f"{args=}")
318
319        func = lambda a, b: f"(A-{a})/({b}-{a})"
320        ds = calc(func, **vars(args), a=a, b=b)
321
322        args.infile = args.outfile
323        args.outfile = keep
324        print(f"{args=}")
325
326        func = lambda: "0*(A<0)+1*(A>1)"
327        ds = calc(func, **vars(args))
328    elif args.method == "bipiecewiselinear_percent":
329        """
330        rela_delta = data.max() - data.min() / 100
331        real_a = rela_delta * a
332        real_b = rela_delta * b
333        data = (data - real_a) / (real_b - real_a)
334        """
335        a = args.params[0]
336        b = args.params[1]
337        del args.params
338
339        keep = args.outfile
340        args.outfile = Path(NamedTemporaryFile(suffix=".tif", delete=False).name)
341        print(f"{args=}")
342
343        func = lambda a, b, r: f"(A-{a*r})/({b*r}-{a*r})"
344        ds = calc(func, **vars(args), a=a, b=b, r=0)
345
346        args.infile = args.outfile
347        args.outfile = keep
348        print(f"{args=}")
349
350        func = lambda: "0*(A<0)+1*(A>1)"
351        ds = calc(func, **vars(args))
352    elif args.method == "stepup_percent":
353        threshold = args.params[0]
354        del args.params
355        func = lambda threshold, r: f"0*(A<{threshold*r})+1*(A>={threshold*r})"
356        ds = calc(func, **vars(args), threshold=threshold, r=0)
357    elif args.method == "stepdown_percent":
358        threshold = args.params[0]
359        del args.params
360        func = lambda threshold, r: f"1*(A<{threshold*r})+0*(A>={threshold*r})"
361        ds = calc(func, **vars(args), threshold=threshold, r=0)
362    """
363    elif args.method == "":
364        del args.method
365        func = lambda : f""
366        ds = maxmin(func, **vars(args))
367    """
368
369    if args.return_dataset:
370        return ds
371
372    if isinstance(ds, Dataset):
373        return 0
374    return 1
375
376
377def get_file_minmax(filename, force=True):
378    # try:
379    dataset = Open(filename, GA_ReadOnly)
380    # except RuntimeError as e:
381    #     if "not recognized as a supported file format" in str(e):
382    #         raise FileNotFoundError("not recognized as a supported file format " + filename)
383    if dataset is None:
384        raise FileNotFoundError(filename)
385    raster_band = dataset.GetRasterBand(1)
386    rmin = raster_band.GetMinimum()
387    rmax = raster_band.GetMaximum()
388    if not rmin or not rmax or force:
389        (rmin, rmax) = raster_band.ComputeRasterMinMax(True)
390    return rmin, rmax
391
392
393if __name__ == "__main__":
394    sys.exit(main(sys.argv))