update _ 03 chạy thành công CNN

.vscode/settings.json

@@ -1,5 +0,0 @@
-{
-    "python-envs.defaultEnvManager": "ms-python.python:conda",
-    "python-envs.defaultPackageManager": "ms-python.python:conda",
-    "python-envs.pythonProjects": []
-}

cnn_model.py

@@ -0,0 +1,274 @@
+"""
+PyTorch 1D CNN Model for Land Use Classification
+Dùng cho 3 channels: NDVI, VH, VV (13 time steps)
+"""
+
+import torch
+import torch.nn as nn
+import torch.optim as optim
+from torch.utils.data import Dataset, DataLoader
+import numpy as np
+from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix
+import matplotlib.pyplot as plt
+
+
+class TimeSeriesDataset(Dataset):
+    """Custom Dataset for time series data"""
+    def __init__(self, X, y):
+        self.X = torch.FloatTensor(X)
+        self.y = torch.LongTensor(y)
+    
+    def __len__(self):
+        return len(self.X)
+    
+    def __getitem__(self, idx):
+        return self.X[idx], self.y[idx]
+
+
+class CNN1D(nn.Module):
+    """1D CNN for time series classification
+    
+    Input: (batch_size, 3, 13)
+        - 3 channels: NDVI, VH, VV
+        - 13 time steps
+    Output: (batch_size, num_classes)
+    """
+    def __init__(self, num_classes=8, dropout_rate=0.3):
+        super(CNN1D, self).__init__()
+        
+        # 1D Convolutional layers
+        self.conv1 = nn.Conv1d(in_channels=3, out_channels=32, kernel_size=3, padding=1)
+        self.bn1 = nn.BatchNorm1d(32)
+        self.relu1 = nn.ReLU()
+        self.pool1 = nn.MaxPool1d(kernel_size=2, stride=2)
+        
+        self.conv2 = nn.Conv1d(in_channels=32, out_channels=64, kernel_size=3, padding=1)
+        self.bn2 = nn.BatchNorm1d(64)
+        self.relu2 = nn.ReLU()
+        self.pool2 = nn.MaxPool1d(kernel_size=2, stride=2)
+        
+        self.conv3 = nn.Conv1d(in_channels=64, out_channels=128, kernel_size=3, padding=1)
+        self.bn3 = nn.BatchNorm1d(128)
+        self.relu3 = nn.ReLU()
+        self.pool3 = nn.MaxPool1d(kernel_size=2, stride=2)
+        
+        # Global average pooling
+        self.global_avg_pool = nn.AdaptiveAvgPool1d(1)
+        
+        # Fully connected layers
+        self.fc1 = nn.Linear(128, 64)
+        self.dropout = nn.Dropout(dropout_rate)
+        self.fc2 = nn.Linear(64, num_classes)
+    
+    def forward(self, x):
+        # Conv block 1
+        x = self.conv1(x)
+        x = self.bn1(x)
+        x = self.relu1(x)
+        x = self.pool1(x)
+        
+        # Conv block 2
+        x = self.conv2(x)
+        x = self.bn2(x)
+        x = self.relu2(x)
+        x = self.pool2(x)
+        
+        # Conv block 3
+        x = self.conv3(x)
+        x = self.bn3(x)
+        x = self.relu3(x)
+        x = self.pool3(x)
+        
+        # Global average pooling
+        x = self.global_avg_pool(x)
+        x = x.view(x.size(0), -1)
+        
+        # FC layers
+        x = self.fc1(x)
+        x = self.dropout(x)
+        x = self.fc2(x)
+        
+        return x
+
+
+class CNNTrainer:
+    """Trainer for PyTorch CNN Model"""
+    
+    def __init__(self, num_classes=8, learning_rate=0.001, device=None):
+        self.device = device or torch.device("cuda" if torch.cuda.is_available() else "cpu")
+        self.model = CNN1D(num_classes=num_classes).to(self.device)
+        self.criterion = nn.CrossEntropyLoss()
+        self.optimizer = optim.Adam(self.model.parameters(), lr=learning_rate)
+        self.history = {'train_loss': [], 'val_loss': [], 'train_acc': [], 'val_acc': []}
+        
+        print(f"🚀 Model initialized on device: {self.device}")
+        print(f"   Total parameters: {sum(p.numel() for p in self.model.parameters()):,}")
+    
+    def train_epoch(self, train_loader):
+        """Train one epoch"""
+        self.model.train()
+        total_loss = 0
+        correct = 0
+        total = 0
+        
+        for X_batch, y_batch in train_loader:
+            X_batch, y_batch = X_batch.to(self.device), y_batch.to(self.device)
+            
+            # Forward pass
+            outputs = self.model(X_batch)
+            loss = self.criterion(outputs, y_batch)
+            
+            # Backward pass
+            self.optimizer.zero_grad()
+            loss.backward()
+            self.optimizer.step()
+            
+            # Metrics
+            total_loss += loss.item()
+            _, predicted = torch.max(outputs.data, 1)
+            correct += (predicted == y_batch).sum().item()
+            total += y_batch.size(0)
+        
+        avg_loss = total_loss / len(train_loader)
+        accuracy = correct / total
+        return avg_loss, accuracy
+    
+    def validate(self, val_loader):
+        """Validate model"""
+        self.model.eval()
+        total_loss = 0
+        correct = 0
+        total = 0
+        
+        with torch.no_grad():
+            for X_batch, y_batch in val_loader:
+                X_batch, y_batch = X_batch.to(self.device), y_batch.to(self.device)
+                
+                outputs = self.model(X_batch)
+                loss = self.criterion(outputs, y_batch)
+                
+                total_loss += loss.item()
+                _, predicted = torch.max(outputs.data, 1)
+                correct += (predicted == y_batch).sum().item()
+                total += y_batch.size(0)
+        
+        avg_loss = total_loss / len(val_loader)
+        accuracy = correct / total
+        return avg_loss, accuracy
+    
+    def fit(self, X_train, y_train, X_val, y_val, epochs=50, batch_size=32, verbose=True):
+        """Train model with validation"""
+        train_dataset = TimeSeriesDataset(X_train, y_train)
+        val_dataset = TimeSeriesDataset(X_val, y_val)
+        
+        train_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True)
+        val_loader = DataLoader(val_dataset, batch_size=batch_size, shuffle=False)
+        
+        print(f"\n📊 Training start:")
+        print(f"   Train samples: {len(X_train)}")
+        print(f"   Val samples: {len(X_val)}")
+        print(f"   Batch size: {batch_size}")
+        print(f"   Epochs: {epochs}\n")
+        
+        for epoch in range(epochs):
+            train_loss, train_acc = self.train_epoch(train_loader)
+            val_loss, val_acc = self.validate(val_loader)
+            
+            self.history['train_loss'].append(train_loss)
+            self.history['val_loss'].append(val_loss)
+            self.history['train_acc'].append(train_acc)
+            self.history['val_acc'].append(val_acc)
+            
+            if verbose and (epoch + 1) % 10 == 0:
+                print(f"Epoch [{epoch+1}/{epochs}] "
+                      f"Train Loss: {train_loss:.4f}, Acc: {train_acc:.4f} | "
+                      f"Val Loss: {val_loss:.4f}, Acc: {val_acc:.4f}")
+        
+        print(f"\n✅ Training completed!")
+        print(f"   Final Train Acc: {train_acc:.4f}")
+        print(f"   Final Val Acc: {val_acc:.4f}")
+    
+    def predict(self, X_test):
+        """Predict on test data"""
+        self.model.eval()
+        X_test = torch.FloatTensor(X_test).to(self.device)
+        
+        with torch.no_grad():
+            outputs = self.model(X_test)
+            _, predictions = torch.max(outputs, 1)
+        
+        return predictions.cpu().numpy()
+    
+    def evaluate(self, X_test, y_test):
+        """Evaluate on test data"""
+        y_pred = self.predict(X_test)
+        
+        accuracy = accuracy_score(y_test, y_pred)
+        precision = precision_score(y_test, y_pred, average='weighted', zero_division=0)
+        recall = recall_score(y_test, y_pred, average='weighted', zero_division=0)
+        f1 = f1_score(y_test, y_pred, average='weighted', zero_division=0)
+        
+        print(f"\n📈 Test Results:")
+        print(f"   Accuracy:  {accuracy:.4f}")
+        print(f"   Precision: {precision:.4f}")
+        print(f"   Recall:    {recall:.4f}")
+        print(f"   F1-Score:  {f1:.4f}")
+        
+        return {
+            'accuracy': accuracy,
+            'precision': precision,
+            'recall': recall,
+            'f1': f1,
+            'predictions': y_pred,
+            'confusion_matrix': confusion_matrix(y_test, y_pred)
+        }
+    
+    def plot_history(self):
+        """Plot training history"""
+        fig, axes = plt.subplots(1, 2, figsize=(12, 4))
+        
+        # Loss
+        axes[0].plot(self.history['train_loss'], label='Train Loss')
+        axes[0].plot(self.history['val_loss'], label='Val Loss')
+        axes[0].set_xlabel('Epoch')
+        axes[0].set_ylabel('Loss')
+        axes[0].set_title('Training and Validation Loss')
+        axes[0].legend()
+        axes[0].grid(True)
+        
+        # Accuracy
+        axes[1].plot(self.history['train_acc'], label='Train Acc')
+        axes[1].plot(self.history['val_acc'], label='Val Acc')
+        axes[1].set_xlabel('Epoch')
+        axes[1].set_ylabel('Accuracy')
+        axes[1].set_title('Training and Validation Accuracy')
+        axes[1].legend()
+        axes[1].grid(True)
+        
+        plt.tight_layout()
+        plt.show()
+    
+    def save(self, filepath):
+        """Save model"""
+        torch.save(self.model.state_dict(), filepath)
+        print(f"✅ Model saved to {filepath}")
+    
+    def load(self, filepath):
+        """Load model"""
+        self.model.load_state_dict(torch.load(filepath, map_location=self.device))
+        print(f"✅ Model loaded from {filepath}")
+
+
+def reshape_for_cnn(X):
+    """Reshape data for CNN
+    
+    Input: (n_samples, n_features) where n_features = 13*3 = 39 (13 timesteps x 3 channels)
+    Output: (n_samples, 3, 13) - (batch, channels, timesteps)
+    """
+    n_samples = X.shape[0]
+    n_timesteps = 13
+    n_channels = 3
+    
+    # Reshape: (n_samples, 39) -> (n_samples, 3, 13)
+    X_cnn = X.reshape(n_samples, n_channels, n_timesteps)
+    return X_cnn

