Commit 8d72a548 by Eva Börgens

remove typos from docstring and export Docstrings as documentation

parent 29e2767f
 src.covariance API documentation

Module src.covariance

Module handling the computation of the covariances and uncertainties

Expand source code
"""
Module handling the computation of the covariances and uncertainties
"""

import math
from typing import Dict

import scipy.special as sp
import numpy as np

rho = np.pi/180
R = 6367.4447

def get_grid_area(lon : np.ndarray, lat : np.ndarray) -> np.ndarray:
"""
Function getting the area weights of a coordinate list

Args:
lon : np.ndarray
array of longitudes
lat : np.ndarray
array of latitudes

Returns:
np.ndarray
area per grid point
"""
area = np.zeros(len(lat))
delta_lon = np.abs(np.unique(lon)[1]-np.unique(lon)[0])
delta_lat = np.abs(np.unique(lat)[1] - np.unique(lat)[0])

for i, lat_i in enumerate(lat):
area[i] = get_area(0, lat_i-delta_lat/2, delta_lon, lat_i+delta_lat/2)

return area

def get_area(lon_ll: float, lat_ll: float, lon_ur: float,lat_ur: float) -> float:
"""
returns area of rectangular on sphere, defined by corner points lower left ll and upper right ur

Args:
lon_ll: float
lon of lower left corner
lat_ll: float
lat of lower left corner
lon_ur: float
lon of upper right corner
lat_ur: float
lat of upper right corner

Returns:
float
area of grid point
"""
area = np.abs((lon_ll-lon_ur)*rho)*np.abs((np.sin(lat_ll*rho)-np.sin(lat_ur*rho)))*R**2
return area

def cov_function_lat(lat1: float, lat2: float,
dist: float,
theta: float,
ny: float,
a_0: float, k_a_0_2: float, k_a_0_4: float, k_a_0_6: float, k_a_0_8: float,
a_1: float, k_a_1_2: float, k_a_1_4: float, k_a_1_6: float, k_a_1_8: float,
c_0: float, k_2: float, k_4: float, k_6: float, k_8: float) -> float:
"""
Function to compute covariance function between two points according to
publication Boergens et al. 2020

Args:
lat1, lat2: float
Latitude of the two points
d: float
distance between points
theta: float
azimuth angle between points
ny: float
order of Bessle function
a0: float
anisotropy parameter
ka0_2, ka0_4, ka0_6, ka0_8: float
Legende polynome Args for a0
a1: float
isotropy shape parameter
ka1_2, ka1_4, ka1_6, ka1_8: float
Legende polynome Args for a1
c0: float
amplitude parameter
k2,k4, k6, k8: float
Legende polynome Args for c0
Returns:
float
Covariance
"""

k = [1, 0, k_2, 0, k_4, 0, k_6, 0, k_8]
P_1 = sum_legendre(len(k) - 1, k, lat1 * rho)
P_2 = sum_legendre(len(k) - 1, k, lat2 * rho)

k_a_0 = [1, 0, k_a_0_2, 0, k_a_0_4, 0, k_a_0_6, 0, k_a_0_8]
P_a_0_1 = sum_legendre(len(k_a_0) - 1, k_a_0, np.mean([lat1, lat2]) * rho)

A_0 = P_a_0_1 * a_0

k_a_1 = [1, 0, k_a_1_2, 0, k_a_1_4, 0, k_a_1_6, 0, k_a_1_8]
P_a_1_1 = sum_legendre(len(k_a_1) - 1, k_a_1, np.mean([lat1, lat2]) * rho)

A_1 = P_a_1_1 * a_1
if dist==0:
dist = 0.001
C_0 = yaglom(dist, theta, A_0, A_1, c_0, ny)

Cov = P_1*P_2*C_0

return Cov

def legendre_polynome(n: int, lat: float) -> float:
"""
Computes  Legendre Polynome of degree n at given latitude lat

Args:
n: int
lat: float

Returns:
float
"""
sin_lat = np.sin(lat)
P_n = sp.legendre(n)
P_n_x = P_n(sin_lat)
return P_n_x

def sum_legendre(n_max: int, leg_weights: np.ndarray, lat:float) -> float:
"""
Computes weighted sum of Legendre Polynomes

Args:
n_max: int
maximum degree of Legendre Polynomes
leg_weights: np.ndarray
arrays of the weights for the sum
lat: float
latitude where the Legendre Polynomes are summed

Returns:
float
"""
if n_max!=len(leg_weights)-1:
raise('N_max coefficients are needed ')
p_sum = 0
for i in range(n_max+1):
p_sum += leg_weights[i] * legendre_polynome(i, lat)
return p_sum

def yaglom(dist: float,
theta: float,
a_0: float,
a_1: float,
c_0: float,
ny: float) -> float:
"""
Function to compute the adapted Yaglom function

Args:
dist: float
spherical distance
theta: float
azimut angel
a_0: float
anisotropic width parameter
a_1: float
isotropic width parameter
c_0: float
global scaling factor
ny: int
Order of Bessel function

Returns:
float

"""
alpha = a_0 * np.sin(theta) + a_1

cov = c_0 / (alpha * dist) ** ny * sp.jv(ny, alpha * dist)
return cov

def distance(lon_0: float, lat_0: float, lon_1: float, lat_1: float):
"""
convert geograpic coordinates to spherical distances

Args:
lon_0: float
[degree]
lat_0: float
[degree]
lon_1: float
[degree]
lat_1: float
[degree]

Returns:
float
"""

hilf = sin_lat*sin_lat0 + cos_lat*cos_lat0*cos_lon
d_i = R * math.acos(hilf)

if lon_0==lon_1 and lat_0==lat_1:
d_i =0

return d_i

def azimut_angle(lon_0: float, lat_0: float, lon_1: float, lat_1: float):
"""
get azimut angle between geograpic coordinates

Args:
lon_0: float
[degree]
lat_0: float
[degree]
lon_1: float
[degree]
lat_1: float
[degree]

Returns:
float
"""

alpha = math.atan2(sin_lon,(cos_lat0*tan_lat-sin_lat0*cos_lon))
if lat_1 == 90:
alpha=0
elif lat_1 ==-90:
alpha = np.pi
return alpha

def compute_covariance(region_coords: np.ndarray,
gridstd: np.ndarray,
flag_uncertainty: bool, flag_matrix: bool) -> Dict[str, np.ndarray]:
"""
Function to compute the covariances for a region

Args:
region_coords: np.ndarray
coordinates of region, size [n,2]
gridstd: np.ndarray
standard deviation for each grid point, size[n]
flag_uncertainty: bool
flag, return uncertainty of mean tws of region
flag_matrix: bool
flag, return covariance matrix of region

Returns:
Dict[str, np.ndarray]
"""

lon = np.array([r[0] for r in region_coords])
lat = np.array([r[1] for r in region_coords])

n = len(lat)
m = len(gridstd)

# Grid weights
area = get_grid_area(lon, lat)
region_size = np.sum(area)

# Parameter of covariance function
x = np.array([1.09e-03, -2.92, -1.90, -0.86, -4.60,
2.66e-03, 0.40, -0.30, 0.07, 0.32,
3.74, -0.2, -0.27, -0.17, -0.22])

if flag_uncertainty:
var_model = np.zeros(m)
else:
var_model = None
if flag_matrix:
cov_model = np.zeros((m,n,n))
else:
cov_model = None

for ii in range(n):
dist = 0
theta = np.pi * 2
corr = cov_function_lat(lat[ii], lat[ii], dist, theta, 2, *x)
w_1 = area[ii] / region_size
w_2 = area[ii] / region_size

for j in range(m):
sigma = gridstd[j]
if np.isnan(sigma):
continue
if flag_uncertainty:
var_model[j] += w_1 * w_2 * sigma **2 * corr
if flag_matrix:
cov_model[j,ii,ii] = sigma **2 * corr

for jj in range(ii+1, n):
dist = distance(lon[ii], lat[ii], lon[jj], lat[jj])
theta = azimut_angle(lon[ii], lat[ii], lon[jj], lat[jj])
corr = cov_function_lat(lat[ii], lat[jj], dist, theta, 2, *x)
w_1 = area[ii] / region_size
w_2 = area[jj] / region_size

for j in range(m):
sigma = gridstd[j]
if np.isnan(sigma):
continue
if flag_uncertainty:
var_model[j] += w_1 * w_2 * sigma**2 * corr
if flag_matrix:
cov_model[j, ii, jj] = sigma ** 2 * corr
cov_model[j, jj, ii] = sigma ** 2 * corr
result = {}
if flag_uncertainty:
std_model = np.sqrt(var_model)
result['uncertainty'] = std_model
if flag_matrix:
result['matrix'] = cov_model
return result

def get_timeseries(grid, lon_grid, lat_grid, region_coords):
"""
Returns mean tws time series of region

Args:
grid: np.ndarray
tws grid, size [t,n,m]
lon_grid: np.ndarray
longitude of grid, size [m]
lat_grid: np.ndarray
latitude of grid, size [n]
region_coords: np.ndarray
coordinates of region of interest, size [l,2]

Returns:
np.ndarray
size [t]
"""
lon_region = np.array([r[0] for r in region_coords])
lat_region = np.array([r[1] for r in region_coords])

if len(np.isin(lon_grid, lon_region).nonzero()[0])< len(np.unique(lon_region)) or \
len(np.isin(lat_grid, lat_region).nonzero()[0])< len(np.unique(lat_region)):
print("Warning: Region coordinates are not located on the given TWS grid, thus no "
"mean tws time series can be computed. Continue without timeseries output")
return None, False

# Grid weights
area = get_grid_area(lon_region, lat_region)

region_id_lat = np.array([np.nonzero(lat_grid==l)[0] for l in lat_region]).squeeze()
region_id_lon = np.array([np.nonzero(lon_grid==l)[0] for l in lon_region]).squeeze()

