Merge pull request #14 from jejjohnson/main

[WIP] MODIS Pipeline

rs_tools/_src/data/modis/bands.py

@@ -3,4 +3,10 @@
         "EV_500_Aggr1km_RefSB":[3,4,5,6,7],
         "EV_1KM_RefSB": [8,9,10,11,12,"13lo","13hi","14lo","14hi",15,16,17,18,19,26],
         "EV_1KM_Emissive": [20,21,22,23,24,25,27,28,29,30,31,32,33,34,35,36],
-}
+}
+
+
+MODIS_CHANNEL_NUMBERS = [
+    str(i) for i in range(1, 39)
+] + ['13lo', '14lo', '13hi', '14hi']
+}

scripts/pipeline/preprocess_modis.py

@@ -1,16 +1,18 @@
 import autoroot
 import numpy as np
 
+import rioxarray
 from pathlib import Path
 from dataclasses import dataclass
 from typing import Optional, List, Union, Tuple
 
-from rs_tools import goes_download, modis_download, MODIS_VARIABLES
+from rs_tools import goes_download, modis_download, MODIS_VARIABLES, MODIS_CHANNEL_NUMBERS
 from rs_tools._src.utils.io import get_list_filenames
 from rs_tools._src.geoprocessing.grid import create_latlon_grid
 import typer
 from loguru import logger
 import xarray as xr
+from satpy import Scene
 
 
 @dataclass
@@ -27,20 +29,61 @@ class GeoProcessingParams:
 
 @dataclass
 class MODISGeoProcessing:
-    resolution: float = 0.25 # in degrees?
+    # resolution: float = 0.25 # in degrees?
+    resolution: float = 1_000 # km
     save_path: str = "./"
+    channels: list[str] = MODIS_CHANNEL_NUMBERS
 
-    def geoprocess(self, params: GeoProcessingParams, files: List[str]):
-        save_path: str = params.save_path.join(self.save_path)
-        target_grid: np.ndarray = create_latlon_grid(params.region, self.resolution) # create using np.linspace?
+    # def geoprocess_files(self, params: GeoProcessingParams, files: List[str]):
+    # resolution: float = 0.25 # in degrees?
+    # save_path: str = "./"
 
-        # loop through files
-        # open dataset
-        # stack variables to channels
-        # resample
-        # convert units (before or after resampling???)
-        # save as netcdf
+    # def geoprocess(self, params: GeoProcessingParams, files: List[str]):
+    #     save_path: str = params.save_path.join(self.save_path)
+    #     target_grid: np.ndarray = create_latlon_grid(params.region, self.resolution) # create using np.linspace?
 
+    #     # loop through files
+    #     # open dataset
+    #     # stack variables to channels
+    #     # resample
+    #     # convert units (before or after resampling???)
+    #     # save as netcdf
+
+    def geoprocess_file(self, file: str, save_name: str="test") -> None:
+
+        # load MODIS SCENE
+        ds, time_stamp = self.load_modis_scene(file)
+        # convert units
+        ds = self.convert_units(file)
+        # TODO: resample
+        # save to raster
+        time_stamp = time_stamp.strftime("%Y%m%d%H%M")
+        save_name =f"modis_{time_stamp}.tif"
+        ds.rio.to_raster(Path(self.save_path).joinpath(save_name))
+
+        return None
+
+    def load_modis_scene(self, file: str) -> xr.Dataset:
+
+        # load modis scene
+        scn = Scene(reader="modis_l1b", filenames=[file])
+
+        # load channels
+        scn.load(self.channels, resolution = self.resolution)
+
+        time_stamp = scn.start_time + (scn.end_time - scn.start_time) / 2
+
+        # convert to xarray
+        ds = scn.to_xarray_dataset()
+
+        return ds, time_stamp
+
+    def convert_units(self, ds: xr.Dataset) -> xr.Dataset:
+
+        ds.rio.write_crs("epsg:4326", inplace=True)
+        ds = ds.set_coords(["latitude", "longitude"])
+        ds = ds.assign_coords({"latitude": ds.latitude, "longitude": ds.longitude})
+        return ds