Commit 3220fa1c authored by Eva Börgens's avatar Eva Börgens
Browse files

add functionality to use non constant grid std per time step

parent 809d1080
......@@ -28,7 +28,9 @@ else:
timeseries = None
if flag_uncertainty or flag_matrix:
results = compute_covariance(region_coords, grid_std, flag_uncertainty, flag_matrix)
std_vec = vectorize_std(grid_std, lon, lat, region_coords)
results = compute_covariance(region_coords, std_vec, flag_uncertainty, flag_matrix)
else:
results = None
......
......@@ -8,8 +8,10 @@ import math
from typing import Dict
import scipy.special as sp
import scipy.interpolate as interp
import numpy as np
rho = np.pi/180
R = 6367.4447
......@@ -42,7 +44,7 @@ def get_grid_area(lon : np.ndarray, lat : np.ndarray) -> np.ndarray:
else:
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)
......@@ -208,7 +210,7 @@ def yaglom(dist: float,
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):
def distance(lon_0: float, lat_0: float, lon_1: float, lat_1: float) -> float:
"""
convert geograpic coordinates to spherical distances
......@@ -247,7 +249,7 @@ def distance(lon_0: float, lat_0: float, lon_1: float, lat_1: float):
return d_i
def azimut_angle(lon_0: float, lat_0: float, lon_1: float, lat_1: float):
def azimut_angle(lon_0: float, lat_0: float, lon_1: float, lat_1: float) -> float:
"""
get azimut angle between geograpic coordinates
......@@ -312,7 +314,7 @@ def compute_covariance(region_coords: np.ndarray,
lat = np.array([r[1] for r in region_coords])
n = len(lat)
m = len(gridstd)
m = gridstd.shape[0]
# Grid weights
area = get_grid_area(lon, lat)
......@@ -340,7 +342,7 @@ def compute_covariance(region_coords: np.ndarray,
w_2 = area[ii] / region_size
for j in range(m):
sigma = gridstd[j]
sigma = gridstd[j,ii]
if np.isnan(sigma):
continue
if flag_uncertainty:
......@@ -356,14 +358,15 @@ def compute_covariance(region_coords: np.ndarray,
w_2 = area[jj] / region_size
for j in range(m):
sigma = gridstd[j]
if np.isnan(sigma):
sigma_ii = gridstd[j, ii]
sigma_jj = gridstd[j, jj]
if np.isnan(sigma_ii) or np.isnan(sigma_jj):
continue
if flag_uncertainty:
var_model[j] += w_1 * w_2 * sigma**2 * corr
var_model[j] += w_1 * w_2 * sigma_ii * sigma_jj * corr
if flag_matrix:
cov_model[j, ii, jj] = sigma ** 2 * corr
cov_model[j, jj, ii] = sigma ** 2 * corr
cov_model[j, ii, jj] = sigma_ii * sigma_jj * corr
cov_model[j, jj, ii] = sigma_ii * sigma_jj * corr
result = {}
if flag_uncertainty:
std_model = np.sqrt(var_model)
......@@ -373,7 +376,8 @@ def compute_covariance(region_coords: np.ndarray,
return result
def get_timeseries(grid, lon_grid, lat_grid, region_coords):
def get_timeseries(grid: np.ndarray, lon_grid: np.ndarray,
lat_grid: np.ndarray, region_coords: np.ndarray) -> np.ndarray:
"""
Returns mean tws time series of region
......@@ -415,3 +419,42 @@ def get_timeseries(grid, lon_grid, lat_grid, region_coords):
timeseries[i] = np.nansum(area * grid_part[region_id_lat, region_id_lon]) / np.sum(area)
return timeseries, True
def vectorize_std(grid_std: np.ndarray, lon: np.ndarray,
lat: np.ndarray, region_coords: np.ndarray) -> np.ndarray:
"""
Returns vectorized grid point std needed for function "compute_covariance"
Parameters
----------
grid_std: np.ndarray
tws std grid, size [t,n,m]
lon: np.ndarray
longitude of grid, size [m]
lat: np.ndarray
latitude of grid, size [n]
region_coords: np.ndarray
coordinates of region of interest, size [l,2]
Returns
-------
np.ndarray
size [t,l]
"""
t = grid_std.shape[0]
l = region_coords.shape[0]
Lon_grid, Lat_grid = np.meshgrid(lon, lat)
std_vec = np.zeros(t,l)
for i in range(t):
if len(np.unique(grid_std[i][~np.isnan(grid_std[i])]))==1:
std_vec[i,:] = np.unique(grid_std[i][~np.isnan(grid_std[i])])
else:
f = interp.interp2d(Lon_grid.ravel(), Lat_grid.ravel(),grid_std[i].ravel() )
std_vec[i,:] = f(region_coords[:,0], region_coords[:,1])
return std_vec
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment