CIMR L2 TSA Top-Level Script

# Stand-alone script to run TSA algorithm and write L2 TSA product, alternatively call through:
# run_CIMR_L2_PolarOceans.py --no-oza-adjust "C:/Users/zschende/OneDrive - Ilmatieteen laitos/Files/Projects/CIMR Devalgo/Data/SCEPS/SCEPS_l1b_sceps_geo_polar_scene_1_unfiltered_tot_minimal_nom_nedt_apc_tot_v2p1.nc"

import os 
import sys
import numpy as np
import xarray as xr
import pyresample as pr
import importlib

# local imports
if '/tools/' not in sys.path:
    sys.path.insert(0, os.path.abspath('../.../') + 'algorithm/tools/')
import l2_format as l2
import TSA_algorithm as algorithm
import l2_tools as tools

# reload local imports
importlib.reload(l2)
importlib.reload(algorithm)
importlib.reload(tools)

# tag 'parameters' for the CLI with papermill, adjust if necessary for stand-alone run
l1b_path = ''
l1x_path = os.path.abspath('../../../..../') + 'Data/SCEPS/SCEPS_l1x@KA_sceps_geo_polar_scene_1_unfiltered_tot_minimal_nom_nedt_apc_tot_v2p1.nc'
aux_dir = '../data/auxiliary/'
l2_dir = os.path.abspath('../../../..../') + 'Data/L2 Files/'
l2_grid = 'ease2-3.125km-nh'

# check input parameters
l1x_scenes = ('devalgo_geometric', 'devalgo_radiometric', 'sceps_polar1')

if not os.path.isfile(l1x_path):
    if l1x_path not in l1x_scenes:
        raise ValueError("The input L1X file does not exist, and is not one of the pre-registered {}".format(l1x_scenes,))

if not os.path.isdir(l2_dir):
    raise ValueError("The L2 output directory {} does not exist.".format(l2_dir))

if not os.path.isdir(aux_dir):
    raise ValueError("The auxiliary directory {} does not exist.".format(aux_dir))

# Handle pre-defined L1X files
if l1x_path in l1x_scenes:
    test_card_name = l1x_path
    if l1x_path == 'devalgo_geometric':    
        # DEVALGO simulated geometric test card
        l1x_path = os.path.abspath('../../../..../') + 'Data/L1C Files/Geometric/'
        l1x_fn = 'W_PT-DME-Lisbon-SAT-CIMR-1X@KA_C_DME_20230417T105425_LD_20280110T114800_20280110T115700_TN.nc'
    elif l1x_path == 'devalgo_radiometric':
        # DEVALGO simulated radiometric test card
        l1x_path = os.path.abspath('../../../..../') + 'Data/L1C Files/Radiometric/'
        l1x_fn = 'W_PT-DME-Lisbon-SAT-CIMR-1X@KA_C_DME_20230420T103323_LD_20280110T114800_20280110T115700_TN.nc'
    elif l1x_path == 'sceps_polar1':
        # SCEPS simulated radiometric test card
        l1x_path = os.path.abspath('../../../..../') + 'Data/SCEPS/'
        l1x_fn = 'SCEPS_l1x@KA_sceps_geo_polar_scene_1_unfiltered_tot_minimal_nom_nedt_apc_tot_v2p1.nc'
    else:
        raise ValueError("Unknown test_card {}".format(l1x_path))
        
    l1x_path = os.path.join(l1x_path, l1x_fn)
else:
    if 'devalgo_test_scene_1' in os.path.basename(l1x_path):
        test_card_name = 'devalgo_radiometric'
    elif 'devalgo_test_scene_2' in os.path.basename(l1x_path):
        test_card_name = 'devalgo_geometric'
    elif 'sceps_geo_polar_scene_1':
        test_card_name = 'sceps_polar1'
    else:
        test_card_name = 'unknown'

area_def = pr.load_area('ease2_adef.yaml',l2_grid)