easi_tools/__init__.py

@@ -1,9 +0,0 @@
-#!python 3
-
-from .deployments import EasiDefaults
-from .notebook_utils import \
-    heading, \
-    initialize_dask, \
-    mostcommon_crs, \
-    unset_cachingproxy, \
-    xarray_object_size

easi_tools/deployments.py

@@ -1,335 +0,0 @@
-#!python3
-
-import sys
-import os
-import logging
-import collections
-
-# A class that provides notebook variables for each of the EASI deployments
-
-# Map an internal deployment name to deployment variables and search parameters.
-# Update to ensure that the product/space/time parameters are available in the respective databases
-deployment_map = {
-    'adias': {
-        'domain': 'adias.aquawatchaus.space',
-        'db_database': 'adias_prod_db',
-        'training_shapefile': '',
-        'scratch': 'adias-prod-user-scratch',
-        'productmap': {'landsat': 'landsat8_c2l2_sr', 'sentinel-2': 's2_l2a', 'sar': 'asf_s1_grd_gamma0', 'dem': 'copernicus_dem_30'},
-        'location': 'Lake Tahoe, California',
-        'latitude': (39.0, 39.3),
-        'longitude': (-120.2, -119.9),
-        'time': ('2022-02-01', '2022-05-01'),
-        'target': {
-            'landsat': {'crs': 'epsg:26911', 'resolution': (-30,30)},
-            'sentinel-2': {'crs': 'epsg:26911', 'resolution': (-10,10)}
-        },
-        'aliases': {
-            'landsat': {'qa_band': 'qa_pixel', 'nir': 'nir08', 'swir1': 'swir16', 'swir2': 'swir22'}
-        },
-        'qa_mask': {
-            'landsat': {'nodata': False, 'water': 'land_or_cloud',
-                        'cloud': 'not_high_confidence', 'cloud_shadow': 'not_high_confidence'}
-        }
-    },
-    'asia': {
-        'domain': 'asia.easi-eo.solutions',
-        'db_database': 'easi_asia_db',
-        'training_shapefile': '',
-        'scratch': 'easi-asia-user-scratch',
-        'productmap': {'landsat': 'landsat8_c2l2_sr', 'sentinel-2': 'sentinel_2_c1_l2a', 'sentinel-1': 'sentinel1_grd_gamma0_20m', 'dem': 'copernicus_dem_30'},
-        'location': 'Lake Tempe, Indonesia',
-        'latitude': (-4.2, -3.9),
-        'longitude': (119.8, 120.1),
-        'time': ('2020-02-01', '2020-04-01'),
-        'proxy': True,
-        'target': {
-            'landsat': {'crs': 'epsg:32650', 'resolution': (-30,30)},
-            'sentinel-2': {'crs': 'epsg:32650', 'resolution': (-10,10)}
-        },
-        'aliases': {
-            'landsat': {'qa_band': 'qa_pixel', 'nir': 'nir08', 'swir1': 'swir16', 'swir2': 'swir22'}
-        },
-        'qa_mask': {
-            'landsat': {'nodata': False, 'water': 'land_or_cloud',
-                        'cloud': 'not_high_confidence', 'cloud_shadow': 'not_high_confidence'}
-        }
-    },
-    'chile': {
-        'domain': 'datacubechile.cl',
-        'db_database': 'easido_prod_db',
-        'training_shapefile': '',
-        'scratch': 'easido-prod-user-scratch',
-        'productmap': {'landsat': 'landsat8_c2l2_sr', 'sentinel-2': 'sentinel_2_c1_l2a', 'sar': 'asf_s1_grd_gamma0', 'dem': 'copernicus_dem_30'},
-        'location': 'La Serena, Chile',
-        'latitude': (-29.95, -29.85),
-        'longitude': (-71.3, -71.2),
-        'latitude_big': (-29.95, -27.95),
-        'longitude_big': (-71.3, -69.3),
-        'time': ('2022-02-01', '2022-05-01'),
-        'target': {
-            'landsat': {'crs': 'epsg:32718', 'resolution': (-30,30)},
-            'sentinel-2': {'crs': 'epsg:32718', 'resolution': (-10,10)}
-        },
-        'aliases': {
-            'landsat': {'qa_band': 'qa_pixel', 'nir': 'nir08', 'swir1': 'swir16', 'swir2': 'swir22'}
-        },
-        'qa_mask': {
-            'landsat': {'nodata': False, 'water': 'land_or_cloud',
-                        'cloud': 'not_high_confidence', 'cloud_shadow': 'not_high_confidence'}
-        }
-    },
-    'cal': {
-        'domain': 'cal.ceos.org',
-        'db_database': 'ceoseail_eail_db',
-        'training_shapefile': './ancillary_data/VA_Counties_Newport_News.shp',
-        'scratch': 'ceoseail-eail-user-scratch',
-        'ows': False,
-        'map': False,
-        'productmap': {'landsat': 'landsat8_c2l2_sr', 'sentinel-2': 's2_l2a', 'sentinel-1': 's1_rtc', 'dem': 'copernicus_dem_30'},
-        'location': 'Newport News, Virginia',
-        'latitude': (37.02, 37.12),
-        'longitude': (-76.55, -76.45),
-        'time': ('2022-01-01', '2022-04-01'),
-        'target': {
-            'landsat': {'crs': 'epsg:32618', 'resolution': (-30,30)},
-            'sentinel-2': {'crs': 'epsg:32618', 'resolution': (-10,10)}
-        },
-        'aliases': {
-            'landsat': {'qa_band': 'qa_pixel', 'nir': 'nir08', 'swir1': 'swir16', 'swir2': 'swir22'}
-        },
-        'qa_mask': {
-            'landsat': {'nodata': False, 'water': 'land_or_cloud',
-                        'cloud': 'not_high_confidence', 'cloud_shadow': 'not_high_confidence'}
-        }
-    },
-    'csiro': {
-        'domain': 'csiro.easi-eo.solutions',
-        'db_database': 'easihub_csiro_db',
-        'training_shapefile': '',
-        'scratch': 'easihub-csiro-user-scratch',
-        'productmap': {'landsat': 'ga_ls8c_ard_3', 'sentinel-2': 'ga_s2am_ard_3', 'sentinel-1': 'sentinel1_grd_gamma0_20m', 'dem': 'copernicus_dem_30'},
-        'location': 'Lake Hume, Australia',
-        'latitude': (-36.3, -35.8),
-        'longitude': (146.8, 147.3),
-        'time': ('2020-02-01', '2020-04-01'),
-        'aliases': {
-            'landsat': {'red': 'nbart_red', 'green': 'nbart_green', 'blue': 'nbart_blue',
-                        'nir': 'nbart_nir', 'swir1': 'nbart_swir_1', 'swir2': 'nbart_swir_2',
-                        'qa_band': 'oa_fmask'}
-        },
-        'qa_mask': {
-            'landsat': {'fmask':'valid'}
-        }
-    },
-    'sub-apse2': {
-        'domain': 'sub-apse2.easi-eo.solutions',
-        'db_database': '',
-        'training_shapefile': '',
-        'scratch': '',
-        'ows': False,
-        'map': False,
-        'productmap': {'landsat': 'ga_ls8c_ard_3', 'sentinel-2': 'ga_s2am_ard_3', 'dem': 'copernicus_dem_30'},
-        'location': 'Lake Hume, Australia',
-        'latitude': (-36.3, -35.8),
-        'longitude': (146.8, 147.3),
-        'time': ('2020-02-01', '2020-04-01'),
-        'aliases': {
-            'landsat': {'red': 'nbart_red', 'green': 'nbart_green', 'blue': 'nbart_blue',
-                        'nir': 'nbart_nir', 'swir1': 'nbart_swir_1', 'swir2': 'nbart_swir_2',
-                        'qa_band': 'oa_fmask'}
-        },
-        'qa_mask': {
-            'landsat': {'fmask':'valid'}
-        }
-    },
-}
-
-
-class EasiDefaults():
-    """Provide deployment-specific default variables for EASI notebooks"""
-
-    def __init__(self, deployment=None):
-        """Initialise"""
-        self._log = _getlogger(self.__class__.__name__)
-        self.name = deployment if deployment else self._find_deployment()
-        self.deployment = self._validate(self.name)
-        self.proxy = None
-        self._aliases = {}
-        if self.deployment and self.deployment.get('proxy', None):
-            self.proxy = EasiCachingProxy()
-        if self.deployment:
-            self._log.info(f'Successfully found configuration for deployment "{self.name}"')
-
-    def _validate(self, deployment) -> dict:
-        """Return the dict associated with the deployment name"""
-        names = deployment_map.keys()
-        if deployment is None or deployment not in names:
-            self._log.error(f'Deployment name not recognised: {deployment}')
-            self._log.error(f'Select one of: {", ".join(names)}')
-            return None
-        return deployment_map[deployment]
-
-    def _find_deployment(self) -> str:
-        """Use the deployment's database environment variable as a lookup into the deployment_map dict"""
-        db_database = os.environ['DB_DATABASE']
-        deployment_name = [item for item in deployment_map if deployment_map[item]["db_database"] == db_database]
-        msg = 'Try specifying one using EasiDefaults(deployment="deployment_name").'
-        if len(deployment_name) == 0:
-            self._log.error(f'Deployment could not be found automatically. {msg}')
-            return None
-        elif len(deployment_name) > 1:
-            self._log.error(f'More than one deployment found. {msg}')
-            return None
-        return deployment_name[0]
-
-
-    @property
-    def domain(self):
-        """Deployment domain"""
-        return self.deployment['domain']
-
-    @property
-    def db_database(self):
-        """Database name"""
-        return self.deployment['db_database']
-
-    @property
-    def training_shapefile(self):
-        """A local shapefile"""
-        return self.deployment['training_shapefile']
-
-    @property
-    def hub(self):
-        """JupyterLab URL"""
-        return f'https://hub.{self.domain}'
-
-    @property
-    def explorer(self):
-        """Explorer URL"""
-        return f'https://explorer.{self.domain}'
-
-    @property
-    def ows(self):
-        """OWS URL"""
-        if not self.deployment.get('ows', True):
-            self._log.warning(f'Deployment does not have an OWS service: {self.name}')
-            return None
-        return f'https://ows.{self.domain}'
-
-    @property
-    def terria(self):
-        """Terria Map URL"""
-        if not self.deployment.get('map', True):
-            self._log.warning(f'Deployment does not have a Map service: {self.name}')
-            return None
-        return f'https://map.{self._domain()}'
-
-    @property
-    def scratch(self):
-        """Scratch bucket"""
-        return self.deployment['scratch']
-
-    @property
-    def location(self):
-        """Default location name"""
-        return self.deployment['location']
-
-    @property
-    def latitude(self):
-        """Default latitude range"""
-        return self.deployment['latitude']
-
-    @property
-    def longitude(self):
-        """Default longitude range"""
-        return self.deployment['longitude']
-
-    @property
-    def latitude_big(self):
-        """Default big latitude range"""
-        if 'latitude_big' in self.deployment:
-            return self.deployment['latitude_big']
-        self._log.warning(f'Default big latitude range not defined for "{self.deployment}". Using default latitude range')
-        return self.latitude
-
-    @property
-    def longitude_big(self):
-        """Default big longitude range"""
-        if 'longitude_big' in self.deployment:
-            return self.deployment['longitude_big']
-        self._log.warning(f'Default big longitude range not defined for "{self.deployment}". Using default longitude range')
-        return self.latitude
-
-    @property
-    def time(self):
-        """Default time range"""
-        return self.deployment['time']
-
-    def product(self, family='landsat'):
-        """Product name. Family loosely describes products from a satellite series or product type."""
-        p = self.deployment['productmap'].get(family, None)
-        if p is None:
-            self._log.warning(f'Product family not defined for "{self.name}": {family}')
-            out = ', '.join([f'{k} > {v}' for k,v in self.deployment['productmap'].items()])
-            self._log.warning(f'{self.name}: {out}')
-            return None
-        return p
-
-    def crs(self, family='landsat'):
-        """Default resolution. Family loosely describes products from a satellite series or product type."""
-        return self.deployment.get('target', {}).get(family, {}).get('crs', None)
-
-    def resolution(self, family='landsat'):
-        """Default resolution. Family loosely describes products from a satellite series or product type."""
-        return self.deployment.get('target', {}).get(family, {}).get('resolution', None)
-
-    def aliases(self, family='landsat') -> collections.UserDict:
-        """Return a dict-like object that maps a common name to a specific measurement/alias name.
-        Family loosely describes products from a satellite series or product type.
-
-        The common name is returned if there is no specific measurement/alias name defined.
-        That is, the common name should work as a measurement/alias name for the family in this deployment.
-        Else, provide a specific measurement/alias name in the defaults above.
-        """
-        if family not in self._aliases:
-            self._aliases[family] = EasiAlias(self.deployment.get('aliases', {}).get(family, {}))
-        return self._aliases[family]
-
-    def qa_mask(self, family='landsat') -> dict:
-        """Default QA mask values. Family loosely describes products from a satellite series or product type."""
-        return self.deployment.get('qa_mask', {}).get(family, {})
-
-
-class EasiAlias(collections.UserDict):
-    """Custom UserDict that returns a default measurement name for a given key if defined.
-    Else returns the key as the value. Items can not be set."""
-    def __init__(self, default:dict = {}):
-        self.data = default
-        self._log = _getlogger(self.__class__.__name__)
-    def __getitem__(self, key):
-        if key in self.data:
-            return self.data[key]
-        return key
-    def __setitem__(self, key, val):
-        self._log.error(f'Error <{self.__class__.__name__}>: Can not set items')
-
-
-class EasiCachingProxy():
-    """Set, unset and return information about the user's caching-proxy configuration"""
-    def __init__(self):
-        pass
-
-
-def _getlogger(name):
-    """Return a logger. Define here to limit external dependencies"""
-    # Default logger
-    #   log.hasHandlers() = False
-    #   log.getEffectiveLevel() = 30 = warning
-    #   log.propagate = True
-    logger = logging.getLogger(name)
-    logger.setLevel(logging.INFO)
-    if not len(logger.handlers):
-        logger.addHandler(logging.StreamHandler(sys.stdout))
-    logger.propagate = False  # Do not propagate up to root logger, which may have other handlers
-    return logger

