47 lines
1.5 KiB
Python
47 lines
1.5 KiB
Python
import numpy as np
|
|
import xarray as xr
|
|
import os
|
|
|
|
# Define default auxiliary data file (OLCI OCP ADF)
|
|
ref_path = os.path.dirname(os.path.realpath(__file__))
|
|
# ADF_OCP = os.path.join(ref_path, '..', 'AuxiliaryData/OCP/S3A_OL_2_OCP_AX_20160216T000000_20991231T235959_20240327T100000___________________EUM_O_AL_008.SEN3/OL_2_OCP_AX.nc')
|
|
ADF_OCP = os.path.join(ref_path, 'BRDF_LUTs','BRDF_%s.nc')
|
|
|
|
def solve_2nd_order_poly(A, B, C):
|
|
""" Solve 2nd order polynomial inversion
|
|
where coefficients are xr dataArray
|
|
Take only positive solution, otherwise provide 0.
|
|
"""
|
|
|
|
# Compute solution according to sign of delta
|
|
delta = B*B - 4*A*C
|
|
mask = delta > 0
|
|
|
|
# By default (and when delta < 0), take value at extremum
|
|
x = - B / (2*A)
|
|
|
|
# When delta > 0, take biggest solutions
|
|
x_1 = (-B.where(mask) - np.sqrt(delta.where(mask))) / (2*A.where(mask))
|
|
x_2 = (-B.where(mask) + np.sqrt(delta.where(mask))) / (2*A.where(mask))
|
|
x_sol = xr.where(x_1 > x_2, x_1, x_2)
|
|
x = xr.where(delta > 0, x_sol, x)
|
|
|
|
# Take only positive value
|
|
x = xr.where(x > 0, x, 0.)
|
|
|
|
return x
|
|
|
|
def drop_unused_coords(var):
|
|
for coord in var.coords:
|
|
if coord not in var.dims:
|
|
var = var.drop(coord)
|
|
return var
|
|
|
|
def squeeze_trivial_dims(ds,dimOrderAccordingTo='Rw'):
|
|
squeezedDims = {}
|
|
for d0,dim in enumerate(ds[dimOrderAccordingTo].dims):
|
|
if ds.sizes[dim] == 1:
|
|
squeezedDims[dim] = d0
|
|
|
|
return ds.squeeze(), squeezedDims
|