| 1 | # SPDX-License-Identifier: Apache-2.0 |
| 2 | # SPDX-FileCopyrightText: Copyright contributors to the vLLM project |
| 3 | import argparse |
| 4 | import datetime |
| 5 | import os |
| 6 | |
| 7 | import albumentations |
| 8 | import numpy as np |
| 9 | import rasterio |
| 10 | import regex as re |
| 11 | import torch |
| 12 | from einops import rearrange |
| 13 | from terratorch.datamodules import Sen1Floods11NonGeoDataModule |
| 14 | |
| 15 | from vllm import LLM |
| 16 | |
| 17 | torch.set_default_dtype(torch.float16) |
| 18 | |
| 19 | NO_DATA = -9999 |
| 20 | NO_DATA_FLOAT = 0.0001 |
| 21 | OFFSET = 0 |
| 22 | PERCENTILE = 99 |
| 23 | |
| 24 | datamodule_config = { |
| 25 | "bands": ["BLUE", "GREEN", "RED", "NIR_NARROW", "SWIR_1", "SWIR_2"], |
| 26 | "batch_size": 16, |
| 27 | "constant_scale": 0.0001, |
| 28 | "data_root": "/dccstor/geofm-finetuning/datasets/sen1floods11", |
| 29 | "drop_last": True, |
| 30 | "no_data_replace": 0.0, |
| 31 | "no_label_replace": -1, |
| 32 | "num_workers": 8, |
| 33 | "test_transform": [ |
| 34 | albumentations.Resize( |
| 35 | always_apply=False, height=448, interpolation=1, p=1, width=448 |
| 36 | ), |
| 37 | albumentations.pytorch.ToTensorV2( |
| 38 | transpose_mask=False, always_apply=True, p=1.0 |
| 39 | ), |
| 40 | ], |
| 41 | } |
| 42 | |
| 43 | |
| 44 | class PrithviMAE: |
| 45 | def __init__(self, model): |
| 46 | self.model = LLM( |
| 47 | model=model, |
| 48 | skip_tokenizer_init=True, |
| 49 | dtype="float16", |
| 50 | enforce_eager=True, |
| 51 | model_impl="terratorch", |
| 52 | enable_mm_embeds=True, |
| 53 | ) |
| 54 | |
| 55 | def run(self, input_data, location_coords): |
| 56 | # merge the inputs into one data structure |
| 57 | if input_data is not None and input_data.dtype == torch.float32: |
| 58 | input_data = input_data.to(torch.float16) |
| 59 | input_data = input_data[0] |
| 60 | |
| 61 | mm_data = { |
| 62 | "image": { |
| 63 | "pixel_values": input_data, |
| 64 | "location_coords": location_coords, |
| 65 | } |
| 66 | } |
| 67 | |
| 68 | prompt = {"prompt_token_ids": [1], "multi_modal_data": mm_data} |
| 69 | outputs = self.model.encode(prompt, pooling_task="plugin", use_tqdm=False) |
| 70 | |
| 71 | return outputs[0].outputs.data |
| 72 | |
| 73 | |
| 74 | def generate_datamodule(): |
| 75 | datamodule = Sen1Floods11NonGeoDataModule( |
| 76 | data_root=datamodule_config["data_root"], |
| 77 | batch_size=datamodule_config["batch_size"], |
| 78 | num_workers=datamodule_config["num_workers"], |
| 79 | bands=datamodule_config["bands"], |
| 80 | drop_last=datamodule_config["drop_last"], |
| 81 | test_transform=datamodule_config["test_transform"], |
| 82 | ) |
| 83 | |
| 84 | return datamodule |
| 85 | |
| 86 | |
| 87 | def process_channel_group(orig_img, channels): |
| 88 | """Args: |
| 89 | orig_img: torch.Tensor representing original image (reference) |
| 90 | with shape = (bands, H, W). |
| 91 | channels: list of indices representing RGB channels. |
| 92 | |
| 93 | Returns: |
| 94 | torch.Tensor with shape (num_channels, height, width) |
| 95 | for original image |
| 96 | |
| 97 | """ |
| 98 | orig_img = orig_img[channels, ...] |
| 99 | valid_mask = torch.ones_like(orig_img, dtype=torch.bool) |
| 100 | valid_mask[orig_img == NO_DATA_FLOAT] = False |
| 101 | |
| 102 | # Rescale (enhancing contrast) |
| 103 | max_value = max(3000, np.percentile(orig_img[valid_mask], PERCENTILE)) |
| 104 | min_value = OFFSET |
| 105 | |
| 106 | orig_img = torch.clamp((orig_img - min_value) / (max_value - min_value), 0, 1) |
| 107 | |
| 108 | # No data as zeros |
| 109 | orig_img[~valid_mask] = 0 |
| 110 | |
| 111 | return orig_img |
| 112 | |
| 113 | |
| 114 | def read_geotiff(file_path: str): |
| 115 | """Read all bands from *file_path* and return image + meta info. |
| 116 | |
| 117 | Args: |
| 118 | file_path: path to image file. |
| 119 | |
| 120 | Returns: |
| 121 | np.ndarray with shape (bands, height, width) |
| 122 | meta info dict |
| 123 | |
| 124 | """ |
| 125 | with rasterio.open(file_path) as src: |
| 126 | img = src.read() |
| 127 | meta = src.meta |
| 128 | try: |
| 129 | coords = src.lnglat() |
| 130 | except Exception: |
| 131 | # Cannot read coords |
| 132 | coords = None |
| 133 | |
| 134 | return img, meta, coords |
| 135 | |
| 136 | |
| 137 | def save_geotiff(image, output_path: str, meta: dict): |
| 138 | """Save multi-band image in Geotiff file. |
| 139 | |
| 140 | Args: |
| 141 | image: np.ndarray with shape (bands, height, width) |
| 142 | output_path: path where to save the image |
| 143 | meta: dict with meta info. |
| 144 | |
| 145 | """ |
| 146 | with rasterio.open(output_path, "w", **meta) as dest: |
| 147 | for i in range(image.shape[0]): |
| 148 | dest.write(image[i, :, :], i + 1) |
| 149 | |
| 150 | return |
| 151 | |
| 152 | |
| 153 | def _convert_np_uint8(float_image: torch.Tensor): |
| 154 | image = float_image.numpy() * 255.0 |
| 155 | image = image.astype(dtype=np.uint8) |
| 156 | |
| 157 | return image |
| 158 | |
| 159 | |
| 160 | def load_example( |
| 161 | file_paths: list[str], |
| 162 | mean: list[float] = None, |
| 163 | std: list[float] = None, |
| 164 | indices: list[int] | None = None, |
| 165 | ): |
| 166 | """Build an input example by loading images in *file_paths*. |
| 167 | |
| 168 | Args: |
| 169 | file_paths: list of file paths . |
| 170 | mean: list containing mean values for each band in the |
| 171 | images in *file_paths*. |
| 172 | std: list containing std values for each band in the |
| 173 | images in *file_paths*. |
| 174 | indices: optional band indices to select from each image. All bands |
| 175 | are kept when None. |
| 176 | |
| 177 | Returns: |
| 178 | np.array containing created example |
| 179 | list of meta info for each image in *file_paths* |
| 180 | |
| 181 | """ |
| 182 | imgs = [] |
| 183 | metas = [] |
| 184 | temporal_coords = [] |
| 185 | location_coords = [] |
| 186 | |
| 187 | for file in file_paths: |
| 188 | img, meta, coords = read_geotiff(file) |
| 189 | |
| 190 | # Rescaling (don't normalize on nodata) |
| 191 | img = np.moveaxis(img, 0, -1) # channels last for rescaling |
| 192 | if indices is not None: |
| 193 | img = img[..., indices] |
| 194 | if mean is not None and std is not None: |
| 195 | img = np.where(img == NO_DATA, NO_DATA_FLOAT, (img - mean) / std) |
| 196 | |
| 197 | imgs.append(img) |
| 198 | metas.append(meta) |
| 199 | if coords is not None: |
| 200 | location_coords.append(coords) |
| 201 | |
| 202 | try: |
| 203 | match = re.search(r"(\d{7,8}T\d{6})", file) |
| 204 | if match: |
| 205 | year = int(match.group(1)[:4]) |
| 206 | julian_day = match.group(1).split("T")[0][4:] |
| 207 | if len(julian_day) == 3: |
| 208 | julian_day = int(julian_day) |
| 209 | else: |
| 210 | julian_day = ( |
| 211 | datetime.datetime.strptime(julian_day, "%m%d") |
| 212 | .timetuple() |
| 213 | .tm_yday |
| 214 | ) |
| 215 | temporal_coords.append([year, julian_day]) |
| 216 | except Exception as e: |
| 217 | print(f"Could not extract timestamp for {file} ({e})") |
| 218 | |
| 219 | imgs = np.stack(imgs, axis=0) # num_frames, H, W, C |
| 220 | imgs = np.moveaxis(imgs, -1, 0).astype("float32") # C, num_frames, H, W |
| 221 | imgs = np.expand_dims(imgs, axis=0) # add batch di |
| 222 | |
| 223 | return imgs, temporal_coords, location_coords, metas |
| 224 | |
| 225 | |
| 226 | def run_model( |
| 227 | input_data, |
| 228 | temporal_coords, |
| 229 | location_coords, |
| 230 | model, |
| 231 | datamodule, |
| 232 | img_size, |
| 233 | lightning_model=None, |
| 234 | ): |
| 235 | # Reflect pad if not divisible by img_size |
| 236 | original_h, original_w = input_data.shape[-2:] |
| 237 | pad_h = (img_size - (original_h % img_size)) % img_size |
| 238 | pad_w = (img_size - (original_w % img_size)) % img_size |
| 239 | input_data = np.pad( |
| 240 | input_data, ((0, 0), (0, 0), (0, 0), (0, pad_h), (0, pad_w)), mode="reflect" |
| 241 | ) |
| 242 | |
| 243 | # Build sliding window |
| 244 | |
| 245 | batch_size = 1 |
| 246 | # batch = torch.tensor(input_data, device="cpu") |
| 247 | batch = torch.tensor(input_data) |
| 248 | windows = batch.unfold(3, img_size, img_size).unfold(4, img_size, img_size) |
| 249 | h1, w1 = windows.shape[3:5] |
| 250 | windows = rearrange( |
| 251 | windows, "b c t h1 w1 h w -> (b h1 w1) c t h w", h=img_size, w=img_size |
| 252 | ) |
| 253 | |
| 254 | # Split into batches if number of windows > batch_size |
| 255 | num_batches = windows.shape[0] // batch_size if windows.shape[0] > batch_size else 1 |
| 256 | windows = torch.tensor_split(windows, num_batches, dim=0) |
| 257 | |
| 258 | if temporal_coords: |
| 259 | temporal_coords = torch.tensor(temporal_coords).unsqueeze(0) |
| 260 | else: |
| 261 | temporal_coords = None |
| 262 | if location_coords: |
| 263 | location_coords = torch.tensor(location_coords[0]).unsqueeze(0) |
| 264 | else: |
| 265 | location_coords = None |
| 266 | |
| 267 | # Run Prithvi-EO-V2-300M-TL-Sen1Floods11 |
| 268 | pred_imgs = [] |
| 269 | for x in windows: |
| 270 | # Apply standardization |
| 271 | x = datamodule.test_transform(image=x.squeeze().numpy().transpose(1, 2, 0)) |
| 272 | x = datamodule.aug(x)["image"] |
| 273 | |
| 274 | with torch.no_grad(): |
| 275 | pred = model.run(x, location_coords=location_coords) |
| 276 | y_hat = pred.argmax(dim=1) |
| 277 | |
| 278 | y_hat = torch.nn.functional.interpolate( |
| 279 | y_hat.unsqueeze(1).float(), size=img_size, mode="nearest" |
| 280 | ) |
| 281 | |
| 282 | pred_imgs.append(y_hat) |
| 283 | |
| 284 | pred_imgs = torch.concat(pred_imgs, dim=0) |
| 285 | |
| 286 | # Build images from patches |
| 287 | pred_imgs = rearrange( |
| 288 | pred_imgs, |
| 289 | "(b h1 w1) c h w -> b c (h1 h) (w1 w)", |
| 290 | h=img_size, |
| 291 | w=img_size, |
| 292 | b=1, |
| 293 | c=1, |
| 294 | h1=h1, |
| 295 | w1=w1, |
| 296 | ) |
| 297 | |
| 298 | # Cut padded area back to original size |
| 299 | pred_imgs = pred_imgs[..., :original_h, :original_w] |
| 300 | |
| 301 | # Squeeze (batch size 1) |
| 302 | pred_imgs = pred_imgs[0] |
| 303 | |
| 304 | return pred_imgs |
| 305 | |
| 306 | |
| 307 | def main( |
| 308 | data_file: str, |
| 309 | model: str, |
| 310 | output_dir: str, |
| 311 | rgb_outputs: bool, |
| 312 | input_indices: list[int] = None, |
| 313 | ): |
| 314 | os.makedirs(output_dir, exist_ok=True) |
| 315 | |
| 316 | model_obj = PrithviMAE(model=model) |
| 317 | datamodule = generate_datamodule() |
| 318 | img_size = 512 # Size of Sen1Floods11 |
| 319 | |
| 320 | input_data, temporal_coords, location_coords, meta_data = load_example( |
| 321 | file_paths=[data_file], |
| 322 | indices=input_indices, |
| 323 | ) |
| 324 | |
| 325 | meta_data = meta_data[0] # only one image |
| 326 | |
| 327 | if input_data.mean() > 1: |
| 328 | input_data = input_data / 10000 # Convert to range 0-1 |
| 329 | |
| 330 | channels = [ |
| 331 | datamodule_config["bands"].index(b) for b in ["RED", "GREEN", "BLUE"] |
| 332 | ] # BGR -> RGB |
| 333 | |
| 334 | pred = run_model( |
| 335 | input_data, temporal_coords, location_coords, model_obj, datamodule, img_size |
| 336 | ) |
| 337 | # Save pred |
| 338 | meta_data.update(count=1, dtype="uint8", compress="lzw", nodata=0) |
| 339 | pred_file = os.path.join( |
| 340 | output_dir, f"pred_{os.path.splitext(os.path.basename(data_file))[0]}.tiff" |
| 341 | ) |
| 342 | save_geotiff(_convert_np_uint8(pred), pred_file, meta_data) |
| 343 | |
| 344 | # Save image + pred |
| 345 | meta_data.update(count=3, dtype="uint8", compress="lzw", nodata=0) |
| 346 | |
| 347 | if input_data.mean() < 1: |
| 348 | input_data = input_data * 10000 # Scale to 0-10000 |
| 349 | |
| 350 | rgb_orig = process_channel_group( |
| 351 | orig_img=torch.Tensor(input_data[0, :, 0, ...]), |
| 352 | channels=channels, |
| 353 | ) |
| 354 | rgb_orig = rgb_orig.to(torch.float32) |
| 355 | |
| 356 | pred[pred == 0.0] = np.nan |
| 357 | img_pred = rgb_orig * 0.7 + pred * 0.3 |
| 358 | img_pred[img_pred.isnan()] = rgb_orig[img_pred.isnan()] |
| 359 | |
| 360 | img_pred_file = os.path.join( |
| 361 | output_dir, f"rgb_pred_{os.path.splitext(os.path.basename(data_file))[0]}.tiff" |
| 362 | ) |
| 363 | save_geotiff( |
| 364 | image=_convert_np_uint8(img_pred), |
| 365 | output_path=img_pred_file, |
| 366 | meta=meta_data, |
| 367 | ) |
| 368 | |
| 369 | # Save image rgb |
| 370 | if rgb_outputs: |
| 371 | name_suffix = os.path.splitext(os.path.basename(data_file))[0] |
| 372 | rgb_file = os.path.join( |
| 373 | output_dir, |
| 374 | f"original_rgb_{name_suffix}.tiff", |
| 375 | ) |
| 376 | save_geotiff( |
| 377 | image=_convert_np_uint8(rgb_orig), |
| 378 | output_path=rgb_file, |
| 379 | meta=meta_data, |
| 380 | ) |
| 381 | |
| 382 | |
| 383 | if __name__ == "__main__": |
| 384 | parser = argparse.ArgumentParser("MAE run inference", add_help=False) |
| 385 | |
| 386 | parser.add_argument( |
| 387 | "--data_file", |
| 388 | type=str, |
| 389 | default="./India_900498_S2Hand.tif", |
| 390 | help="Path to the file.", |
| 391 | ) |
| 392 | parser.add_argument( |
| 393 | "--model", |
| 394 | type=str, |
| 395 | default="ibm-nasa-geospatial/Prithvi-EO-2.0-300M-TL-Sen1Floods11", |
| 396 | help="Path to a checkpoint file to load from.", |
| 397 | ) |
| 398 | parser.add_argument( |
| 399 | "--output_dir", |
| 400 | type=str, |
| 401 | default="output", |
| 402 | help="Path to the directory where to save outputs.", |
| 403 | ) |
| 404 | parser.add_argument( |
| 405 | "--input_indices", |
| 406 | default=[1, 2, 3, 8, 11, 12], |
| 407 | type=int, |
| 408 | nargs="+", |
| 409 | help=""" |
| 410 | 0-based indices of the six Prithvi channels to be selected from the input. |
| 411 | By default selects [1,2,3,8,11,12] for S2L1C data. |
| 412 | """, |
| 413 | ) |
| 414 | parser.add_argument( |
| 415 | "--rgb_outputs", |
| 416 | action="store_true", |
| 417 | help="If present, output files will only contain RGB channels. " |
| 418 | "Otherwise, all bands will be saved.", |
| 419 | ) |
| 420 | args = parser.parse_args() |
| 421 | |
| 422 | main(**vars(args)) |