easi_tools/load_s2l2a.py

@@ -1,293 +0,0 @@
-#!python
-
-# Sentinel-2 L2A Collection 0 scaling and offset corrections.
-# - Applies to data indexed from https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a
-# - The newer https://earth-search.aws.element84.com/v1/collections/sentinel-2-c1-l2a (Collection 1) may not be affected in the same way
-#
-# TL;DR:
-# DN values in COG files have different definitions depending on the processing baseline version
-# and whether the offset change has been pre-applied by the cloud data custodian.
-#
-# Background:
-#
-# ESA has undertaken a reprocessing of the Sentinel-2 L2A product that includes
-# a change to the offset value used to convert digital numbers (in file) to
-# scientific values (reflectances).
-#
-# https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/level-2a-algorithms-products
-#
-# L2A algorithm and products: Starting with the PB 04.00 (25th January 2022), the dynamic
-# range of the Level-2A products is shifted by a band-dependent constant: BOA_ADD_OFFSET.
-# This offset will allow encoding negative surface reflectances that may occur over very 
-# dark surfaces.
-#
-# L2A_SRi = (L2A_DNi + BOA_ADD_OFFSETi) / QUANTIFICATION_VALUEi
-#
-# QUANTIFICATION_VALUEi = 10000
-# BOA_ADD_OFFSETi = -1000
-#
-# refl = (dn -1000) / 10000
-# refl = dn/10000 - 1000/10000
-# refl = dn * 0.0001 - 0.1
-# 
-# These are the values in the EASI product definition, e.g.
-# https://explorer.asia.easi-eo.solutions/products/s2_l2a.odc-product.yaml
-#
-# Example workflow:
-#
-# ESA's reprocessing is flowing through to the AWS open data repository of S2 L2A but
-# while this stabilises we may see inconsistencies in time series queries due to:
-# - More than one processed version of a dataset (scene) in the AWS bucket and indexed in an EASI database
-# - Datasets (scenes) that indicate they have an offset applied by ESA but the offset correction
-#   has been not been applied to the COG
-#
-# Element-84 discussion:
-# https://github.com/Element84/earth-search/issues/23#issuecomment-1834674853
-
-
-import xarray as xr
-import pandas as pd
-import logging
-from pathlib import Path
-import sys, re
-
-import datacube
-from datacube.api.core import output_geobox
-from datacube.api.query import SPATIAL_KEYS, CRS_KEYS, OTHER_KEYS
-from datacube.utils import masking
-
-
-# Set logger
-log = logging.getLogger(Path(__file__).stem)
-log.setLevel(logging.INFO)
-log.addHandler(logging.StreamHandler(sys.stdout))
-
-# Constants
-search_keys = (
-    'product',
-    'time',
-    'geopolygon',
-    'like',
-    'limit',
-    'ensure_location',
-    'dataset_predicate',
-) + SPATIAL_KEYS + CRS_KEYS + OTHER_KEYS
-
-# TODO: get measurement aliases from the ODC product record
-refl_bands = {
-    'coastal','band_01','B01','coastal_aerosol',
-    'blue','band_02','B02',
-    'green','band_03','B03',
-    'red','band_04','B04',
-    'rededge1','band_05','B05','red_edge_1',
-    'rededge2','band_06','B06','red_edge_2',
-    'rededge3','band_07','B07','red_edge_3',
-    'nir','band_08','B08','nir_1',
-    'nir08','band_8a','B8A','nir_2',
-    'nir09','band_09','B09','nir_3',
-    'swir16','band_11','B11','swir_1','swir_16',
-    'swir22','band_12','B12','swir_2','swir_22',
-}
-scale_factor = 0.0001
-add_offset = -0.1
-
-
-def highest_sequence_number(matches: list) -> dict:
-    """Filter for the highest element84 processing sequence number per scene (scene label excluding the sequence number)
-    
-    : return : { scene_id_excluding_sequence_number : { highest_sequence_number : datacube.model.Dataset }}
-    """
-    p = re.compile(r'(S2.+)_([0-9]+)_(L2A)')
-    sorter = {}
-    for ds in matches:
-        # Separate the scene label from the sequence number
-        label = ds.metadata_doc['label']
-        m = p.match(label)
-        if not m:
-            log.warning(f'Dataset label does not match expected pattern: {label}')
-            continue
-        key = f'{m.group(1)}_{m.group(3)}'
-        seq = int(m.group(2))
-        # Retain the highest sequence number
-        if key in sorter:
-            if list(sorter[key])[0] < seq:
-                sorter[key] = {seq: ds}
-        else:
-            sorter[key] = {seq: ds}
-    return sorter
-
-
-def ds_requires_offset(ds: datacube.model.Dataset) -> bool:
-    """Return True if a dataset's metadata indicates that the offset correction should be applied"""
-    props = ds.metadata_doc['properties']
-    
-    # If baseline is less than '04.00' then offset correction does not apply
-    baseline = props.get('s2:processing_baseline', '0.0')
-    p = re.compile(r'(\d+)\.(\d+)')
-    m = p.match(baseline)
-    if not m:
-        log.warning(f'Dataset processing_baseline does not match expected pattern: {baseline}')
-        return None
-    if int(m.group(1)) < 4:
-        return False
-
-    # If the boa_offset_applied has been applied then offset correction is not required
-    boa_offset_applied = props.get('earthsearch:boa_offset_applied', False)
-    return not boa_offset_applied
-
-
-def apply_correction_to_data(ds: xr.Dataset, offset: float = 0) -> xr.Dataset:
-    """Apply the scale and offset correction to each reflectance band where there is valid data (not nodata)"""
-    refl_vars = [x for x in ds.data_vars if x in refl_bands]
-    mask = masking.valid_data_mask(ds[refl_vars])
-    # Save on a dask step?
-    if offset == 0:
-        ds[refl_vars] = ds[refl_vars].where(mask) * scale_factor
-    else:
-        ds[refl_vars] = ds[refl_vars].where(mask) * scale_factor + offset
-    return ds
-    
-
-def load_s2l2a_with_offset(
-    dc: datacube.Datacube,
-    query: dict,
-) -> xr.Dataset:
-    """
-    Replaces datacube.load(**query) for s2_l2a products.
-
-    Method:
-    - Find all datasets matching the query (dc.find_datasets)
-    - Filter for the highest element84 processing sequence number per scene (scene label excluding the sequence number)
-    - Filter into two lists for datasets that have
-      - "s2:processing_baseline" >= "04.00" and "earthsearch:boa_offset_applied" == False (offset correction required)
-      - everything else (no correction required)
-    - If either list is empty then load the non-empty list, apply scale (and offset if required), and return the xarray Dataset
-    - Load and combine the two lists of datasets
-      - Load each list, apply scale (and offset if required)
-      - Concat on time dimension and sort by time
-      - Return the combined xarray Dataset
-    
-    Notes:
-    - Any 'groupby' function is applied to each of the xarray Datasets prior to them being combined.
-      This could create "extra" (non-grouped) time layers in the combined Dataset if the groupby function
-      would have grouped datasets (scenes) from both lists.
-    - Scale and offset are applied to the reflectance bands where there is valid data (not `nodata`).
-      This includes applying the "scale_factor" even if no datasets require the offset correction.
-      Other masks can be applied by the user (e.g. pixel quality or cloud masking).
-    """
-
-    product = query.get('product', '<all products>')
-    if product != 's2_l2a':
-        log.error(f'This function only applies to the "s2_l2a" product, not: {product}')
-        return None
-    
-    # Find all datasets matching the query
-    matches = None
-    if 'datasets' in query:
-        matches = query['datasets']
-        del query['datasets']
-    if matches is None:
-        search_params = {k:v for k,v in query.items() if k in search_keys}
-        matches = dc.find_datasets(**search_params)
-    if 'skip_broken_datasets' not in query:
-        # This helps to avoid data loading error messages
-        query['skip_broken_datasets'] = True
-    
-    # Filter for the highest element84 processing sequence number
-    sorter = highest_sequence_number(matches)
-
-    # Filter into two lists
-    offset_applied, offset_required = [], []
-    for key in sorter.keys():
-        ds = list(sorter[key].values())[0]
-        isrequired = ds_requires_offset(ds)
-        if isrequired is None:
-            continue
-        elif isrequired:
-            offset_required.append(ds)
-        else:
-            offset_applied.append(ds)
-    matches_combined = offset_applied + offset_required
-
-    # If either list is empty then no separation and merge is required
-    this_offset = None
-    if len(offset_applied) == 0:
-        log.info('All datasets require offset correction')
-        msg = 'The valid_data_mask, scale and offset have been applied to the reflectance bands'
-        this_offset = add_offset
-    if len(offset_required) == 0:
-        log.info('No datasets require offset correction')
-        msg = 'The valid_data_mask and scale (no offset) have been applied to the reflectance bands'
-        this_offset = 0
-    if this_offset is not None:
-        data = dc.load(
-            datasets = matches_combined,
-            **query
-        )
-        xx = apply_correction_to_data(data, this_offset)
-        log.info(msg)
-        return xx
-
-    # DEBUG: What do we have
-    # def func(s):
-    #     p = re.compile('(\d{8})')
-    #     m = p.search(s[0])
-    #     if m:
-    #         return m.group(1)
-    # log.info(f'Number of datasets in initial query: {len(matches)}')
-    # log.info(f'{sorted([(x.metadata_doc["label"],x.id) for x in matches], key=func)}')
-    # log.info(f'Number of datasets with offset applied: {len(offset_applied)}')
-    # log.info(f'{sorted([(x.metadata_doc["label"],x.id) for x in offset_applied], key=func)}')
-    # log.info(f'Number of datasets without offset applied: {len(offset_required)}')
-    # log.info(f'{sorted( [(x.metadata_doc["label"],x.id) for x in offset_required], key=func)}')
-    # return
-
-    # Else, load data into two Datasets
-    log.info('Mix of datasets found with either offset required or not.')
-    log.info('We will load two xarrays, apply offset where required, and merge into one xarray.')
-
-    # 1. Ensure the target geobox covers all datasets
-    target_geobox = output_geobox(
-        datasets = matches_combined,
-        **query,
-    )
-    
-    # 2. Edit the query for our needs
-    # Ensure that dask time chunking = 1
-    dask_input = None
-    if 'dask_chunks' in query:
-        dask_input = query['dask_chunks']  # Save
-        if dask_input.get('time', 1) != 1:
-            query['dask_chunks'].update({'time': 1})
-    # Remove keys that are not compatible with 'like'
-    for x in ('output_crs', 'resolution', 'align'):
-        if x in query:
-            del query[x]
-    
-    # 3. Load two xarrays
-    data_offset_applied = dc.load(
-        datasets = offset_applied,
-        like = target_geobox,
-        **query
-    )
-    data_offset_required = dc.load(
-        datasets = offset_required,
-        like = target_geobox,
-        **query
-    )
-
-    # 4. Apply respective scale and offsets
-    data_offset_applied = apply_correction_to_data(data_offset_applied)
-    data_offset_required = apply_correction_to_data(data_offset_required, add_offset)
-    
-    # 5. Combine the two xarrays
-    combined = xr.concat([data_offset_applied, data_offset_required], dim='time')
-    combined = combined.sortby('time')
-    
-    # 6. Reapply any time > 1 chunking
-    if dask_input is not None:
-        if dask_input.get('time', 1) != 1:
-            combined = combined.chunk(dask_input)
-
-    log.info('The valid_data_mask, scale and offset have been applied to the reflectance bands')
-    return combined