# load landmask if available
if test_card_name == 'devalgo_geometric':
    dem_ref_scenario = os.path.abspath('../../../..../') + 'Data/Test_scenes_downscaled_projected/test_scene_2_compressed_lowres.nc'
    ds = xr.open_dataset(dem_ref_scenario)
    lons = ds['Longitude'].values
    lats = ds['Latitude'].values

    swath_def = pr.geometry.SwathDefinition(lons=lons, lats=lats)
    land_mask = pr.kd_tree.resample_nearest(swath_def, np.float64(ds['landflag'].values), area_def,
                            radius_of_influence=20000/4, fill_value=np.nan)

elif test_card_name == 'devalgo_radiometric':
    dem_ref_scenario = os.path.abspath('../../../..../') + 'Data/Test_scenes_downscaled_projected/test_scene_1_compressed_lowres.nc'
    ds= xr.open_dataset(dem_ref_scenario)
    lons = ds['Longitude'].values
    lats = ds['Latitude'].values

    swath_def = pr.geometry.SwathDefinition(lons=lons, lats=lats)
    land_mask = pr.kd_tree.resample_nearest(swath_def, np.float64(ds['landflag'].values), area_def,
                            radius_of_influence=20000/4, fill_value=np.nan)

elif test_card_name == 'sceps_polar1':
    dem_ref_scenario = os.path.abspath('../../../..../') + 'Data/SCEPS/cimr_sceps_geo_card_devalgo_polarscene_1_20161217_harmonised_v2p0_surface.nc'
    ds = xr.open_dataset(dem_ref_scenario)
    lons = ds['longitude'].values[0,:,:]
    lats = ds['latitude'].values[0,:,:]

    swath_def = pr.geometry.SwathDefinition(lons=lons, lats=lats)
    land_mask = pr.kd_tree.resample_nearest(swath_def, np.float64(ds['land_sea_ice_mask'].values[0,:,:]), area_def,
                            radius_of_influence=20000/4, fill_value=np.nan)

    land_mask[(land_mask == 1) | (land_mask == 9)] = 0      # set water (incl. sea ice)
    land_mask[land_mask == 2] = 1                           # set land

ds.close()

algo = 'Pulliainen2010'
algo_version = '0.9.0'

# read L1X data
data_fwd,geo_fwd,data_bck,geo_bck = tools.read_l1x(l1x_path)

# detection forward and backward
TSA_fwd = algorithm.dry_snow_detection(data_fwd,tsa_algorithm=algo)
TSA_bck = algorithm.dry_snow_detection(data_bck,tsa_algorithm=algo)

C:\Users\zschende\OneDrive - Ilmatieteen laitos\Files\Projects\CIMR Devalgo\ATBD\TerrestrialSnowArea_ATBD_v2\algorithm\TSA_algorithm.py:27: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  shp = tuple(data['KA'].dims[d] for d in ['n_scans_interleave_feed','n_samples_earth'])

# individual reprojection
TSA_fwd_proj = tools.reproject_to_grid(TSA_fwd,geo_fwd,area_def=area_def)
TSA_bck_proj = tools.reproject_to_grid(TSA_bck,geo_bck,area_def=area_def)

# intermediate combined data
TSA_comb, TSA_comb_uncert = tools.combine_proj(TSA_fwd_proj,TSA_bck_proj)

# status_flag
TSA_status_flag = np.zeros(TSA_comb.shape, dtype='i8')

TSA_status_flag[TSA_status_flag == 0] = 8                    # 8: no data, out of grid
TSA_status_flag[land_mask == 0] = 0                          # 0: water
TSA_status_flag[(TSA_comb == 0) & (land_mask == 1)] = 1      # 1: land
TSA_status_flag[(TSA_comb == 1) & (land_mask == 1)] = 2      # 2: valid snow
# TSA_status_flag[] = 3                                      # 3: [placeholder]

status_flag_comment = '0: water; 1: land; 2: snow (valid); 3-7: [placeholders]; 8: no data, out of grid.'

# TSA_uncertainty
TSA_uncertainty = np.zeros(TSA_comb.shape, dtype='i8') 

