fire2a.raster.raster

  1#!python3
  2"""👋🌎 🌲🔥
  3Here are some raster utilities:
  4
  5Sanity check your system's raster format support:
  6
  7`$ python -c "from osgeo import gdal;print('\n'.join(sorted([gdal.GetDriver(i).GetDescription() for i in range(gdal.GetDriverCount())])))"`
  8"""
  9__author__ = "Fernando Badilla"
 10__revision__ = "$Format:%H$"
 11
 12import logging
 13import warnings
 14from pathlib import Path
 15from typing import Any, Dict, List, Optional, Tuple, Union
 16
 17import numpy as np
 18from osgeo import gdal, ogr
 19from qgis.core import QgsRasterLayer, QgsRectangle
 20
 21from fire2a.utils import fprint, qgis2numpy_dtype
 22
 23logger = logging.getLogger(__name__)
 24"""@private"""
 25
 26
 27def id2xy(idx: int, w: int, h: int) -> tuple[int, int]:
 28    """Transform a pixel or cell index, into x,y coordinates.
 29
 30    In GIS, the origin is at the top-left corner, read from left to right, top to bottom.  
 31    If your're used to matplotlib, the y-axis is inverted.  
 32    Also as numpy array, the index of the pixel is [y, x].
 33
 34    Args:
 35
 36        param idx: index of the pixel or cell (0,..,w*h-1)  
 37        param w: width of the image or grid  
 38        param h: height of the image or grid (not really used!)
 39
 40    Returns:
 41
 42        tuple: (x, y) coordinates of the pixel or cell  
 43    """  # fmt: skip
 44    return idx % w, idx // w
 45
 46
 47def xy2id(x: int, y: int, w: int) -> int:
 48    """Transform a x,y coordinates into a pixel or cell index.
 49
 50    In GIS, the origin is at the top-left corner, read from left to right, top to bottom.
 51    If your're used to matplotlib, the y-axis is inverted.
 52    Also as numpy array, the index of the pixel is [y, x].
 53
 54    Args:
 55
 56        param x: width or horizontal coordinate of the pixel or cell
 57        param y: height or vertical coordinate of the pixel or cell
 58        param w: width of the image or grid
 59
 60    Returns:
 61
 62        int: index of the pixel or cell (0,..,w*h-1)
 63    """
 64    return y * w + x
 65
 66
 67def read_raster_band(filename: str, band: int = 1) -> tuple[np.ndarray, int, int]:
 68    """Read a raster file and return the data as a numpy array, along width and height.
 69
 70    Args:
 71
 72        param filename: name of the raster file
 73        param band: band number to read (default 1)
 74
 75    Returns:
 76
 77        tuple: (data, width, height)
 78
 79    Raises:
 80
 81        FileNotFoundError: if the file is not found
 82    """
 83    dataset = gdal.Open(filename, gdal.GA_ReadOnly)
 84    if dataset is None:
 85        raise FileNotFoundError(filename)
 86    return dataset.GetRasterBand(band).ReadAsArray(), dataset.RasterXSize, dataset.RasterYSize
 87
 88
 89def read_raster(
 90    filename: str, band: int = 1, data: bool = True, info: bool = True
 91) -> tuple[Optional[np.ndarray], Optional[dict]]:
 92    """Reads a raster file gets the data as a numpy array along useful raster info: transform, projection, raster count, raster width, raster height.
 93
 94    Args:
 95
 96        param filename: name of the raster file
 97        param band: band number to read (default 1)
 98        param data: if True, return the data as a numpy array (default True)
 99        param info: if True, return the raster info (default True)
100
101    Return tuple: (data, info)
102
103        data: numpy 2d array with the raster data
104        info: dictionary with keys:
105            - Transform: geotransform parameters
106            - Projection: projection string
107            - RasterCount: number of bands
108            - RasterXSize: width of the raster
109            - RasterYSize: height of the raster
110            - DataType: data type of the raster
111            - NoDataValue: no data value of the raster
112            - Minimum: minimum value of the raster
113            - Maximum: maximum value of the raster
114
115    Raises:
116
117        FileNotFoundError: if the file is not found
118    """
119    dataset = gdal.Open(filename, gdal.GA_ReadOnly)
120    if dataset is None:
121        raise FileNotFoundError(filename)
122    raster_band = dataset.GetRasterBand(band)
123    data_output = raster_band.ReadAsArray() if data else None
124
125    if info:
126        rmin = raster_band.GetMinimum()
127        rmax = raster_band.GetMaximum()
128        if not rmin or not rmax:
129            (rmin, rmax) = raster_band.ComputeRasterMinMax(True)
130
131    info_output = (
132        {
133            "Transform": dataset.GetGeoTransform(),
134            "Projection": dataset.GetProjection(),
135            "RasterCount": dataset.RasterCount,
136            "RasterXSize": dataset.RasterXSize,
137            "RasterYSize": dataset.RasterYSize,
138            "DataType": gdal.GetDataTypeName(raster_band.DataType),
139            "NoDataValue": raster_band.GetNoDataValue(),
140            "Minimum": rmin,
141            "Maximum": rmax,
142        }
143        if info
144        else None
145    )
146    return data_output, info_output
147
148
149def get_geotransform(raster_filename: str) -> tuple[float, float, float, float, float, float]:
150    """Get geotransform from raster file.
151
152    Args:
153
154        raster_filename (str):
155
156    Returns: tuple: geotransform
157
158        GT[0] x-coordinate of the upper-left corner of the upper-left pixel.
159        GT[1] w-e pixel resolution / pixel width.
160        GT[2] row rotation (typically zero).
161        GT[3] y-coordinate of the upper-left corner of the upper-left pixel.
162        GT[4] column rotation (typically zero).
163        GT[5] n-s pixel resolution / pixel height (negative value for a north-up image).
164
165    reference: https://gdal.org/tutorials/geotransforms_tut.html
166    """
167    dataset = gdal.Open(raster_filename, gdal.GA_ReadOnly)
168    if dataset is None:
169        raise Exception(f"Data set is None, could not open {raster_filename}")
170    return dataset.GetGeoTransform()
171
172
173def transform_coords_to_georef(x_pixel: int, y_line: int, GT: tuple) -> tuple[float, float]:
174    """ Transform pixel coordinates to georeferenced coordinates.
175
176    Args:
177
178        x_pixel (int): x pixel coordinate.
179        y_line (int): y pixel coordinate.
180        GT (tuple): geotransform, see get_geotransform(filename)
181
182    Returns:
183
184        tuple: x_geo, y_geo.
185
186    reference: https://gdal.org/tutorials/geotransforms_tut.html
187    """  # fmt: skip
188    x_geo = GT[0] + x_pixel * GT[1] + y_line * GT[2]
189    y_geo = GT[3] + x_pixel * GT[4] + y_line * GT[5]
190    return x_geo, y_geo
191
192
193def transform_georef_to_coords(x_geo: int, y_geo: int, GT: tuple) -> tuple[float, float]:
194    """Inverse of transform_coords_to_georef.
195
196    Made with symbolic-py package, to solve the system of equations:
197
198        import sympy
199        a, b, c, d, e, f, g, i, j, x, y = sympy.symbols('a, b, c, d, e, f, g, i, j, x, y', real=True)
200        sympy.linsolve([a+i*b+j*c - x,d+i*e+j*f-y],(i,j))
201        {((-a*f + c*d - c*y + f*x)/(b*f - c*e), (a*e - b*d + b*y - e*x)/(b*f - c*e))}
202
203    Args:
204
205        x_geo (int): x georeferenced coordinate.
206        y_geo (int): y georeferenced coordinate.
207        GT (tuple): geotransform, see get_geotransform(filename)
208
209    Returns:
210
211        tuple: x_pixel, y_line.
212
213    Help!
214
215        Implement Raise ValueError Exception ?
216        if x_pixel or y_line are not integer coordinates. by setting a tolerance?
217
218    reference: https://gdal.org/tutorials/geotransforms_tut.html
219    """
220    a, b, c, d, e, f = GT
221    x, y = x_geo, y_geo
222    i, j = (-a * f + c * d - c * y + f * x) / (b * f - c * e), (a * e - b * d + b * y - e * x) / (b * f - c * e)
223    # if i % 1 != 0 or j % 1 != 0:
224    #     raise Exception("Not integer coordinates!")
225    return int(i), int(j)
226
227
228def get_rlayer_info(layer: QgsRasterLayer) -> Dict[str, Any]:
229    """Using QGIS, Get raster layer information
230
231    Args:
232
233        layer (QgsRasterLayer): A raster layer
234
235    Returns: raster info dictionary
236
237        width: Raster width
238        height: Raster height
239        extent: Raster extent (QgsRectangle)
240        crs: Raster CRS
241        cellsize_x: Raster cell size in x
242        cellsize_y: Raster cell size in y
243        nodata: No data value (could be a list)
244        bands: Number of bands
245        file: Raster file path (str)
246    """
247    provider = layer.dataProvider()
248    ndv = []
249    for band in range(1, layer.bandCount() + 1):
250        ndv += [None]
251        if provider.sourceHasNoDataValue(band):
252            ndv[-1] = provider.sourceNoDataValue(band)
253    return {
254        "width": layer.width(),
255        "height": layer.height(),
256        "extent": layer.extent(),
257        "crs": layer.crs(),
258        "cellsize_x": layer.rasterUnitsPerPixelX(),
259        "cellsize_y": layer.rasterUnitsPerPixelY(),
260        "nodata": ndv,
261        "bands": layer.bandCount(),
262        "file": layer.publicSource(),
263    }
264
265
266def get_rlayer_data(layer: QgsRasterLayer) -> np.ndarray:
267    """Using QGIS, Get raster layer data (EVERY BAND) as numpy array; Also returns nodata value, width and height
268
269    The user should check the shape of the data to determine if it is a single band or multiband raster.
270
271    len(data.shape) == 2 for single band, len(data.shape) == 3 for multiband.
272
273    Args:
274
275        layer (QgsRasterLayer): A raster layer
276
277    Returns:
278
279        data (np.array): Raster data as numpy array
280
281    Help! Can a multiband raster have different nodata values and/or data types for each band?
282
283    TODO? make a band list as input
284
285    """
286    provider = layer.dataProvider()
287    if layer.bandCount() == 1:
288        block = provider.block(1, layer.extent(), layer.width(), layer.height())
289        nodata = None
290        if block.hasNoDataValue():
291            nodata = block.noDataValue()
292        np_dtype = qgis2numpy_dtype(provider.dataType(1))
293        data = np.frombuffer(block.data(), dtype=np_dtype).reshape(layer.height(), layer.width())
294        # return data, nodata, np_dtype
295    else:
296        data = []
297        nodata = []
298        np_dtypel = []
299        for i in range(layer.bandCount()):
300            block = provider.block(i + 1, layer.extent(), layer.width(), layer.height())
301            nodata += [None]
302            if block.hasNoDataValue():
303                nodata[-1] = block.noDataValue()
304            np_dtypel += [qgis2numpy_dtype(provider.dataType(i + 1))]
305            data += [np.frombuffer(block.data(), dtype=np_dtypel[-1]).reshape(layer.height(), layer.width())]
306        # would different data types bug this next line?
307        data = np.array(data)
308        # return data, nodata, np_dtypl
309    return data
310
311
312def get_cell_sizeV2(filename: str, band: int = 1) -> tuple[float, float]:
313    """This function is going to be deprecated.
314
315    Get the cell size(s) of a raster.
316    """
317    warnings.warn("This function is going to be deprecated", DeprecationWarning)
318    _, info = read_raster(filename, band=band, data=False, info=True)
319    return info["RasterXSize"], info["RasterYSize"]
320
321
322def get_cell_size(raster: gdal.Dataset) -> tuple[float, float]:
323    """This function is going to be deprecated.
324
325    Get the cell size(s) of a raster.
326
327    Args:
328
329        raster (gdal.Dataset | str): The GDAL dataset or path to the raster.
330
331    Returns:
332
333        float | tuple[float, float]: The cell size(s) as a single float or a tuple (x, y).
334    """
335    warnings.warn("This function is going to be deprecated", DeprecationWarning)
336    if isinstance(raster, str):
337        ds = gdal.Open(raster, gdal.GA_ReadOnly)
338    elif isinstance(raster, gdal.Dataset):
339        ds = raster
340    else:
341        raise ValueError("Invalid input type for raster")
342
343    # Get the affine transformation parameters
344    affine = ds.GetGeoTransform()
345
346    if affine[1] != -affine[5]:
347        # If x and y cell sizes are not equal
348        cell_size = (affine[1], -affine[5])  # Return as a tuple
349    else:
350        cell_size = affine[1]  # Return as a single float
351
352    return cell_size
353
354
355def mask_raster(raster_ds: gdal.Dataset, band: int, polygons: list[ogr.Geometry]) -> np.ndarray:
356    """This function is going to be deprecated.
357
358    Mask a raster with polygons using GDAL.
359
360    Args:
361
362        raster_ds (gdal.Dataset): GDAL dataset of the raster.
363        band (int): Band index of the raster.
364        polygons (list[ogr.Geometry]): List of OGR geometries representing polygons for masking.
365
366    Returns:
367
368        np.array: Masked raster data as a NumPy array.
369    """
370    warnings.warn("This function is going to be deprecated", DeprecationWarning)
371
372    # Get the mask as a NumPy boolean array
373    mask_array = rasterize_polygons(polygons, raster_ds.RasterXSize, raster_ds.RasterYSize)
374
375    # Read the original raster data
376    original_data = band.ReadAsArray()  #  FIXME: wrong type hint : int has no attribute ReadAsArray
377
378    # Apply the mask
379    masked_data = np.where(mask_array, original_data, np.nan)
380
381    return masked_data
382
383
384def rasterize_polygons(polygons: list[ogr.Geometry], width: int, height: int) -> np.ndarray:
385    """Rasterize polygons to a boolean array.
386
387    Args:
388
389        polygons (list[ogr.Geometry]): List of OGR geometries representing polygons for rasterization.
390        geo_transform (tuple): GeoTransform parameters for the raster.
391        width (int): Width of the raster.
392        height (int): Height of the raster.
393
394    Returns:
395
396        mask_array (np.array): Rasterized mask as a boolean array.
397    """
398
399    mask_array = np.zeros((height, width), dtype=bool)
400
401    # Create an in-memory layer to hold the polygons
402    mem_driver = ogr.GetDriverByName("Memory")
403    mem_ds = mem_driver.CreateDataSource("memData")
404    mem_layer = mem_ds.CreateLayer("memLayer", srs=None, geom_type=ogr.wkbPolygon)
405
406    for geometry in polygons:
407        mem_feature = ogr.Feature(mem_layer.GetLayerDefn())
408        mem_feature.SetGeometry(geometry.Clone())
409        mem_layer.CreateFeature(mem_feature)
410
411    # Rasterize the in-memory layer and update the mask array
412    gdal.RasterizeLayer(mask_array, [1], mem_layer, burn_values=[1])
413
414    mem_ds = None  # Release the in-memory dataset
415
416    return mask_array
417
418
419def stack_rasters(
420    file_list: list[Path], mask_polygon: Union[list[ogr.Geometry], None] = None
421) -> tuple[np.ndarray, list[str]]:
422    """This function is going to be deprecated.
423
424    Stack raster files from a list into a 3D NumPy array.
425
426    Args:
427
428        file_list (list[Path]): List of paths to raster files.
429        mask_polygon (list[ogr.Geometry], optional): List of OGR geometries for masking. Defaults to None.
430
431    Returns:
432
433        np.array: Stacked raster array.
434        list: List of layer names corresponding to the rasters.
435    """
436    warnings.warn("This function is going to be deprecated", DeprecationWarning)
437    array_list = []
438    cell_sizes = set()
439    layer_names = []
440
441    for raster_path in file_list:
442        layer_name = raster_path.stem
443        layer_names.append(layer_name)
444
445        ds = gdal.Open(str(raster_path))
446        if ds is None:
447            raise ValueError(f"Failed to open raster file: {raster_path}")
448
449        band = ds.GetRasterBand(1)
450
451        if mask_polygon:
452            flatten_array = mask_raster(ds, band, mask_polygon)
453        else:
454            flatten_array = band.ReadAsArray()
455
456        array_list.append(flatten_array)
457        cell_sizes.add(get_cell_size(ds))
458
459    assert len(cell_sizes) == 1, f"There are rasters with different cell sizes: {cell_sizes}"
460    stacked_array = np.stack(array_list, axis=0)  #  type: np.array
461    print(stacked_array.shape)
462    return stacked_array, layer_names
463
464
465def write_raster(
466    data,
467    outfile="output.tif",
468    driver_name="GTiff",
469    authid="EPSG:3857",
470    geotransform=(0, 1, 0, 0, 0, -1),
471    nodata: Optional[int] = None,
472    feedback=None,
473    logger=None,  # logger default ?
474) -> bool:
475    """Write a raster file from a numpy array.
476
477    To spatially match another raster, get authid and geotransform using:
478
479        from fire2a.raster import read_raster
480        _,info = read_raster(filename, data=False, info=True).
481        authid = info["Transform"]
482        geotransform = info["Projection"].
483
484    Args:
485
486        data (np.array): numpy array to write as raster
487        outfile (str, optional): output raster filename. Defaults to "output.tif".
488        driver_name (str, optional): GDAL driver name. Defaults to "GTiff".
489        authid (str, optional): EPSG code. Defaults to "EPSG:3857".
490        geotransform (tuple, optional): geotransform parameters. Defaults to (0, 1, 0, 0, 0, 1).
491        feedback (Optional, optional): qgsprocessing.feedback object. Defaults to None.
492        logger ([type], optional): logging.logger object. Defaults to None.
493
494    Returns:
495
496        bool: True if the raster was written successfully, False otherwise.
497    """
498    try:
499        from fire2a.processing_utils import get_output_raster_format
500
501        driver_name = get_output_raster_format(outfile, feedback=feedback)
502    except Exception as e:
503        fprint(
504            f"Couln't get output raster format: {e}, defaulting to GTiff",
505            level="warning",
506            feedback=feedback,
507            logger=logger,
508        )
509        driver_name = "GTiff"
510    H, W = data.shape
511    ds = gdal.GetDriverByName(driver_name).Create(outfile, W, H, 1, gdal.GDT_Float32)
512    ds.SetGeoTransform(geotransform)
513    ds.SetProjection(authid)
514    band = ds.GetRasterBand(1)
515    if 0 != band.WriteArray(data):
516        fprint("WriteArray failed", level="warning", feedback=feedback, logger=logger)
517        return False
518    if nodata and data[data == nodata].size > 0:
519        band.SetNoDataValue(nodata)
520        # TBD : always returns 1?
521        # if 0 != band.SetNoDataValue(nodata):
522        #     fprint("Set NoData failed", level="warning", feedback=feedback, logger=logger)
523        #     return False
524    ds.FlushCache()
525    ds = None
526    return True
527
528
529def get_projwin(
530    transform: Tuple[float, float, float, float, float, float],
531    width: int,
532    height: int,
533) -> Tuple[float, float, float, float]:
534    """Calculate the projwin from the raster transform and size.
535
536    Args:
537
538        transform: geotransform parameters
539        width :  of the raster
540        height : of the raster
541
542    Returns:
543
544        projwin: (min_x, max_y, max_x, min_y)
545
546    Example:
547
548        transform = (325692.3826, 20.0, 0.0, 4569655.6528, 0.0, -20.0)
549        raster_x_size = 658
550        raster_y_size = 597
551        projwin = get_projwin(transform, raster_x_size, raster_y_size)
552    """
553
554    min_x = transform[0]
555    max_x = transform[0] + width * transform[1]
556    max_y = transform[3]
557    min_y = transform[3] + height * transform[5]
558
559    projwin = (min_x, max_y, max_x, min_y)
560    # print(projwin)
561    return projwin
562
563
564def extent_to_projwin(extent: QgsRectangle) -> Tuple[float, float, float, float]:
565    """Transform a QgsRectangle extent to a projwin format. Scrambling the order"""
566    # Extract the coordinates
567    min_x = extent.xMinimum()
568    max_x = extent.xMaximum()
569    min_y = extent.yMinimum()
570    max_y = extent.yMaximum()
571
572    # Convert to projwin format (min_x, max_y, max_x, min_y)
573    projwin = (min_x, max_y, max_x, min_y)
574    return projwin
575
576
577if __name__ == "__main__":
578    file_list = list(Path().cwd().glob("*.asc"))
579    print(file_list)
580    array = stack_rasters(file_list)
581    print(array)