easi_tools/notebook_utils.py

@@ -1,171 +0,0 @@
-#!python3
-
-# A collection of utilities that can be used in Python notebooks.
-#
-# License: Apache 2.0
-
-# Created for EASI Hub training notebooks, https://dev.azure.com/csiro-easi/easi-hub-public/_git/hub-notebooks
-
-# Data tools
-import numpy as np
-import xarray as xr
-import pandas as pd
-import geopandas as gpd
-import datacube
-from datacube.utils import masking
-from datetime import datetime
-
-# hvPlot, Holoviews, Datashader and Bokeh
-import hvplot.pandas
-import hvplot.xarray
-import panel as pn
-import holoviews as hv
-# hv.extension("bokeh", logo=False)  # Its likely set from in the notebooks
-
-# Jupyter Lab
-from IPython.display import HTML
-
-# Python
-import sys, os, re
-import logging
-from pathlib import Path
-from collections import Counter
-import contextlib
-
-# Dask
-import dask
-from dask.distributed import Client, LocalCluster
-from dask_gateway import Gateway
-
-# EASIDefaults
-from . import EasiDefaults
-
-# Set logger
-logger = logging.getLogger(Path(__file__).stem)
-logger.setLevel(logging.INFO)
-if not len(logger.handlers):
-    logger.addHandler(logging.StreamHandler(sys.stdout))
-
-
-def display_table(
-    df: pd.DataFrame,
-    panel: bool = False,
-):
-    """Display the full pandas dataframe. If panel is True use a panel object"""
-    table = None
-    if panel:
-        #  Dicts are rendered as "[object Object]". Need to set a formatter, I guess.
-        table = pn.widgets.DataFrame(df,
-            # sizing_mode='stretch_width',  # equal column widths, full screen
-            autosize_mode='fit_viewport',  # fitted columns, about 90-95% width
-            # reorderable=True,  # didn't work first try
-        )
-    else:
-        with pd.option_context("display.max_rows", None,
-                               "display.max_columns", None,
-                               "display.max_colwidth", -1):
-            table = HTML( df.to_html().replace(r"\n", "<br>") )
-    display(table)
-    
-        
-def heading(txt: str):
-    """Print a simple HTML heading"""
-    display(HTML( f"<h4>{txt}</h4>" ))
-    
-    
-def hv_table_hook(plot, element):
-    """Selected options for hv.table() formatting
-    
-    Use: df.hv.table().opts(hooks=[hv_table_hook])
-    """
-    plot.handles["table"].autosize_mode="fit_viewport"
-    # Other examples
-    # plot.handles['table'].row_height = 40
-    # from bokeh.models.widgets import DateFormatter
-    # plot.handles['table'].columns[6].formatter = DateFormatter(format='%Y-%m-%d')
-    
-
-def xarray_object_size(data):
-    """Return a formatted string"""
-    val, unit = data.nbytes / (1024 ** 2), "MB"
-    if val > 1024:
-        val, unit = data.nbytes / (1024 ** 3), "GB"
-    return f"Dataset size: {val:.2f} {unit}"
-
-
-def mostcommon_crs(dc, query):
-    """Adapted from https://github.com/GeoscienceAustralia/dea-notebooks/blob/develop/Tools/dea_tools/datahandling.py"""
-    matching_datasets = dc.find_datasets(**query)
-    crs_list = [str(i.crs) for i in matching_datasets]
-    crs_mostcommon = None
-    if len(crs_list) > 0:
-        # Identify most common CRS
-        crs_counts = Counter(crs_list)
-        crs_mostcommon = crs_counts.most_common(1)[0][0]
-    else:
-        logger.warning("No data was found for the supplied product query")
-    return crs_mostcommon
-
-
-def initialize_dask(use_gateway=False, workers=(1,2), wait=False, local_port=8786, **kwargs):
-    """Initialize a Dask Gateway or Local cluster"""
-    # Check inputs
-    if isinstance(workers, (int, float)):
-        workers = (int(workers), int(workers))
-    if len(workers) != 2:
-        logger.error("Require workers to be a single integer or a 2-element tuple/list")
-        return None, None
-    if isinstance(local_port, (str, float)):
-        local_port = int(local_port)
-    
-    # Dask gateway
-    if use_gateway:
-        gateway = Gateway()
-        clusters = gateway.list_clusters()
-        if not clusters:
-            logger.info("Starting new cluster")
-            cluster = gateway.new_cluster(**kwargs)
-        else:
-            logger.info(f"An existing cluster was found. Connecting to: {clusters[0].name}")
-            cluster = gateway.connect(clusters[0].name)
-        client = cluster.get_client()
-        cluster.adapt(minimum=workers[0], maximum=workers[1])
-        if wait:
-            logger.info("Waiting for at least one cluster worker")
-            # client.wait_for_workers(n_workers=1)  # Before release 2023.10.0
-            client.sync(client._wait_for_workers,n_workers=1) # Since release 2023.10.0
-
-    # Local cluster
-    else:
-        cluster = LocalCluster(n_workers=4)
-        client = Client(cluster)
-        server = f'https://hub.{EasiDefaults().domain}'  # Or replace if not using EasiDefaults
-        user = os.environ.get('JUPYTERHUB_SERVICE_PREFIX')  # Current user
-        dask.config.set({"distributed.dashboard.link": f'{server}{user}' + "proxy/{port}/status"})  # port is evaluated by dask
-        
-    return cluster, client
-
-
-def localcluster_dashboard(client, server="https://hub.csiro.easi-eo.solutions"):
-    """Return a dashboard link using jupyter proxy"""
-    dashboard_link = client.dashboard_link
-    for host in ("127.0.0.1", "localhost"):
-        if host in dashboard_link:
-            port = re.search(r":(\d+)\/status", dashboard_link).group(1)
-            dashboard_link = f'{server}{os.environ["JUPYTERHUB_SERVICE_PREFIX"]}proxy/{port}/status'
-            break
-    return dashboard_link
-
-
-@contextlib.contextmanager
-def unset_cachingproxy():
-    """Unset the EASI caching proxy with a context manager"""
-    # Inspired by https://stackoverflow.com/a/34333710
-    env = os.environ
-    remove = ("AWS_HTTPS", "GDAL_HTTP_PROXY")
-    update_after = {k: env[k] for k in remove}
-    try:
-        [env.pop(k, None) for k in remove]
-        yield
-    finally:
-        env.update(update_after)

new_import_ODC.py

@@ -81,187 +81,48 @@ from sklearn.metrics import mean_squared_error, r2_score
 import joblib
 
 
-def load_data_from_rasterio(dc, date_range, longtitude_range, latitude_range):
-    """
-    Load Sentinel-2 L2A data directly from S3 COGs using rasterio.
-    Returns a xarray Dataset with 10980x10980 resolution data.
-    
-    This approach:
-    - Loads ALL available data without spatial filtering
-    - Uses direct S3 COG access (rasterio) for reliability
-    - Returns data at native 10m resolution
-    - Matches the pipeline's downstream processing requirements
-    """
-    
-    print(f'Loading Sentinel-2 data from S3 COGs (rasterio)...')
-    print(f'  Date range: {date_range}')
-    print(f'  Target area: Lon {longtitude_range}, Lat {latitude_range}')
-    
-    try:
-        # Get first matching scene
-        datasets = list(dc.find_datasets(
-            product='s2_l2a',
-            time=date_range
-        ))
-        
-        if not datasets:
-            print(f'❌ No datasets found for date range {date_range}')
-            return None
-            
-        selected = datasets[0]
-        print(f'\n📦 Using scene: {selected.metadata.label}')
-        
-        # Load measurements from S3 COGs
-        measurements_to_load = ['red', 'green', 'blue', 'nir', 'scl']
-        data_dict = {}
-        
-        print(f'\n⏳ Loading bands from S3 COGs...')
-        for band_name in measurements_to_load:
-            if band_name in selected.measurements:
-                band_path = selected.measurements[band_name]['path']
-                
-                try:
-                    with rasterio.open(band_path) as src:
-                        data = src.read(1)
-                        data_dict[band_name] = data
-                        print(f'   ✅ {band_name}: {data.shape}, dtype={data.dtype}')
-                except Exception as e:
-                    print(f'   ⚠️ Could not load {band_name}: {e}')
-        
-        if not data_dict:
-            print('❌ Could not load any bands')
-            return None
-        
-        # Create xarray Dataset
-        print(f'\n🔄 Converting to xarray Dataset...')
-        
-        # Get dimensions from red band (highest resolution)
-        red_data = data_dict['red']
-        y_size, x_size = red_data.shape
-        
-        # Create coordinate arrays (placeholder - real georeferencing would come from rasterio metadata)
-        y_coords = np.arange(y_size)
-        x_coords = np.arange(x_size)
-        
-        # Create data arrays for each variable
-        data_vars = {}
-        for band_name, band_data in data_dict.items():
-            if band_data.shape == red_data.shape:
-                # Same resolution - direct assignment
-                data_vars[band_name] = (['y', 'x'], band_data)
-            else:
-                # Different resolution (e.g., SCL at 20m) - resample to match red
-                from scipy import ndimage
-                scale_factor = red_data.shape[0] // band_data.shape[0]
-                resampled = ndimage.zoom(band_data, scale_factor, order=0)
-                data_vars[band_name] = (['y', 'x'], resampled)
-        
-        # Create xarray Dataset
-        data = xr.Dataset(
-            data_vars,
-            coords={
-                'x': x_coords,
-                'y': y_coords
-            }
-        )
-        
-        print(f'\n✅ Data converted successfully!')
-        print(f'   Dimensions: {dict(data.sizes)}')
-        print(f'   Variables: {list(data.data_vars)}')
-        print(f'   Shape: {red_data.shape}')
-        print(f'   Data type: numpy arrays (in-memory)')
-        
-        return data
-        
-    except Exception as e:
-        print(f'❌ Error loading data: {e}')
-        import traceback
-        traceback.print_exc()
-        return None
-
-
 def load_data(dc, date_range, longtitude_range, latitude_range):