TSA_uncertainty[TSA_comb_uncert == 0] = 0                    # 0: very likely snow-free
TSA_uncertainty[TSA_comb_uncert == 0.5] = 1                  # 1: likely snow
TSA_uncertainty[TSA_comb_uncert == 1] = 2                    # 2: very likely snow

# TSA
TSA_l2 = TSA_comb
# TSA_l2 = np.zeros(TSA_status_flag.shape)

# TSA_l2[TSA_l2 == 0] = np.nan                                 # nan: no data, out of grid
# TSA_l2[TSA_comb == 0] = 0                                    # 0: snow free
# TSA_l2[TSA_comb == 1] = 1                                    # 1: snow covered (dry)
# TSA_l2[TSA_comb == 2] = 2                                    # 2: snow covered (wet) [placeholder]

# get template L2 format (netCDF/CF) from the Tools module
ds_l2 = l2.get_CIMR_L2_template('grid', geo_def=area_def, add_time=None)

# create DataArray for TSA from template
da_tsa = xr.DataArray(TSA_l2, coords=ds_l2['template'].coords, dims=ds_l2['template'].dims,
                        attrs=ds_l2['template'].attrs, name='tsa')
da_tsa.attrs['long_name'] = 'Terrestrial Snow Area ({} algorithm)'.format(algo)
da_tsa.attrs['standard_name'] = 'terrestrial_snow_area'
da_tsa.attrs['units'] = 1
da_tsa.attrs['coverage_content_type'] = 'NA'
da_tsa.attrs['auxiliary_variables'] = 'terrestrial_snow_area_uncertainty,status_flag'

# create DataArray for TSA_uncertainty from template
da_uncert = xr.DataArray(TSA_uncertainty, coords=ds_l2['template'].coords, dims=ds_l2['template'].dims,
                        attrs=ds_l2['template'].attrs, name='tsa_uncertainty')
da_uncert.attrs['long_name'] = 'Qualitative uncertainty for Terrestrial Snow Area'
da_uncert.attrs['standard_name'] = 'terrestrial_snow_area_uncertainty'
da_uncert.attrs['coverage_content_type'] = 'QualityInformation'
da_uncert.attrs['units'] = 1

# create DataArray for status_flag from template
da_flag = xr.DataArray(TSA_status_flag, coords=ds_l2['template'].coords, dims=ds_l2['template'].dims,
                        attrs=ds_l2['template'].attrs, name='status_flag')
da_flag.attrs['long_name'] = 'Status flag for Terrestrial Snow Area'
da_flag.attrs['coverage_content_type'] = 'AuxiliaryInformation'
da_flag.attrs['comment'] = status_flag_comment

# add data arrays to ds_l2 object
ds_l2 = ds_l2.merge(da_tsa)
ds_l2 = ds_l2.merge(da_uncert)
ds_l2 = ds_l2.merge(da_flag)

# customize global attributes
ds_l2.attrs['title'] = 'CIMR L2 NRT3H Terrestrial Snow Area'
ds_l2.attrs['summary'] = 'Terrestrial Snow Area computed with the prototype algorithm developed in the ESA CIMR DEVALGO study. The algorithm combines brightness temperatures from the Ku and Ka channels. The product file contains the TSA, its uncertainties, and status flag.'
ds_l2.attrs['l1b_file'] = l1b_path
ds_l2.attrs['l1x_file'] = l1x_path
ds_l2.attrs['algorithm_version'] = algo_version
ds_l2.attrs['creator_name'] = 'Lina Zschenderlein'
ds_l2.attrs['creator_email'] = 'lina.zschenderlein@fmi.fi'
ds_l2.attrs['institution'] = 'Finnish Meteorological Institute'

# remove 'template' variable
ds_l2 = ds_l2.drop('template')

# write to file
l2_n = 'cimr_devalgo_l2_tsa_{}_{}.nc'.format(l2_grid, test_card_name.replace('_','-'), )
l2_n = os.path.join(l2_dir, l2_n)
ds_l2.to_netcdf(l2_n, format='NETCDF4_CLASSIC')
print(l2_n)