|
| 1 | +from ..abc_loader import LoaderBase |
| 2 | +import rasterio |
| 3 | +from typing import Any |
| 4 | +import numpy as np |
| 5 | +import geopandas as gpd |
| 6 | +from shapely.geometry import Point |
| 7 | +from shapely.geometry import Polygon |
| 8 | +from rasterio.transform import xy |
| 9 | +from pyproj import CRS, Transformer |
| 10 | + |
| 11 | +class RasterLoader: |
| 12 | + """ |
| 13 | + Loader for raster files (GeoTIFF, PNG+world file, etc.) with block-wise downsampling (average pooling) and polygons as geometry. |
| 14 | + Returns a GeoDataFrame where each row corresponds to an aggregated pixel (block). |
| 15 | + |
| 16 | + It allows fast preview of raster properties, pixel-wise spatialization, and direct integration with the UrbanMapper factory. |
| 17 | +
|
| 18 | + Attributes: |
| 19 | + file_path (str): Path to the raster file to load. |
| 20 | + gdf (geopandas.GeoDataFrame): GeoDataFrame where each row is a pixel (with geometry, area, coordinates, and value). |
| 21 | + meta (dict): Raster metadata (dimensions, CRS, etc.). |
| 22 | + bounds (tuple): Geographic extent of the raster as (left, bottom, right, top). |
| 23 | + block_size (int): Size of the blocks for downsampling (default is 10, meaning 10x10 pixels). |
| 24 | +
|
| 25 | + Example: |
| 26 | + >>> rst_loader = ( |
| 27 | + mapper |
| 28 | + .loader # From the loader module |
| 29 | + .from_file("file_path.tif") # To update with your own path |
| 30 | + ) |
| 31 | + >>> gdf = rst_loader.load() # Load the data and return data |
| 32 | + >>> gdf |
| 33 | + >>> meta = rst_loader._instance.meta |
| 34 | + >>> bounds = rst_loader._instance.bounds |
| 35 | + |
| 36 | + |
| 37 | + """ |
| 38 | + |
| 39 | + def __init__(self, file_path, block_size=10, **kwargs): # block_size est le facteur de downsampling (4x4 par défaut) |
| 40 | + self.file_path = file_path |
| 41 | + self.block_size = block_size |
| 42 | + self.meta = None |
| 43 | + self.bounds = None |
| 44 | + |
| 45 | + def _downsample_band(self, band): |
| 46 | + """ |
| 47 | + Downsamples the raster band by averaging non-overlapping blocks of pixels.Effectue le downsampling par blocs et calcule la moyenne pour chaque bloc non-recouvrant. |
| 48 | + """ |
| 49 | + h, w = band.shape |
| 50 | + bs = self.block_size |
| 51 | + |
| 52 | + # Découpe l'image aux dimensions multiples du block_size pour éviter les bords incomplets |
| 53 | + h_ds = h // bs |
| 54 | + w_ds = w // bs |
| 55 | + band_cropped = band[:h_ds * bs, :w_ds * bs] |
| 56 | + |
| 57 | + # Remodèle pour avoir un 4D (h_blocks, block_size, w_blocks, block_size) puis moyenne par bloc |
| 58 | + band_blocks = band_cropped.reshape(h_ds, bs, w_ds, bs) |
| 59 | + band_ds = band_blocks.mean(axis=(1, 3)) # moyenne sur les axes blocs internes |
| 60 | + |
| 61 | + return band_ds |
| 62 | + |
| 63 | + def _load_data_from_file(self) -> gpd.GeoDataFrame: |
| 64 | + """ |
| 65 | + Loads raster data and returns a GeoDataFrame where each row represents an aggregated pixel (bloc) by downsampling (with geometry, area, coordinates, and value). |
| 66 | + NoData pixels are included with value set to None. |
| 67 | +
|
| 68 | + Returns : |
| 69 | + ------- |
| 70 | + gpd.GeoDataFrame |
| 71 | + A GeoDataFrame with columns: pixel_id, row, col, area, latitude, longitude, value, geometry. |
| 72 | + Raises: |
| 73 | + RuntimeError: If there is an error while loading the raster file. |
| 74 | + NB : the loader doesn't return metadata and bounds, but they are stored in the instance attributes (cf docstring example). |
| 75 | + """ |
| 76 | + try: |
| 77 | + with rasterio.open(self.file_path) as src: |
| 78 | + band = src.read(1) |
| 79 | + transform = src.transform |
| 80 | + crs = src.crs |
| 81 | + nodata = src.nodata |
| 82 | + |
| 83 | + self.meta = src.meta |
| 84 | + self.bounds = src.bounds |
| 85 | + |
| 86 | + # Handle NoData: |
| 87 | + if nodata is not None: |
| 88 | + band = np.where(band == nodata, np.nan, band) |
| 89 | + |
| 90 | + # Downsampling |
| 91 | + band_ds = self._downsample_band(band) |
| 92 | + h_ds, w_ds = band_ds.shape |
| 93 | + |
| 94 | + # Generate indices for the downsampled raster |
| 95 | + rows, cols = np.indices((h_ds, w_ds)) |
| 96 | + |
| 97 | + # Calcul coordinates of the center of each block |
| 98 | + bs = self.block_size |
| 99 | + center_rows = rows * bs + bs // 2 |
| 100 | + center_cols = cols * bs + bs // 2 |
| 101 | + |
| 102 | + # Transform raster to world coordinates for block centers |
| 103 | + xs, ys = rasterio.transform.xy(transform, center_rows, center_cols) |
| 104 | + xs = np.array(xs).flatten() |
| 105 | + ys = np.array(ys).flatten() |
| 106 | + |
| 107 | + # Flatten values for the downsampled band |
| 108 | + values = band_ds.flatten() |
| 109 | + |
| 110 | + # Geometry creation : Polygon for each aggregated pixel |
| 111 | + # Each pixel is represented as a polygon with 4 corners |
| 112 | + geoms = [] |
| 113 | + for r, c in zip(center_rows.flatten(), center_cols.flatten()): |
| 114 | + # (r, c) = ligne et colonne du bloc (dans la grille agrégée) |
| 115 | + min_row = r - bs // 2 |
| 116 | + min_col = c - bs // 2 |
| 117 | + max_row = min_row + bs |
| 118 | + max_col = min_col + bs |
| 119 | + |
| 120 | + # Collect coordinates of the 4 corners of the aggregated pixel : |
| 121 | + corners = [ |
| 122 | + rasterio.transform.xy(transform, min_row, min_col, offset='ul'), # haut-gauche |
| 123 | + rasterio.transform.xy(transform, min_row, max_col, offset='ur'), # haut-droit |
| 124 | + rasterio.transform.xy(transform, max_row, max_col, offset='lr'), # bas-droit |
| 125 | + rasterio.transform.xy(transform, max_row, min_col, offset='ll') # bas-gauche |
| 126 | + ] |
| 127 | + poly = Polygon(corners) |
| 128 | + geoms.append(poly) |
| 129 | + |
| 130 | + |
| 131 | + # Latitude/longitude per transformation from CRS to WGS84 |
| 132 | + transformer = Transformer.from_crs(crs, 4326, always_xy=True) |
| 133 | + lon, lat = transformer.transform(xs, ys) |
| 134 | + |
| 135 | + # Area for each block (in projected CRS units) — area of the block |
| 136 | + if CRS.from_user_input(crs).is_projected: |
| 137 | + pixel_width = abs(transform.a) |
| 138 | + pixel_height = abs(transform.e) |
| 139 | + area = (pixel_width * pixel_height) * (bs * bs) |
| 140 | + areas = np.full(values.shape, area) |
| 141 | + else: |
| 142 | + areas = [None] * values.size # Hors CRS projeté, zone complexe : à raffiner si utile |
| 143 | + |
| 144 | + |
| 145 | + # Create GeoDataFrame with pixel_id, row, col, area, value, latitude, longitude, and geometry |
| 146 | + gdf = gpd.GeoDataFrame({ |
| 147 | + 'pixel_id': np.arange(len(values)), |
| 148 | + 'row': rows.flatten(), |
| 149 | + 'col': cols.flatten(), |
| 150 | + 'area': areas, |
| 151 | + 'value': values, |
| 152 | + 'latitude': lat, |
| 153 | + 'longitude': lon, |
| 154 | + 'geometry': geoms |
| 155 | + }, crs=crs) |
| 156 | + |
| 157 | + return gdf |
| 158 | + |
| 159 | + except Exception as e: |
| 160 | + raise RuntimeError(f"Error while loading downsampled raster: {e}") |
| 161 | + |
| 162 | + |
| 163 | + def preview(self, format: str = "ascii") -> Any: |
| 164 | + """ |
| 165 | + Generates a preview of the loaded raster information. |
| 166 | +
|
| 167 | + Args: |
| 168 | + format (str): Output format ("ascii" for text display, "json" for dictionary). |
| 169 | +
|
| 170 | + Returns: |
| 171 | + str or dict: A summary of the raster properties. |
| 172 | +
|
| 173 | + Raises: |
| 174 | + ValueError: If the requested format is not supported. |
| 175 | +
|
| 176 | + """ |
| 177 | + # If metadata is not loaded, try to load it |
| 178 | + if self.meta is None: |
| 179 | + try: |
| 180 | + with rasterio.open(self.file_path) as src: |
| 181 | + self.meta = src.meta |
| 182 | + except Exception as e: |
| 183 | + return f"Unable to open raster: {e}" |
| 184 | + |
| 185 | + # Get main raster information |
| 186 | + shape = ( |
| 187 | + self.meta.get("count", "?"), |
| 188 | + self.meta.get("height", "?"), |
| 189 | + self.meta.get("width", "?") |
| 190 | + ) |
| 191 | + dtype = self.meta.get("dtype", "?") |
| 192 | + crs = self.meta.get("crs", "?") |
| 193 | + |
| 194 | + # Return preview according to requested format |
| 195 | + if format == "ascii": |
| 196 | + return ( |
| 197 | + f"Loader: RasterLoader\n" |
| 198 | + f" File: {self.file_path}\n" |
| 199 | + f" Dimensions (bands, height, width): {shape}\n" |
| 200 | + f" Data type: {dtype}\n" |
| 201 | + f" CRS: {crs}" |
| 202 | + ) |
| 203 | + elif format == "json": |
| 204 | + return { |
| 205 | + "loader": "RasterLoader", |
| 206 | + "file": self.file_path, |
| 207 | + "shape": shape, |
| 208 | + "dtype": str(dtype), |
| 209 | + "crs": str(crs) |
| 210 | + } |
| 211 | + else: |
| 212 | + raise ValueError(f"Unsupported format: {format}") |
| 213 | + |
| 214 | + |
| 215 | + |
| 216 | + |
0 commit comments