-    """
-    Load Sentinel-2 L2A data using direct datacube.load() 
-    without spatial filtering (which was causing 0 results).
-    
-    Note: Data is loaded in UTM (EPSG:32648) to avoid CRS issues.
-    Spatial filtering on lat/lon is skipped to return maximum data.
-    """
     product = 's2_l2a'
-    native_crs = 'EPSG:32648'  # UTM Zone 48N for Vietnam
+    query = {
+        'product': product,                     # Product name
+        'x': longtitude_range,    # "x" axis bounds
+        'y': latitude_range,      # "y" axis bounds
+        'time': date_range,           # Any parsable date strings
+    }
+    native_crs = notebook_utils.mostcommon_crs(dc, query)
+    print(f'Most common native CRS: {native_crs}')
     measurements = ['red', 'nir', 'scl']
-    
+
-    print(f'Loading Sentinel-2 data (EPSG:32648)...')
+    load_params = {
-    print(f'  Time range: {date_range}')
+        'measurements': measurements,                   # Selected measurement or alias names
-    print(f'  Measurements: {measurements}')
+        'output_crs': native_crs,                       # Target EPSG code
-    
+        'resolution': (-10, 10),                        # Target resolution
-    try:
+        'group_by': 'solar_day',                        # Scene grouping
-        # Load ALL available data WITHOUT dask_chunks (forces immediate load)
+        'dask_chunks': {'x': 2048, 'y': 2048},          # Dask chunks
-        # This avoids the metadata issue with dc.load() when using dask_chunks
+    }
-        data = dc.load(
+    data = load_s2l2a_with_offset(
-            product=product,
+        dc,
-            time=date_range,
+        query | load_params   # Combine the two dicts that contain our search and load parameters
-            measurements=measurements,
+    )
-            output_crs=native_crs,
+    return data
-            resolution=(-10, 10),
-            group_by='solar_day',
-            skip_broken_datasets=True
-        )
-        
-        print(f'✅ Data loaded successfully!')
-        print(f'   Dimensions: {dict(data.sizes)}')
-        print(f'   Time steps: {len(data.time)}')
-        print(f'   Spatial extent: x={len(data.x)}, y={len(data.y)}')
-        print(f'   Data type: numpy arrays (not Dask)')
-        
-        return data
-        
-    except Exception as e:
-        print(f'❌ Error loading data: {e}')
-        import traceback
-        traceback.print_exc()
-        return None
 
 
 def mask_clean(data):
