buteo.raster.grid

  1"""### Create grids from rasters.  ###
  2
  3Cut rasters to grids. Use vectors or rasters as grids.
  4"""
  5
  6# TODO: Verify this function and create more grid options.
  7
  8# Standard library
  9from uuid import uuid4
 10from typing import Union, Optional, List
 11
 12# External
 13from osgeo import gdal, ogr
 14
 15# Internal
 16from buteo.utils import utils_base, utils_gdal, utils_bbox, utils_path
 17from buteo.raster import core_raster, core_stack
 18from buteo.raster.clip import _raster_clip
 19from buteo.vector import core_vector
 20from buteo.vector.intersect import _vector_intersect
 21from buteo.vector.reproject import _vector_reproject
 22from buteo.vector.metadata import _vector_to_metadata
 23
 24
 25
 26def raster_to_grid(
 27    raster: Union[str, gdal.Dataset],
 28    grid: Union[str, ogr.DataSource],
 29    out_dir: str,
 30    *,
 31    use_field: bool = None,
 32    generate_vrt: bool = True,
 33    overwrite: bool = True,
 34    process_layer: int = 0,
 35    creation_options: Optional[List[str]] = None,
 36    verbose: int = 0,
 37) -> str:
 38    """Clips a raster to a grid. Generates .vrt.
 39
 40    Parameters
 41    ----------
 42    raster : Union[str, gdal.Dataset]
 43        The input raster.
 44
 45    grid : Union[str, ogr.DataSource]
 46        The grid to use.
 47
 48    out_dir : str
 49        The output directory.
 50
 51    use_field : bool, optional
 52        A field to use to name the grid cells, default: None.
 53
 54    generate_vrt : bool, optional
 55        If True, the output raster will be a .vrt, default: True.
 56
 57    overwrite : bool, optional
 58        If True, the output raster will be overwritten, default: True.
 59
 60    process_layer : int, optional
 61        The layer from the grid to process, default: 0.
 62
 63    creation_options : Optional[List[str]], optional
 64        Creation options for the output raster, default: None.
 65
 66    verbose : int, optional
 67        The verbosity level, default: 0.
 68
 69    Returns
 70    -------
 71    str
 72        The file path for the newly created raster.
 73    """
 74    utils_base._type_check(raster, [str, gdal.Dataset], "raster")
 75    utils_base._type_check(grid, [str, ogr.DataSource], "grid")
 76    utils_base._type_check(out_dir, [str], "out_dir")
 77    utils_base._type_check(overwrite, [bool], "overwrite")
 78    utils_base._type_check(process_layer, [int], "process_layer")
 79    utils_base._type_check(creation_options, [[str], None], "creation_options")
 80    utils_base._type_check(verbose, [int], "verbose")
 81
 82    use_grid = core_vector.vector_open(grid)
 83    grid_metadata = _vector_to_metadata(use_grid)
 84    raster_metadata = core_raster._get_basic_metadata_raster(raster)
 85
 86    # Reproject raster if necessary.
 87    if not raster_metadata["projection_osr"].IsSame(grid_metadata["projection_osr"]):
 88        use_grid = _vector_reproject(grid, raster_metadata["projection_osr"])
 89        grid_metadata = _vector_to_metadata(use_grid)
 90
 91        if not isinstance(grid_metadata, dict):
 92            raise RuntimeError("Error while parsing metadata.")
 93
 94    # Only use the polygons in the grid that intersect the extent of the raster.
 95    use_grid = _vector_intersect(use_grid, raster_metadata["extent_datasource"]())
 96
 97    ref = core_raster._raster_open(raster)
 98    use_grid = core_vector.vector_open(use_grid)
 99
100    layer = use_grid.GetLayer(process_layer)
101    feature_count = layer.GetFeatureCount()
102    raster_extent = raster_metadata["bbox"]
103    filetype = utils_path._get_ext_from_path(raster)
104    name = raster_metadata["name"]
105    geom_type = grid_metadata["layers"][process_layer]["geom_type_ogr"]
106
107    if use_field is not None:
108        if use_field not in grid_metadata["layers"][process_layer]["field_names"]:
109            names = grid_metadata["layers"][process_layer]["field_names"]
110            raise ValueError(
111                f"Requested field not found. Fields available are: {names}"
112            )
113
114    generated = []
115
116    # For the sake of good reporting - lets first establish how many features intersect
117    # the raster.
118
119    if verbose:
120        print("Finding intersections.")
121
122    intersections = 0
123    for _ in range(feature_count):
124        feature = layer.GetNextFeature()
125        geom = feature.GetGeometryRef()
126
127        if not utils_bbox._check_bboxes_intersect(raster_extent, geom.GetEnvelope()):
128            continue
129
130        intersections += 1
131
132    layer.ResetReading()
133
134    if verbose:
135        print(f"Found {intersections} intersections.")
136
137    if intersections == 0:
138        print("Warning: Found 0 intersections. Returning empty list.")
139        return ([], None)
140
141    driver = ogr.GetDriverByName("GPKG")
142
143    clipped = 0
144    for _ in range(feature_count):
145
146        feature = layer.GetNextFeature()
147        geom = feature.GetGeometryRef()
148
149        if not utils_bbox._check_bboxes_intersect(raster_extent, geom.GetEnvelope()):
150            continue
151
152        if verbose == 1:
153            utils_base.progress(clipped, intersections - 1, "clip_grid")
154
155        fid = feature.GetFID()
156
157        test_ds_path = f"/vsimem/grid_{uuid4().int}.gpkg"
158        test_ds = driver.CreateDataSource(test_ds_path)
159        test_ds_lyr = test_ds.CreateLayer(
160            "mem_layer_grid",
161            geom_type=geom_type,
162            srs=raster_metadata["projection_osr"],
163        )
164        test_ds_lyr.CreateFeature(feature.Clone())
165        test_ds.FlushCache()
166
167        out_name = None
168
169        if use_field is not None:
170            out_name = f"{out_dir}{feature.GetField(use_field)}{filetype}"
171        else:
172            out_name = f"{out_dir}{name}_{fid}{filetype}"
173
174        _raster_clip(
175            ref,
176            test_ds_path,
177            out_path=out_name,
178            adjust_bbox=True,
179            crop_to_geom=True,
180            all_touch=False,
181            suffix="",
182            prefix="",
183            creation_options=utils_gdal._get_default_creation_options(creation_options),
184            verbose=0,
185        )
186
187        generated.append(out_name)
188        clipped += 1
189
190    if generate_vrt:
191        vrt_name = f"{out_dir}{name}.vrt"
192
193        core_stack.raster_stack_vrt_list(generated, vrt_name, separate=False)
194
195        return (generated, vrt_name)
196
197    return generated