ARicci commited on
Commit
ffee91f
1 Parent(s): 0510723

Upload 7 files

Browse files
Prithvi.py ADDED
@@ -0,0 +1,319 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ # Copyright (c) Meta Platforms, Inc. and affiliates.
2
+ # All rights reserved.
3
+
4
+ # This source code is licensed under the license found in the
5
+ # LICENSE file in the root directory of this source tree.
6
+ # --------------------------------------------------------
7
+ # References:
8
+ # timm: https://github.com/rwightman/pytorch-image-models/tree/master/timm
9
+ # DeiT: https://github.com/facebookresearch/deit
10
+ # --------------------------------------------------------
11
+
12
+ from functools import partial
13
+
14
+ import torch
15
+ import torch.nn as nn
16
+
17
+ from timm.models.vision_transformer import Block
18
+ from timm.models.layers import to_2tuple
19
+
20
+ import numpy as np
21
+
22
+ from einops import rearrange
23
+
24
+ def get_1d_sincos_pos_embed_from_grid(embed_dim, pos):
25
+ """
26
+ embed_dim: output dimension for each position
27
+ pos: a list of positions to be encoded: size (M,)
28
+ out: (M, D)
29
+ """
30
+ assert embed_dim % 2 == 0
31
+ omega = np.arange(embed_dim // 2, dtype=np.float32)
32
+ omega /= embed_dim / 2.
33
+ omega = 1. / 10000**omega # (D/2,)
34
+
35
+ pos = pos.reshape(-1) # (M,)
36
+ out = np.einsum('m,d->md', pos, omega) # (M, D/2), outer product
37
+
38
+ emb_sin = np.sin(out) # (M, D/2)
39
+ emb_cos = np.cos(out) # (M, D/2)
40
+
41
+ emb = np.concatenate([emb_sin, emb_cos], axis=1) # (M, D)
42
+ return emb
43
+
44
+ def get_2d_sincos_pos_embed_from_grid(embed_dim, grid):
45
+ assert embed_dim % 2 == 0
46
+
47
+ # use half of dimensions to encode grid_h
48
+ emb_h = get_1d_sincos_pos_embed_from_grid(embed_dim // 2, grid[0]) # (H*W, D/2)
49
+ emb_w = get_1d_sincos_pos_embed_from_grid(embed_dim // 2, grid[1]) # (H*W, D/2)
50
+
51
+ emb = np.concatenate([emb_h, emb_w], axis=1) # (H*W, D)
52
+ return emb
53
+
54
+ def get_3d_sincos_pos_embed(embed_dim, grid_size, cls_token=False):
55
+ """
56
+ grid_size: 3d tuple of grid size: t, h, w
57
+ return:
58
+ pos_embed: L, D
59
+ """
60
+
61
+ assert embed_dim % 16 == 0
62
+
63
+ t_size, h_size, w_size = grid_size
64
+
65
+ w_embed_dim = embed_dim // 16 * 6
66
+ h_embed_dim = embed_dim // 16 * 6
67
+ t_embed_dim = embed_dim // 16 * 4
68
+
69
+ w_pos_embed = get_1d_sincos_pos_embed_from_grid(w_embed_dim, np.arange(w_size))
70
+ h_pos_embed = get_1d_sincos_pos_embed_from_grid(h_embed_dim, np.arange(h_size))
71
+ t_pos_embed = get_1d_sincos_pos_embed_from_grid(t_embed_dim, np.arange(t_size))
72
+
73
+ w_pos_embed = np.tile(w_pos_embed, (t_size * h_size, 1))
74
+ h_pos_embed = np.tile(np.repeat(h_pos_embed, w_size, axis=0), (t_size, 1))
75
+ t_pos_embed = np.repeat(t_pos_embed, h_size * w_size, axis=0)
76
+
77
+ pos_embed = np.concatenate((w_pos_embed, h_pos_embed, t_pos_embed), axis=1)
78
+
79
+ if cls_token:
80
+ pos_embed = np.concatenate([np.zeros([1, embed_dim]), pos_embed], axis=0)
81
+ return pos_embed
82
+
83
+
84
+ class PatchEmbed(nn.Module):
85
+ """ Frames of 2D Images to Patch Embedding
86
+ The 3D version of timm.models.vision_transformer.PatchEmbed
87
+ """
88
+ def __init__(
89
+ self,
90
+ img_size=224,
91
+ patch_size=16,
92
+ num_frames=3,
93
+ tubelet_size=1,
94
+ in_chans=3,
95
+ embed_dim=768,
96
+ norm_layer=None,
97
+ flatten=True,
98
+ bias=True,
99
+ ):
100
+ super().__init__()
101
+ img_size = to_2tuple(img_size)
102
+ patch_size = to_2tuple(patch_size)
103
+ self.img_size = img_size
104
+ self.patch_size = patch_size
105
+ self.num_frames = num_frames
106
+ self.tubelet_size = tubelet_size
107
+ self.grid_size = (num_frames // tubelet_size, img_size[0] // patch_size[0], img_size[1] // patch_size[1])
108
+ self.num_patches = self.grid_size[0] * self.grid_size[1] * self.grid_size[2]
109
+ self.flatten = flatten
110
+
111
+ self.proj = nn.Conv3d(in_chans, embed_dim,
112
+ kernel_size=(tubelet_size, patch_size[0], patch_size[1]),
113
+ stride=(tubelet_size, patch_size[0], patch_size[1]), bias=bias)
114
+ self.norm = norm_layer(embed_dim) if norm_layer else nn.Identity()
115
+
116
+ def forward(self, x):
117
+ B, C, T, H, W = x.shape
118
+ x = self.proj(x)
119
+ if self.flatten:
120
+ x = x.flatten(2).transpose(1, 2) # B,C,T,H,W -> B,C,L -> B,L,C
121
+ x = self.norm(x)
122
+ return x
123
+
124
+
125
+ class MaskedAutoencoderViT(nn.Module):
126
+ """ Masked Autoencoder with VisionTransformer backbone
127
+ """
128
+ def __init__(self, img_size=224, patch_size=16,
129
+ num_frames=3, tubelet_size=1,
130
+ in_chans=3, embed_dim=1024, depth=24, num_heads=16,
131
+ decoder_embed_dim=512, decoder_depth=8, decoder_num_heads=16,
132
+ mlp_ratio=4., norm_layer=nn.LayerNorm, norm_pix_loss=False):
133
+ super().__init__()
134
+
135
+ # --------------------------------------------------------------------------
136
+ # MAE encoder specifics
137
+ self.patch_embed = PatchEmbed(img_size, patch_size,num_frames, tubelet_size, in_chans, embed_dim)
138
+ num_patches = self.patch_embed.num_patches
139
+
140
+ self.cls_token = nn.Parameter(torch.zeros(1, 1, embed_dim))
141
+ self.pos_embed = nn.Parameter(torch.zeros(1, num_patches + 1, embed_dim), requires_grad=False) # fixed sin-cos embedding
142
+
143
+ self.blocks = nn.ModuleList([
144
+ Block(embed_dim, num_heads, mlp_ratio, qkv_bias=True, norm_layer=norm_layer)
145
+ for i in range(depth)])
146
+ self.norm = norm_layer(embed_dim)
147
+ # --------------------------------------------------------------------------
148
+
149
+ # --------------------------------------------------------------------------
150
+ # MAE decoder specifics
151
+ self.decoder_embed = nn.Linear(embed_dim, decoder_embed_dim, bias=True)
152
+
153
+ self.mask_token = nn.Parameter(torch.zeros(1, 1, decoder_embed_dim))
154
+
155
+ self.decoder_pos_embed = nn.Parameter(torch.zeros(1, num_patches + 1, decoder_embed_dim), requires_grad=False) # fixed sin-cos embedding
156
+
157
+ self.decoder_blocks = nn.ModuleList([
158
+ Block(decoder_embed_dim, decoder_num_heads, mlp_ratio, qkv_bias=True, norm_layer=norm_layer)
159
+ for i in range(decoder_depth)])
160
+
161
+ self.decoder_norm = norm_layer(decoder_embed_dim)
162
+ self.decoder_pred = nn.Linear(decoder_embed_dim, tubelet_size * patch_size * patch_size * in_chans, bias=True) # decoder to patch
163
+ # --------------------------------------------------------------------------
164
+
165
+ self.norm_pix_loss = norm_pix_loss
166
+
167
+ self.initialize_weights()
168
+
169
+ def initialize_weights(self):
170
+ # initialization
171
+ # initialize (and freeze) pos_embed by sin-cos embedding
172
+ pos_embed = get_3d_sincos_pos_embed(self.pos_embed.shape[-1], self.patch_embed.grid_size, cls_token=True)
173
+ self.pos_embed.data.copy_(torch.from_numpy(pos_embed).float().unsqueeze(0))
174
+
175
+ decoder_pos_embed = get_3d_sincos_pos_embed(self.decoder_pos_embed.shape[-1], self.patch_embed.grid_size, cls_token=True)
176
+ self.decoder_pos_embed.data.copy_(torch.from_numpy(decoder_pos_embed).float().unsqueeze(0))
177
+
178
+ # initialize patch_embed like nn.Linear (instead of nn.Conv2d)
179
+ w = self.patch_embed.proj.weight.data
180
+ torch.nn.init.xavier_uniform_(w.view([w.shape[0], -1]))
181
+
182
+ # timm's trunc_normal_(std=.02) is effectively normal_(std=0.02) as cutoff is too big (2.)
183
+ torch.nn.init.normal_(self.cls_token, std=.02)
184
+ torch.nn.init.normal_(self.mask_token, std=.02)
185
+
186
+ # initialize nn.Linear and nn.LayerNorm
187
+ self.apply(self._init_weights)
188
+
189
+ def _init_weights(self, m):
190
+ if isinstance(m, nn.Linear):
191
+ # we use xavier_uniform following official JAX ViT:
192
+ torch.nn.init.xavier_uniform_(m.weight)
193
+ if isinstance(m, nn.Linear) and m.bias is not None:
194
+ nn.init.constant_(m.bias, 0)
195
+ elif isinstance(m, nn.LayerNorm):
196
+ nn.init.constant_(m.bias, 0)
197
+ nn.init.constant_(m.weight, 1.0)
198
+
199
+ def patchify(self, imgs):
200
+ """
201
+ imgs: B, C, T, H, W
202
+ x: B, L, D
203
+ """
204
+ p = self.patch_embed.patch_size[0]
205
+ tub = self.patch_embed.tubelet_size
206
+ x = rearrange(imgs, 'b c (t tub) (h p) (w q) -> b (t h w) (tub p q c)', tub=tub, p=p, q=p)
207
+
208
+ return x
209
+
210
+ def unpatchify(self, x):
211
+ """
212
+ x: B, L, D
213
+ imgs: B, C, T, H, W
214
+ """
215
+ p = self.patch_embed.patch_size[0]
216
+ num_p = self.patch_embed.img_size[0] // p
217
+ tub = self.patch_embed.tubelet_size
218
+ imgs = rearrange(x, 'b (t h w) (tub p q c) -> b c (t tub) (h p) (w q)', h=num_p, w=num_p, tub=tub, p=p, q=p)
219
+ return imgs
220
+
221
+ def random_masking(self, x, mask_ratio):
222
+ """
223
+ Perform per-sample random masking by per-sample shuffling.
224
+ Per-sample shuffling is done by argsort random noise.
225
+ x: [N, L, D], sequence
226
+ """
227
+ N, L, D = x.shape # batch, length, dim
228
+ len_keep = int(L * (1 - mask_ratio))
229
+
230
+ noise = torch.rand(N, L, device=x.device) # noise in [0, 1]
231
+
232
+ # sort noise for each sample
233
+ ids_shuffle = torch.argsort(noise, dim=1) # ascend: small is keep, large is remove
234
+ ids_restore = torch.argsort(ids_shuffle, dim=1)
235
+
236
+ # keep the first subset
237
+ ids_keep = ids_shuffle[:, :len_keep]
238
+ x_masked = torch.gather(x, dim=1, index=ids_keep.unsqueeze(-1).repeat(1, 1, D))
239
+
240
+ # generate the binary mask: 0 is keep, 1 is remove
241
+ mask = torch.ones([N, L], device=x.device)
242
+ mask[:, :len_keep] = 0
243
+ # unshuffle to get the binary mask
244
+ mask = torch.gather(mask, dim=1, index=ids_restore)
245
+
246
+ return x_masked, mask, ids_restore
247
+
248
+ def forward_encoder(self, x, mask_ratio):
249
+ # embed patches
250
+ x = self.patch_embed(x)
251
+
252
+ # add pos embed w/o cls token
253
+ x = x + self.pos_embed[:, 1:, :]
254
+
255
+ # masking: length -> length * mask_ratio
256
+ x, mask, ids_restore = self.random_masking(x, mask_ratio)
257
+
258
+ # append cls token
259
+ cls_token = self.cls_token + self.pos_embed[:, :1, :]
260
+ cls_tokens = cls_token.expand(x.shape[0], -1, -1)
261
+ x = torch.cat((cls_tokens, x), dim=1)
262
+
263
+ # apply Transformer blocks
264
+ for blk in self.blocks:
265
+ x = blk(x)
266
+ x = self.norm(x)
267
+
268
+ return x, mask, ids_restore
269
+
270
+ def forward_decoder(self, x, ids_restore):
271
+ # embed tokens
272
+ x = self.decoder_embed(x)
273
+
274
+ # append mask tokens to sequence
275
+ mask_tokens = self.mask_token.repeat(x.shape[0], ids_restore.shape[1] + 1 - x.shape[1], 1)
276
+ x_ = torch.cat([x[:, 1:, :], mask_tokens], dim=1) # no cls token
277
+ x_ = torch.gather(x_, dim=1, index=ids_restore.unsqueeze(-1).repeat(1, 1, x.shape[2])) # unshuffle
278
+ x = torch.cat([x[:, :1, :], x_], dim=1) # append cls token
279
+
280
+ # add pos embed
281
+ x = x + self.decoder_pos_embed
282
+
283
+ # apply Transformer blocks
284
+ for blk in self.decoder_blocks:
285
+ x = blk(x)
286
+ x = self.decoder_norm(x)
287
+
288
+ # predictor projection
289
+ x = self.decoder_pred(x)
290
+
291
+ # remove cls token
292
+ x = x[:, 1:, :]
293
+
294
+ return x
295
+
296
+ def forward_loss(self, imgs, pred, mask):
297
+ """
298
+ imgs: B, C, T, H, W
299
+ target: B, L, D
300
+ pred: B, L, D
301
+ mask: B, L. 0 is keep, 1 is remove,
302
+ """
303
+ target = self.patchify(imgs)
304
+ if self.norm_pix_loss:
305
+ mean = target.mean(dim=-1, keepdim=True)
306
+ var = target.var(dim=-1, keepdim=True)
307
+ target = (target - mean) / (var + 1.e-6)**.5
308
+
309
+ loss = (pred - target) ** 2
310
+ loss = loss.mean(dim=-1) # [N, L], mean loss per patch
311
+
312
+ loss = (loss * mask).sum() / mask.sum() # mean loss on removed patches
313
+ return loss
314
+
315
+ def forward(self, imgs, mask_ratio=0.75):
316
+ latent, mask, ids_restore = self.forward_encoder(imgs, mask_ratio)
317
+ pred = self.forward_decoder(latent, ids_restore)
318
+ loss = self.forward_loss(imgs, pred, mask)
319
+ return loss, pred, mask
Prithvi_100M.pt ADDED
@@ -0,0 +1,3 @@
 
 
 
 
1
+ version https://git-lfs.github.com/spec/v1
2
+ oid sha256:69f8ac286f649d1bbed520f5c8560a60eba91d688f74e1a0f9aa8203b6fd62ab
3
+ size 453672901
Prithvi_100M_config.yaml ADDED
@@ -0,0 +1,36 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ model_args:
2
+ decoder_depth: 8
3
+ decoder_embed_dim: 512
4
+ decoder_num_heads: 16
5
+ depth: 12
6
+ embed_dim: 768
7
+ img_size: 224
8
+ in_chans: 6
9
+ num_frames: 3
10
+ num_heads: 12
11
+ patch_size: 16
12
+ tubelet_size: 1
13
+ train_params:
14
+ bands:
15
+ - B02
16
+ - B03
17
+ - B04
18
+ - B05
19
+ - B06
20
+ - B07
21
+ data_mean:
22
+ - 775.2290211032589
23
+ - 1080.992780391705
24
+ - 1228.5855250417867
25
+ - 2497.2022620507532
26
+ - 2204.2139147975554
27
+ - 1610.8324823273745
28
+ data_std:
29
+ - 1281.526139861424
30
+ - 1270.0297974547493
31
+ - 1399.4802505642526
32
+ - 1368.3446143747644
33
+ - 1291.6764008585435
34
+ - 1154.505683480695
35
+ mask_ratio: 0.75
36
+ random_cropping: true
Prithvi_run_inference.py ADDED
@@ -0,0 +1,517 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ import argparse
2
+ import functools
3
+ import os
4
+ from typing import List, Union
5
+
6
+ import numpy as np
7
+ import rasterio
8
+ import torch
9
+ import yaml
10
+ from einops import rearrange
11
+
12
+ from Prithvi import MaskedAutoencoderViT
13
+
14
+ NO_DATA = -9999
15
+ NO_DATA_FLOAT = 0.0001
16
+ PERCENTILES = (0.1, 99.9)
17
+
18
+
19
+ def process_channel_group(orig_img, new_img, channels, data_mean, data_std):
20
+ """Process *orig_img* and *new_img* for RGB visualization. Each band is rescaled back to the
21
+ original range using *data_mean* and *data_std* and then lowest and highest percentiles are
22
+ removed to enhance contrast. Data is rescaled to (0, 1) range and stacked channels_first.
23
+
24
+ Args:
25
+ orig_img: torch.Tensor representing original image (reference) with shape = (bands, H, W).
26
+ new_img: torch.Tensor representing image with shape = (bands, H, W).
27
+ channels: list of indices representing RGB channels.
28
+ data_mean: list of mean values for each band.
29
+ data_std: list of std values for each band.
30
+
31
+ Returns:
32
+ torch.Tensor with shape (num_channels, height, width) for original image
33
+ torch.Tensor with shape (num_channels, height, width) for the other image
34
+ """
35
+
36
+ stack_c = [], []
37
+
38
+ for c in channels:
39
+ orig_ch = orig_img[c, ...]
40
+ valid_mask = torch.ones_like(orig_ch, dtype=torch.bool)
41
+ valid_mask[orig_ch == NO_DATA_FLOAT] = False
42
+
43
+ # Back to original data range
44
+ orig_ch = (orig_ch * data_std[c]) + data_mean[c]
45
+ new_ch = (new_img[c, ...] * data_std[c]) + data_mean[c]
46
+
47
+ # Rescale (enhancing contrast)
48
+ min_value, max_value = np.percentile(orig_ch[valid_mask], PERCENTILES)
49
+
50
+ orig_ch = torch.clamp((orig_ch - min_value) / (max_value - min_value), 0, 1)
51
+ new_ch = torch.clamp((new_ch - min_value) / (max_value - min_value), 0, 1)
52
+
53
+ # No data as zeros
54
+ orig_ch[~valid_mask] = 0
55
+ new_ch[~valid_mask] = 0
56
+
57
+ stack_c[0].append(orig_ch)
58
+ stack_c[1].append(new_ch)
59
+
60
+ # Channels first
61
+ stack_orig = torch.stack(stack_c[0], dim=0)
62
+ stack_rec = torch.stack(stack_c[1], dim=0)
63
+
64
+ return stack_orig, stack_rec
65
+
66
+
67
+ def read_geotiff(file_path: str):
68
+ """Read all bands from *file_path* and return image + meta info.
69
+
70
+ Args:
71
+ file_path: path to image file.
72
+
73
+ Returns:
74
+ np.ndarray with shape (bands, height, width)
75
+ meta info dict
76
+ """
77
+
78
+ with rasterio.open(file_path) as src:
79
+ img = src.read()
80
+ meta = src.meta
81
+
82
+ return img, meta
83
+
84
+
85
+ def save_geotiff(image, output_path: str, meta: dict):
86
+ """Save multi-band image in Geotiff file.
87
+
88
+ Args:
89
+ image: np.ndarray with shape (bands, height, width)
90
+ output_path: path where to save the image
91
+ meta: dict with meta info.
92
+ """
93
+
94
+ with rasterio.open(output_path, "w", **meta) as dest:
95
+ for i in range(image.shape[0]):
96
+ dest.write(image[i, :, :], i + 1)
97
+
98
+ return
99
+
100
+
101
+ def _convert_np_uint8(float_image: torch.Tensor):
102
+ image = float_image.numpy() * 255.0
103
+ image = image.astype(dtype=np.uint8)
104
+
105
+ return image
106
+
107
+
108
+ def load_example(
109
+ file_paths: List[str],
110
+ mean: List[float],
111
+ std: List[float],
112
+ indices: Union[list[int], None] = None,
113
+ ):
114
+ """Build an input example by loading images in *file_paths*.
115
+
116
+ Args:
117
+ file_paths: list of file paths .
118
+ mean: list containing mean values for each band in the images in *file_paths*.
119
+ std: list containing std values for each band in the images in *file_paths*.
120
+
121
+ Returns:
122
+ np.array containing created example
123
+ list of meta info for each image in *file_paths*
124
+ """
125
+
126
+ imgs = []
127
+ metas = []
128
+
129
+ for file in file_paths:
130
+ img, meta = read_geotiff(file)
131
+
132
+ # Rescaling (don't normalize on nodata)
133
+ img = np.moveaxis(img, 0, -1) # channels last for rescaling
134
+ if indices is not None:
135
+ img = img[..., indices]
136
+ img = np.where(img == NO_DATA, NO_DATA_FLOAT, (img - mean) / std)
137
+
138
+ imgs.append(img)
139
+ metas.append(meta)
140
+
141
+ imgs = np.stack(imgs, axis=0) # num_frames, H, W, C
142
+ imgs = np.moveaxis(imgs, -1, 0).astype("float32") # C, num_frames, H, W
143
+ imgs = np.expand_dims(imgs, axis=0) # add batch dim
144
+
145
+ return imgs, metas
146
+
147
+
148
+ def run_model(
149
+ model: torch.nn.Module,
150
+ input_data: torch.Tensor,
151
+ mask_ratio: float,
152
+ device: torch.device,
153
+ ):
154
+ """Run *model* with *input_data* and create images from output tokens (mask, reconstructed + visible).
155
+
156
+ Args:
157
+ model: MAE model to run.
158
+ input_data: torch.Tensor with shape (B, C, T, H, W).
159
+ mask_ratio: mask ratio to use.
160
+ device: device where model should run.
161
+
162
+ Returns:
163
+ 3 torch.Tensor with shape (B, C, T, H, W).
164
+ """
165
+
166
+ with torch.no_grad():
167
+ x = input_data.to(device)
168
+
169
+ _, pred, mask = model(x, mask_ratio)
170
+
171
+ # Create mask and prediction images (un-patchify)
172
+ mask_img = (
173
+ model.unpatchify(mask.unsqueeze(-1).repeat(1, 1, pred.shape[-1])).detach().cpu()
174
+ )
175
+ pred_img = model.unpatchify(pred).detach().cpu()
176
+
177
+ # Mix visible and predicted patches
178
+ rec_img = input_data.clone()
179
+ rec_img[mask_img == 1] = pred_img[
180
+ mask_img == 1
181
+ ] # binary mask: 0 is keep, 1 is remove
182
+
183
+ # Switch zeros/ones in mask images so masked patches appear darker in plots (better visualization)
184
+ mask_img = (~(mask_img.to(torch.bool))).to(torch.float)
185
+
186
+ return rec_img, mask_img
187
+
188
+
189
+ def save_rgb_imgs(
190
+ input_img, rec_img, mask_img, channels, mean, std, output_dir, meta_data
191
+ ):
192
+ """Wrapper function to save Geotiff images (original, reconstructed, masked) per timestamp.
193
+
194
+ Args:
195
+ input_img: input torch.Tensor with shape (C, T, H, W).
196
+ rec_img: reconstructed torch.Tensor with shape (C, T, H, W).
197
+ mask_img: mask torch.Tensor with shape (C, T, H, W).
198
+ channels: list of indices representing RGB channels.
199
+ mean: list of mean values for each band.
200
+ std: list of std values for each band.
201
+ output_dir: directory where to save outputs.
202
+ meta_data: list of dicts with geotiff meta info.
203
+ """
204
+
205
+ for t in range(input_img.shape[1]):
206
+ rgb_orig, rgb_pred = process_channel_group(
207
+ orig_img=input_img[:, t, :, :],
208
+ new_img=rec_img[:, t, :, :],
209
+ channels=channels,
210
+ data_mean=mean,
211
+ data_std=std,
212
+ )
213
+
214
+ rgb_mask = mask_img[channels, t, :, :] * rgb_orig
215
+
216
+ # Saving images
217
+
218
+ save_geotiff(
219
+ image=_convert_np_uint8(rgb_orig),
220
+ output_path=os.path.join(output_dir, f"original_rgb_t{t}.tiff"),
221
+ meta=meta_data[t],
222
+ )
223
+
224
+ save_geotiff(
225
+ image=_convert_np_uint8(rgb_pred),
226
+ output_path=os.path.join(output_dir, f"predicted_rgb_t{t}.tiff"),
227
+ meta=meta_data[t],
228
+ )
229
+
230
+ save_geotiff(
231
+ image=_convert_np_uint8(rgb_mask),
232
+ output_path=os.path.join(output_dir, f"masked_rgb_t{t}.tiff"),
233
+ meta=meta_data[t],
234
+ )
235
+
236
+
237
+ def save_imgs(rec_img, mask_img, mean, std, output_dir, meta_data):
238
+ """Wrapper function to save Geotiff images (reconstructed, mask) per timestamp.
239
+
240
+ Args:
241
+ rec_img: reconstructed torch.Tensor with shape (C, T, H, W).
242
+ mask_img: mask torch.Tensor with shape (C, T, H, W).
243
+ mean: list of mean values for each band.
244
+ std: list of std values for each band.
245
+ output_dir: directory where to save outputs.
246
+ meta_data: list of dicts with geotiff meta info.
247
+ """
248
+
249
+ mean = torch.tensor(np.asarray(mean)[:, None, None]) # C H W
250
+ std = torch.tensor(np.asarray(std)[:, None, None])
251
+
252
+ for t in range(rec_img.shape[1]):
253
+ # Back to original data range
254
+ rec_img_t = ((rec_img[:, t, :, :] * std) + mean).to(torch.int16)
255
+
256
+ mask_img_t = mask_img[:, t, :, :].to(torch.int16)
257
+
258
+ # Saving images
259
+
260
+ save_geotiff(
261
+ image=rec_img_t,
262
+ output_path=os.path.join(output_dir, f"predicted_t{t}.tiff"),
263
+ meta=meta_data[t],
264
+ )
265
+
266
+ save_geotiff(
267
+ image=mask_img_t,
268
+ output_path=os.path.join(output_dir, f"mask_t{t}.tiff"),
269
+ meta=meta_data[t],
270
+ )
271
+
272
+
273
+ def main(
274
+ data_files: List[str],
275
+ yaml_file_path: str,
276
+ checkpoint: str,
277
+ output_dir: str,
278
+ rgb_outputs: bool,
279
+ img_size: int,
280
+ mask_ratio: float = None,
281
+ input_indices: list[int] = None,
282
+ ):
283
+ os.makedirs(output_dir, exist_ok=True)
284
+
285
+ # Get parameters --------
286
+
287
+ with open(yaml_file_path, "r") as f:
288
+ params = yaml.safe_load(f)
289
+
290
+ # data related
291
+ train_params = params["train_params"]
292
+ num_frames = len(data_files)
293
+ bands = train_params["bands"]
294
+ mean = train_params["data_mean"]
295
+ std = train_params["data_std"]
296
+
297
+ # model related
298
+ model_params = params["model_args"]
299
+ img_size = model_params["img_size"] if img_size is None else img_size
300
+ depth = model_params["depth"]
301
+ patch_size = model_params["patch_size"]
302
+ embed_dim = model_params["embed_dim"]
303
+ num_heads = model_params["num_heads"]
304
+ tubelet_size = model_params["tubelet_size"]
305
+ decoder_embed_dim = model_params["decoder_embed_dim"]
306
+ decoder_num_heads = model_params["decoder_num_heads"]
307
+ decoder_depth = model_params["decoder_depth"]
308
+
309
+ batch_size = 1
310
+
311
+ mask_ratio = train_params["mask_ratio"] if mask_ratio is None else mask_ratio
312
+
313
+ print(
314
+ f"\nTreating {len(data_files)} files as {len(data_files)} time steps from the same location\n"
315
+ )
316
+ if len(data_files) != 3:
317
+ print(
318
+ "The original model was trained for 3 time steps (expecting 3 files). \nResults with different numbers of timesteps may vary"
319
+ )
320
+
321
+ if torch.cuda.is_available():
322
+ device = torch.device("cuda")
323
+ else:
324
+ device = torch.device("cpu")
325
+
326
+ print(f"Using {device} device.\n")
327
+
328
+ # Loading data ---------------------------------------------------------------------------------
329
+
330
+ input_data, meta_data = load_example(
331
+ file_paths=data_files, indices=input_indices, mean=mean, std=std
332
+ )
333
+
334
+ # Create model and load checkpoint -------------------------------------------------------------
335
+
336
+ model = MaskedAutoencoderViT(
337
+ img_size=img_size,
338
+ patch_size=patch_size,
339
+ num_frames=num_frames,
340
+ tubelet_size=tubelet_size,
341
+ in_chans=len(bands),
342
+ embed_dim=embed_dim,
343
+ depth=depth,
344
+ num_heads=num_heads,
345
+ decoder_embed_dim=decoder_embed_dim,
346
+ decoder_depth=decoder_depth,
347
+ decoder_num_heads=decoder_num_heads,
348
+ mlp_ratio=4.0,
349
+ norm_layer=functools.partial(torch.nn.LayerNorm, eps=1e-6),
350
+ norm_pix_loss=False,
351
+ )
352
+
353
+ total_params = sum(p.numel() for p in model.parameters() if p.requires_grad)
354
+ print(f"\n--> Model has {total_params:,} parameters.\n")
355
+
356
+ model.to(device)
357
+
358
+ state_dict = torch.load(checkpoint, map_location=device)
359
+ # discard fixed pos_embedding weight
360
+ del state_dict["pos_embed"]
361
+ del state_dict["decoder_pos_embed"]
362
+ model.load_state_dict(state_dict, strict=False)
363
+ print(f"Loaded checkpoint from {checkpoint}")
364
+
365
+ # Running model --------------------------------------------------------------------------------
366
+
367
+ model.eval()
368
+ channels = [bands.index(b) for b in ["B04", "B03", "B02"]] # BGR -> RGB
369
+
370
+ # Reflect pad if not divisible by img_size
371
+ original_h, original_w = input_data.shape[-2:]
372
+ pad_h = img_size - (original_h % img_size)
373
+ pad_w = img_size - (original_w % img_size)
374
+ input_data = np.pad(
375
+ input_data, ((0, 0), (0, 0), (0, 0), (0, pad_h), (0, pad_w)), mode="reflect"
376
+ )
377
+
378
+ # Build sliding window
379
+ batch = torch.tensor(input_data, device="cpu")
380
+ windows = batch.unfold(3, img_size, img_size).unfold(4, img_size, img_size)
381
+ h1, w1 = windows.shape[3:5]
382
+ windows = rearrange(
383
+ windows, "b c t h1 w1 h w -> (b h1 w1) c t h w", h=img_size, w=img_size
384
+ )
385
+
386
+ # Split into batches if number of windows > batch_size
387
+ num_batches = windows.shape[0] // batch_size if windows.shape[0] > batch_size else 1
388
+ windows = torch.tensor_split(windows, num_batches, dim=0)
389
+
390
+ # Run model
391
+ rec_imgs = []
392
+ mask_imgs = []
393
+ for x in windows:
394
+ rec_img, mask_img = run_model(model, x, mask_ratio, device)
395
+ rec_imgs.append(rec_img)
396
+ mask_imgs.append(mask_img)
397
+
398
+ rec_imgs = torch.concat(rec_imgs, dim=0)
399
+ mask_imgs = torch.concat(mask_imgs, dim=0)
400
+
401
+ # Build images from patches
402
+ rec_imgs = rearrange(
403
+ rec_imgs,
404
+ "(b h1 w1) c t h w -> b c t (h1 h) (w1 w)",
405
+ h=img_size,
406
+ w=img_size,
407
+ b=1,
408
+ c=len(bands),
409
+ t=num_frames,
410
+ h1=h1,
411
+ w1=w1,
412
+ )
413
+ mask_imgs = rearrange(
414
+ mask_imgs,
415
+ "(b h1 w1) c t h w -> b c t (h1 h) (w1 w)",
416
+ h=img_size,
417
+ w=img_size,
418
+ b=1,
419
+ c=len(bands),
420
+ t=num_frames,
421
+ h1=h1,
422
+ w1=w1,
423
+ )
424
+
425
+ # Cut padded images back to original size
426
+ rec_imgs_full = rec_imgs[..., :original_h, :original_w]
427
+ mask_imgs_full = mask_imgs[..., :original_h, :original_w]
428
+ batch_full = batch[..., :original_h, :original_w]
429
+
430
+ # Build output images
431
+ if rgb_outputs:
432
+ for d in meta_data:
433
+ d.update(count=3, dtype="uint8", compress="lzw", nodata=0)
434
+
435
+ save_rgb_imgs(
436
+ batch_full[0, ...],
437
+ rec_imgs_full[0, ...],
438
+ mask_imgs_full[0, ...],
439
+ channels,
440
+ mean,
441
+ std,
442
+ output_dir,
443
+ meta_data,
444
+ )
445
+ else:
446
+ for d in meta_data:
447
+ d.update(compress="lzw", nodata=0)
448
+
449
+ save_imgs(
450
+ rec_imgs_full[0, ...],
451
+ mask_imgs_full[0, ...],
452
+ mean,
453
+ std,
454
+ output_dir,
455
+ meta_data,
456
+ )
457
+
458
+ print("Done!")
459
+
460
+
461
+ if __name__ == "__main__":
462
+ parser = argparse.ArgumentParser("MAE run inference", add_help=False)
463
+
464
+ parser.add_argument(
465
+ "--data_files",
466
+ required=True,
467
+ type=str,
468
+ nargs="+",
469
+ help="Path to the data files. Assumes multi-band files.",
470
+ )
471
+ parser.add_argument(
472
+ "--yaml_file_path",
473
+ type=str,
474
+ required=True,
475
+ help="Path to yaml file containing model training parameters.",
476
+ )
477
+ parser.add_argument(
478
+ "--checkpoint",
479
+ required=True,
480
+ type=str,
481
+ help="Path to a checkpoint file to load from.",
482
+ )
483
+ parser.add_argument(
484
+ "--output_dir",
485
+ required=True,
486
+ type=str,
487
+ help="Path to the directory where to save outputs.",
488
+ )
489
+ parser.add_argument(
490
+ "--mask_ratio",
491
+ default=None,
492
+ type=float,
493
+ help="Masking ratio (percentage of removed patches). "
494
+ "If None (default) use same value used for pretraining.",
495
+ )
496
+ parser.add_argument(
497
+ "--img_size",
498
+ default=224,
499
+ type=int,
500
+ help="Image size to be used with model. Defaults to 224",
501
+ )
502
+ parser.add_argument(
503
+ "--input_indices",
504
+ default=None,
505
+ type=int,
506
+ nargs="+",
507
+ help="0-based indices of channels to be selected from the input. By default takes all.",
508
+ )
509
+ parser.add_argument(
510
+ "--rgb_outputs",
511
+ action="store_true",
512
+ help="If present, output files will only contain RGB channels. "
513
+ "Otherwise, all bands will be saved.",
514
+ )
515
+ args = parser.parse_args()
516
+
517
+ main(**vars(args))
README.md CHANGED
@@ -1,3 +1,72 @@
1
  ---
2
  license: apache-2.0
 
 
 
 
 
3
  ---
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
  ---
2
  license: apache-2.0
3
+ tags:
4
+ - Pytorch
5
+ - Geospatial
6
+ - Temporal ViT
7
+ - Vit
8
  ---
9
+
10
+ ### Model and Inputs
11
+ Prithvi is a first-of-its-kind temporal Vision transformer pre-trained by the IBM and NASA team on contiguous US Harmonised Landsat Sentinel 2 (HLS) data. The model adopts a self-supervised encoder developed with a ViT architecture and Masked AutoEncoder (MAE) learning strategy, with an MSE loss function. The model includes spatial attention across multiple patches and also temporal attention for each patch.
12
+
13
+ ![](GFM.png)
14
+
15
+ The model accepts remote sensing data in a video format (B, C, T, H, W). Note that the temporal dimension (T) is very important in this application and not present in most other works around remote sensing modeling. The ability to handle a time series of remote sensing images can benefit a variety of downstream tasks (e.g. Burn Scars segmentation, Flood Segmentation, Land Cover Classification). The model can also handle static imagery which can be fed into the model with T=1.
16
+
17
+ ### Pre-training
18
+ The model was pre-trained with NASA's HLS V2 L30 product (30m granularity) from the contiguous United States. The bands that were used are the following:
19
+
20
+ 1. Blue
21
+ 2. Green
22
+ 3. Red
23
+ 4. Narrow NIR
24
+ 5. SWIR 1
25
+ 6. SWIR 2
26
+
27
+ ### Code
28
+ The model follows the [original MAE repo](https://github.com/facebookresearch/mae) with some modifications including:
29
+
30
+ 1. replace 2D patch embed with 3D patch embed;
31
+ 2. replace 2D positional embed with 3D positional embed;
32
+ 3. replace 2D patchify and unpatchify with 3D.
33
+ 4. adding infrared bands besides RGB
34
+
35
+ ### Inference and demo
36
+ There is an inference script (`Prithvi_run_inference.py`) that allows to run the image reconstruction on a set of HLS images assumed to be from the same location at different time steps(see example below). These should be provided in chronological order in geotiff format, including the channels described above (Blue, Green, Red, Narrow NIR, SWIR 1, SWIR 2) in reflectance units. There is also a **demo** that leverages the same code [here](https://huggingface.co/spaces/ibm-nasa-geospatial/Prithvi-100M-demo).
37
+
38
+ ```
39
+ python Prithvi_run_inference.py --data_files t1.tif t2.tif t3.tif --yaml_file_path /path/to/yaml/Prithvi_100.yaml --checkpoint /path/to/checkpoint/Prithvi_100.pth --output_dir /path/to/out/dir/ --input_indices <space separated 0-based indices of channels to select from input> --mask_ratio 0.5 --img_size <length of one side of square input shape>
40
+ ```
41
+
42
+ This demo is a starting point that can be used as a starting point to generalize to different input shapes / types.
43
+
44
+ ### Finetuning examples
45
+ Examples of finetuning the model for image segmentation using the mmsegmentation library are available through Hugging Face (e.g. [burn scars segmentation](https://huggingface.co/ibm-nasa-geospatial/Prithvi-100M-burn-scar), [flood mapping](https://huggingface.co/ibm-nasa-geospatial/Prithvi-100M-sen1floods11), and [multi temporal crop classification](https://huggingface.co/ibm-nasa-geospatial/Prithvi-100M-multi-temporal-crop-classification)), with the code used for the experiments available on [github](https://github.com/NASA-IMPACT/hls-foundation-os/tree/main/fine-tuning-examples). This also contains instructions to finetune the model for flood detection on the popular open access [sen1floods11 dataset](https://github.com/cloudtostreet/Sen1Floods11).
46
+
47
+ ### Feedback
48
+
49
+ Your feedback is invaluable to us. If you have any feedback about the model, please feel free to share it with us. You can do this by submitting issues on our open-source repository, [hls-foundation-os](https://github.com/NASA-IMPACT/hls-foundation-os/issues), on GitHub.
50
+
51
+ ### Citation
52
+
53
+ If this model helped your research, please cite `Prithvi-100M` in your publications. Here are two BibTeX entries as examples:
54
+
55
+ ```
56
+ @article{Prithvi-100M-preprint,
57
+ author = {Jakubik, Johannes and Roy, Sujit and Phillips, C. E. and Fraccaro, Paolo and Godwin, Denys and Zadrozny, Bianca and Szwarcman, Daniela and Gomes, Carlos and Nyirjesy, Gabby and Edwards, Blair and Kimura, Daiki and Simumba, Naomi and Chu, Linsong and Mukkavilli, S. Karthik and Lambhate, Devyani and Das, Kamal and Bangalore, Ranjini and Oliveira, Dario and Muszynski, Michal and Ankur, Kumar and Ramasubramanian, Muthukumaran and Gurung, Iksha and Khallaghi, Sam and Li, Hanxi (Steve) and Cecil, Michael and Ahmadi, Maryam and Kordi, Fatemeh and Alemohammad, Hamed and Maskey, Manil and Ganti, Raghu and Weldemariam, Kommy and Ramachandran, Rahul},
58
+ month = oct,
59
+ title = {{Foundation Models for Generalist Geospatial Artificial Intelligence}},
60
+ journal = {Preprint Available on arxiv:2310.18660},
61
+ year = {2023}
62
+ }
63
+
64
+ @misc{Prithvi-100M,
65
+ author = {Jakubik, Johannes and Chu, Linsong and Fraccaro, Paolo and Gomes, Carlos and Nyirjesy, Gabby and Bangalore, Ranjini and Lambhate, Devyani and Das, Kamal and Oliveira Borges, Dario and Kimura, Daiki and Simumba, Naomi and Szwarcman, Daniela and Muszynski, Michal and Weldemariam, Kommy and Zadrozny, Bianca and Ganti, Raghu and Costa, Carlos and Edwards, Blair & Watson, Campbell and Mukkavilli, Karthik and Schmude, Johannes & Hamann, Hendrik and Robert, Parkin and Roy, Sujit and Phillips, Christopher and Ankur, Kumar and Ramasubramanian, Muthukumaran and Gurung, Iksha and Leong, Wei Ji and Avery, Ryan and Ramachandran, Rahul and Maskey, Manil and Olofossen, Pontus and Fancher, Elizabeth and Lee, Tsengdar and Murphy, Kevin and Duffy, Dan and Little, Mike and Alemohammad, Hamed and Cecil, Michael and Li, Steve and Khallaghi, Sam and Godwin, Denys and Ahmadi, Maryam and Kordi, Fatemeh and Saux, Bertrand and Pastick, Neal and Doucette, Peter and Fleckenstein, Rylie and Luanga, Dalton and Corvin, Alex and Granger, Erwan},
66
+ doi = {10.57967/hf/0952},
67
+ month = aug,
68
+ title = {{Prithvi-100M}},
69
+ repository-code = {https://github.com/NASA-IMPACT/hls-foundation-os},
70
+ year = {2023}
71
+ }
72
+ ```
config.json ADDED
@@ -0,0 +1,43 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ {
2
+ "model_args": {
3
+ "decoder_depth": 8,
4
+ "decoder_embed_dim": 512,
5
+ "decoder_num_heads": 16,
6
+ "depth": 12,
7
+ "embed_dim": 768,
8
+ "img_size": 224,
9
+ "in_chans": 6,
10
+ "num_frames": 3,
11
+ "num_heads": 12,
12
+ "patch_size": 16,
13
+ "tubelet_size": 1
14
+ },
15
+ "train_params": {
16
+ "bands": [
17
+ "B02",
18
+ "B03",
19
+ "B04",
20
+ "B05",
21
+ "B06",
22
+ "B07"
23
+ ],
24
+ "data_mean": [
25
+ 775.2290211032589,
26
+ 1080.992780391705,
27
+ 1228.5855250417867,
28
+ 2497.2022620507532,
29
+ 2204.2139147975554,
30
+ 1610.8324823273745
31
+ ],
32
+ "data_std": [
33
+ 1281.526139861424,
34
+ 1270.0297974547493,
35
+ 1399.4802505642526,
36
+ 1368.3446143747644,
37
+ 1291.6764008585435,
38
+ 1154.505683480695
39
+ ],
40
+ "mask_ratio": 0.75,
41
+ "random_cropping": true
42
+ }
43
+ }
requirements.txt ADDED
@@ -0,0 +1,5 @@
 
 
 
 
 
 
1
+ torch
2
+ torchvision
3
+ timm
4
+ einops
5
+ rasterio