-    """
+    flag_name = 'scl'
-    Clean data by masking clouds and bad pixels using the SCL (Scene Classification Layer).
+    flag_desc = masking.describe_variable_flags(data[flag_name])  # Pandas dataframe
-    
+    display(flag_desc)
-    SCL classes:
+    display(flag_desc.loc['qa'].values[1])
-    - 0: No Data
+    # Create a "data quality" Mask layer
-    - 1: Saturated/Defective
+    flags_def = flag_desc.loc['qa'].values[1]
-    - 2: Dark Area Pixels
+    good_pixel_flags = [flags_def[str(i)] for i in [2, 4, 5, 6]]  # To pass strings to enum_to_bool()
-    - 3: Cloud Shadows
+
-    - 4: Vegetation ✓ GOOD
+    # enum_to_bool calculates the pixel-wise "or" of each set of pixels given by good_pixel_flags
-    - 5: Not Vegetated ✓ GOOD
+    # 1 = good data
-    - 6: Water ✓ GOOD
+    # 0 = "bad" data
-    - 7: Unclassified ✓ GOOD
+    good_pixel_mask = enum_to_bool(data[flag_name], good_pixel_flags)
-    - 8: Cloud Medium Probability ✗ BAD
-    - 9: Cloud High Probability ✗ BAD
-    - 10: Thin Cirrus ✗ BAD
-    - 11: Snow/Ice ✗ BAD
-    """
-    
-    # Good pixel classes (keep these)
-    good_pixel_classes = [4, 5, 6, 7]
-    
-    # Create mask: 1 where SCL is in good_pixel_classes, 0 otherwise
-    good_pixel_mask = data['scl'].isin(good_pixel_classes)
-    
-    print(f'✅ Cloud masking applied')
-    print(f'   Good pixel classes: {good_pixel_classes}')
-    print(f'   Mask created (dask-backed, not yet computed)')
-    
-    # Get all variables except SCL
     data_layer_names = [x for x in data.data_vars if x != 'scl']
