fire2a.raster.gdal_calc_sum

  1#!python
  2"""
  3<pre><code>
  4<!-- BEGIN_ARGPARSE_DOCSTRING -->
  5usage: gdal_calc_sum.py [-h] [-o OUTFILE] [-w [WEIGHTS ...]] [-f FORMAT]
  6                        [-t {Byte,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64,CInt16,CInt32,CFloat32,CFloat64}]
  7                        [-p min_x max_y max_x min_y] [-n [NODATAVALUE]] [-r]
  8                        infiles [infiles ...]
  9
 10Raster(s) (weighted) summation utility, wrapping osgeo_utils.gdal_calc for sum(weights*rasters). Run `gdal_calc.py --help` for more
 11information.
 12
 13positional arguments:
 14  infiles               List of rasters to sum up to 52
 15
 16options:
 17  -h, --help            show this help message and exit
 18  -o OUTFILE, --outfile OUTFILE
 19                        Output file (default: outfile.tif)
 20  -w [WEIGHTS ...], --weights [WEIGHTS ...]
 21                        An optional list of weights to ponder the summation (else 1's) (default: None)
 22  -f FORMAT, --format FORMAT
 23                        Output format (default: GTiff)
 24  -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}
 25                        Output datatype (default: Float32)
 26  -p min_x max_y max_x min_y, --projwin min_x max_y max_x min_y
 27                        An optional list of 4 coordinates defining the projection window, if not provided the 1st raster projwin is
 28                        used (default: None)
 29  -n [NODATAVALUE], --NoDataValue [NODATAVALUE]
 30                        output nodata value (send empty for default datatype specific, see `from osgeo_utils.gdal_calc import
 31                        DefaultNDVLookup`) (default: -9999)
 32  -r, --return_dataset  Return dataset (for scripting -additional keyword arguments are passed to gdal_calc.Calc) instead of return
 33                        code (default: False)
 34
 35documentation at https://fire2a.github.io/fire2a-lib/fire2a/raster/gdal_calc_sum.html
 36<!-- END_ARGPARSE_DOCSTRING -->
 37
 38Sample script usage:
 39    from fire2a.raster.gdal_calc_sum import main
 40    ds = main(["-r", ... other keyword arguments are passed to gdal_calc.Calc
 41
 42function osgeo_utils.gdal_calc.Calc(
 43    calc: Union[str, Sequence[str]],
 44    outfile: Union[str, os.PathLike, NoneType] = None,
 45    NoDataValue: Optional[numbers.Number] = None,
 46    type: Union[int, str, NoneType] = None,
 47    format: Optional[str] = None,
 48    creation_options: Optional[Sequence[str]] = None,
 49    allBands: str = '',
 50    overwrite: bool = False,
 51    hideNoData: bool = False,
 52    projectionCheck: bool = False,
 53    color_table: Union[osgeo.gdal.ColorTable,
 54    ForwardRef('ColorPalette'), str, os.PathLike, Sequence[str], NoneType] = None,
 55    extent: Optional[osgeo_utils.auxiliary.extent_util.Extent] = None,
 56    projwin: Union[Tuple, osgeo_utils.auxiliary.rectangle.GeoRectangle, NoneType] = None,
 57    user_namespace: Optional[Dict] = None,
 58    debug: bool = False,
 59    quiet: bool = False, **infile_files)
 60
 61/usr/bin/gdal_calc.py
 62/usr/lib/python3/dist-packages/osgeo_utils/gdal_calc.py
 63</code></pre>
 64"""
 65import sys
 66
 67if sys.version_info < (3, 10):
 68    raise ImportError("This module requires Python 3.10 or higher")
 69
 70import string
 71from pathlib import Path
 72
 73from osgeo.gdal import Dataset
 74
 75try:
 76    from osgeo_utils.auxiliary.util import GetOutputDriverFor
 77    from osgeo_utils.gdal_calc import Calc, GDALDataTypeNames
 78except ImportError:
 79    from osgeo.utils.gdal_calc import Calc, GetOutputDriverFor
 80
 81    GDALDataTypeNames = (
 82        "Byte",
 83        "UInt16",
 84        "Int16",
 85        "UInt32",
 86        "Int32",
 87        "UInt64",
 88        "Int64",
 89        "Float32",
 90        "Float64",
 91        "CInt16",
 92        "CInt32",
 93        "CFloat32",
 94        "CFloat64",
 95    )
 96
 97
 98def calc(
 99    outfile="outfile.tif",
100    infiles=["infile.tif", "infile2.tif"],
101    weights=None,
102    NoDataValue=None,
103    overwrite=True,
104    type="Float32",
105    format="GTiff",
106    projwin=None,
107    **kwargs,
108) -> Dataset:
109    """This is the wrapper function for the gdal_calc.Calc utility.
110
111    Creates the string symbolizing a weighted sum of rasters.
112
113    All extra keyword arguments are passed to the gdal_calc.Calc function."""
114    for i, infile in enumerate(infiles):
115        if isinstance(infile, Path):
116            infiles[i] = str(infile)
117    if isinstance(outfile, Path):
118        outfile = str(outfile)
119    # if not projwin:
120    #     info = locals().get("info", read_raster(infiles[0], data=False)[1])
121    #     projwin = get_projwin(info["Transform"], info["RasterXSize"], info["RasterYSize"])
122    #     print(f"{projwin=}")
123
124    letter_file = {}
125    letter_calc = ""
126
127    AlphaList = list(string.ascii_letters)
128    if not weights:
129        for alpha, afile in zip(AlphaList, infiles):
130            letter_file[alpha] = afile
131            letter_calc += alpha + "+"
132    else:
133        for alpha, afile, weight in zip(AlphaList, infiles, weights):
134            letter_file[alpha] = afile
135            letter_calc += str(weight) + "*" + alpha + "+"
136
137    letter_calc = letter_calc[:-1]
138
139    dataset = Calc(
140        calc=letter_calc,
141        outfile=outfile,
142        NoDataValue=NoDataValue,
143        overwrite=overwrite,
144        type=type,
145        format=format,
146        projwin=projwin,
147        **kwargs,
148        **letter_file,
149    )
150    dataset.FlushCache()
151
152    return dataset
153
154
155def arg_parser(argv=None):
156    """Parse arguments list"""
157    from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
158
159    parser = ArgumentParser(
160        description="Raster(s) (weighted) summation utility, wrapping osgeo_utils.gdal_calc for sum(weights*rasters). Run `gdal_calc.py --help` for more information.",
161        formatter_class=ArgumentDefaultsHelpFormatter,
162        epilog="documentation at https://fire2a.github.io/fire2a-lib/fire2a/raster/gdal_calc_sum.html",
163    )
164    parser.add_argument(
165        "infiles",
166        nargs="+",
167        type=Path,
168        help="List of rasters to sum up to 52",
169    )
170    parser.add_argument("-o", "--outfile", help="Output file", type=Path, default="outfile.tif")
171    parser.add_argument(
172        "-w",
173        "--weights",
174        nargs="*",
175        type=float,
176        help="An optional list of weights to ponder the summation (else 1's)",
177    )
178    parser.add_argument("-f", "--format", help="Output format", type=str, default="GTiff")
179    parser.add_argument(
180        "-t", "--type", help="Output datatype", type=str, default="Float32", choices=list(map(str, GDALDataTypeNames))
181    )
182    parser.add_argument(
183        "-p",
184        "--projwin",
185        nargs=4,
186        type=float,
187        metavar=("min_x", "max_y", "max_x", "min_y"),
188        help="An optional list of 4 coordinates defining the projection window, if not provided the 1st raster projwin is used",
189    )
190    parser.add_argument(
191        "-n",
192        "--NoDataValue",
193        help="output nodata value (send empty for default datatype specific, see `from osgeo_utils.gdal_calc import DefaultNDVLookup`)",
194        type=float,
195        nargs="?",
196        default=-9999,
197    )
198    parser.add_argument(
199        "-r",
200        "--return_dataset",
201        help="Return dataset (for scripting -additional keyword arguments are passed to gdal_calc.Calc) instead of return code",
202        action="store_true",
203    )
204    args = parser.parse_args(argv)
205    args.projwin = tuple(args.projwin) if args.projwin else None
206    if len(args.infiles) > 52:
207        parser.error("Number of input rasters must be less than 53")
208    for infile in args.infiles:
209        if not infile.exists():
210            parser.error(f"Input raster {infile} does not exist")
211    if args.weights:
212        if len(args.weights) != len(args.infiles):
213            parser.error("Number of weights must match the number of input rasters")
214    if args.format is None:
215        args.format = GetOutputDriverFor(args.outfile)
216    return args
217
218
219def main(argv=None):
220    """
221    <pre><code>
222    All arguments passed as:
223    ds = calc(**vars(args))
224    </code></pre>
225    """
226    if argv is sys.argv:
227        argv = sys.argv[1:]
228    args = arg_parser(argv)
229
230    print(f"{args=}")
231
232    ds = calc(**vars(args))
233
234    if args.return_dataset:
235        return ds
236
237    if isinstance(ds, Dataset):
238        return 0
239    return 1
240
241
242if __name__ == "__main__":
243    sys.exit(main(sys.argv))