timeseries = np.zeros(grid.shape[0])
for i in range(grid.shape[0]):
grid_part = grid[i,:,:]
timeseries[i] = np.nansum(area * grid_part[region_id_lat, region_id_lon]) / np.sum(area)

return timeseries, True

Functions

def azimut_angle(lon_0: float, lat_0: float, lon_1: float, lat_1: float)

get azimut angle between geograpic coordinates

Args: lon_0: float [degree] lat_0: float [degree] lon_1: float [degree] lat_1: float [degree]

Returns

float

Expand source code
def azimut_angle(lon_0: float, lat_0: float, lon_1: float, lat_1: float):
"""
get azimut angle between geograpic coordinates

Args:
lon_0: float
[degree]
lat_0: float
[degree]
lon_1: float
[degree]
lat_1: float
[degree]

Returns:
float
"""

alpha = math.atan2(sin_lon,(cos_lat0*tan_lat-sin_lat0*cos_lon))
if lat_1 == 90:
alpha=0
elif lat_1 ==-90:
alpha = np.pi
return alpha
def compute_covariance(region_coords: numpy.ndarray, gridstd: numpy.ndarray, flag_uncertainty: bool, flag_matrix: bool) ‑> Dict[str, numpy.ndarray]

Function to compute the covariances for a region

Args: region_coords: np.ndarray coordinates of region, size [n,2] gridstd: np.ndarray standard deviation for each grid point, size[n] flag_uncertainty: bool flag, return uncertainty of mean tws of region flag_matrix: bool flag, return covariance matrix of region

Returns

Dict[str, np.ndarray]

Expand source code
def compute_covariance(region_coords: np.ndarray,
gridstd: np.ndarray,
flag_uncertainty: bool, flag_matrix: bool) -> Dict[str, np.ndarray]:
"""
Function to compute the covariances for a region

Args:
region_coords: np.ndarray
coordinates of region, size [n,2]
gridstd: np.ndarray
standard deviation for each grid point, size[n]
flag_uncertainty: bool
flag, return uncertainty of mean tws of region
flag_matrix: bool
flag, return covariance matrix of region

Returns:
Dict[str, np.ndarray]
"""

lon = np.array([r[0] for r in region_coords])
lat = np.array([r[1] for r in region_coords])

n = len(lat)
m = len(gridstd)

# Grid weights
area = get_grid_area(lon, lat)
region_size = np.sum(area)

# Parameter of covariance function
x = np.array([1.09e-03, -2.92, -1.90, -0.86, -4.60,
2.66e-03, 0.40, -0.30, 0.07, 0.32,
3.74, -0.2, -0.27, -0.17, -0.22])

if flag_uncertainty:
var_model = np.zeros(m)
else:
var_model = None
if flag_matrix:
cov_model = np.zeros((m,n,n))
else:
cov_model = None

for ii in range(n):
dist = 0
theta = np.pi * 2
corr = cov_function_lat(lat[ii], lat[ii], dist, theta, 2, *x)
w_1 = area[ii] / region_size
w_2 = area[ii] / region_size

for j in range(m):
sigma = gridstd[j]
if np.isnan(sigma):
continue
if flag_uncertainty:
var_model[j] += w_1 * w_2 * sigma **2 * corr
if flag_matrix:
cov_model[j,ii,ii] = sigma **2 * corr

for jj in range(ii+1, n):
dist = distance(lon[ii], lat[ii], lon[jj], lat[jj])
theta = azimut_angle(lon[ii], lat[ii], lon[jj], lat[jj])
corr = cov_function_lat(lat[ii], lat[jj], dist, theta, 2, *x)
w_1 = area[ii] / region_size
w_2 = area[jj] / region_size

for j in range(m):
sigma = gridstd[j]
if np.isnan(sigma):
continue
if flag_uncertainty:
var_model[j] += w_1 * w_2 * sigma**2 * corr
if flag_matrix:
cov_model[j, ii, jj] = sigma ** 2 * corr
cov_model[j, jj, ii] = sigma ** 2 * corr
result = {}
if flag_uncertainty:
std_model = np.sqrt(var_model)
result['uncertainty'] = std_model
if flag_matrix:
result['matrix'] = cov_model
return result
def cov_function_lat(lat1: float, lat2: float, dist: float, theta: float, ny: float, a_0: float, k_a_0_2: float, k_a_0_4: float, k_a_0_6: float, k_a_0_8: float, a_1: float, k_a_1_2: float, k_a_1_4: float, k_a_1_6: float, k_a_1_8: float, c_0: float, k_2: float, k_4: float, k_6: float, k_8: float) ‑> float

Function to compute covariance function between two points according to publication Boergens et al. 2020

Args

lat1, lat2: float
Latitude of the two points
d
float distance between points
theta
float azimuth angle between points
ny
float order of Bessle function
a0
float anisotropy parameter
ka0_2, ka0_4, ka0_6, ka0_8: float
Legende polynome Args for a0
a1
float isotropy shape parameter
ka1_2, ka1_4, ka1_6, ka1_8: float
Legende polynome Args for a1
c0
float amplitude parameter

k2,k4, k6, k8: float