-    
+    # Apply good pixel mask to blue, green, red and nir.
-    # Apply mask to all layers
     result = data[data_layer_names].where(good_pixel_mask).persist()
-    
-    print(f'   Data variables masked: {data_layer_names}')
-    print(f'   Result persisted to workers')
-    
     return result
 
 
@@ -499,17 +360,7 @@ def load_data_sen2(dc, date_range, coordinates):
         'y': latitude_range,      # "y" axis bounds
         'time': date_range,           # Any parsable date strings
     }
-    
+    native_crs = notebook_utils.mostcommon_crs(dc, query)
-    # Try to get native CRS, default to EPSG:32648 (UTM Zone 48N) for Vietnam
-    try:
-        native_crs = notebook_utils.mostcommon_crs(dc, query)
-        if native_crs is None:
-            print('⚠️ Could not determine native CRS, using EPSG:32648 (UTM Zone 48N)')
-            native_crs = 'EPSG:32648'
-    except Exception as e:
-        print(f'⚠️ Error determining CRS: {e}, using EPSG:32648')
-        native_crs = 'EPSG:32648'
-    
     print(f'Most common native CRS: {native_crs}')
     
     # measurements = ['red','green', 'blue', 'nir', 'scl']
@@ -522,27 +373,10 @@ def load_data_sen2(dc, date_range, coordinates):
         'group_by': 'solar_day',                        # Scene grouping
         'dask_chunks': {'x': 2048, 'y': 2048},          # Dask chunks
     }
-    
+    data = load_s2l2a_with_offset(
-    try:
+        dc,
-        data = load_s2l2a_with_offset(
+        query | load_params   # Combine the two dicts that contain our search and load parameters
-            dc,
+    )
-            query | load_params   # Combine the two dicts that contain our search and load parameters
-        )
-    except Exception as e:
-        print(f'❌ Error loading data: {e}')
-        print('Attempting direct dc.load without offset correction...')
-        data = dc.load(
-            product=product,
-            x=longtitude_range,
-            y=latitude_range,
-            time=date_range,
-            measurements=measurements,
-            output_crs=native_crs,
-            resolution=(-10, 10),
-            group_by='solar_day',
-            dask_chunks={'x': 2048, 'y': 2048},
-            skip_broken_datasets=True
-        )
     return data
 
 def mask_cloud(data):