diff --git a/rs_tools/_src/data/modis/bands.py b/rs_tools/_src/data/modis/bands.py index 2d24c39..d4e22c8 100644 --- a/rs_tools/_src/data/modis/bands.py +++ b/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], -} \ No newline at end of file +} + + +MODIS_CHANNEL_NUMBERS = [ + str(i) for i in range(1, 39) +] + ['13lo', '14lo', '13hi', '14hi'] +} diff --git a/scripts/pipeline/preprocess_modis.py b/scripts/pipeline/preprocess_modis.py index 2f8d803..b3fdc17 100644 --- a/scripts/pipeline/preprocess_modis.py +++ b/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