#!/usr/bin/python
#Code: verif_modules.py (Python Module Library)
#Author: Nate Snook (CASA/CAPS/SoM)
#Written: Dec. 2009
#Last modified: 20 Oct. 2015
#---------------------------------------------------------------------------------#
# List of contents:
# define_kernel(float, float) -- defines a kernel suitable for NEP or other verification methods
# define_kernel_ellipse(float, float, float) -- defines elliptical (rather than circular) kernels
# NEP_2D(float[:,:], float, float[:,:]) -- calculates neighborhood ensemble probability of a
# given 2D (x-y) field using a specified kernel to perform a convolution (SciPy)
# pm_mean_2D(float[:,:,:]) -- calculates the probability-matched mean of a field given an ensemble
# of 2D slices [nens, ny, nx] of that field
#---------------------------------------------------------------------------------#
[docs]def define_kernel(radius, dropoff=0.0):
"""
Defines a 2D kernel suitable for use in neighborhood ensemble probability calculations
(or other verification methods requiring the use of a 2D convolution).
Args:
radius: The radius of the kernel, in gridpoints (can be any whole number)
dropoff: Set this to a value between 0.0 and 1.0 to have points near the edge of the radius recieve a lower weight. By default, dropoff is set to 0.0, giving all points within the radius equal weight.
Returns:
Kernel: a 2D array containing the kernel needed for performing convolutions
"""
from sys import exit
from numpy import zeros, floor
#Check for valid positive radius
if radius < 0:
exit('ERROR: Invalid neighborhood radius -- radius cannot be negative!')
kernel = zeros(((2*floor(radius))+1, 2*(floor(radius))+1 ))
for i in range(0,kernel.shape[0]):
for j in range(0,kernel.shape[1]):
icoord = abs(i - (floor(radius)))
jcoord = abs(j - (floor(radius)))
if ((icoord*icoord) + (jcoord*jcoord) <= (radius*radius)):
kernel[i, j] = 1.0 - ( ((icoord*icoord) + (jcoord*jcoord)) / (radius*radius) ) * dropoff
#Return the valid kernel array
return kernel
#-----------------------------------------------------------------------------------#
[docs]def define_kernel_ellipse(xradius, yradius, dropoff=0.0):
"""
Defines an elliptical 2D kernel suitable for use in neighborhood ensemble probability
calculations (or other verification methods requiring the use of a 2D convolution).
Args:
xradius: The semi-major (or semi-minor) axis in the east-west direction, in gridpoints (can be any whole number)
yradius: The semi-major (or semi-minor) axis in the north-south direction, in gridpoints (can be any whole number)
dropoff: Set this to a value between 0.0 and 1.0 to have points near the edge of the radius recieve a lower weight. By default, dropoff is set to 0.0, giving all points within the radius equal weight.
Returns:
Kernel: a 2D array containing the kernel needed for performing convolutions
"""
from sys import exit
from numpy import zeros, array, shape, floor
#Check for valid positive radius
if xradius < 0 or yradius < 0:
exit('ERROR: Invalid ellipse settings -- semi-minor/semi-major axes cannot be negative!')
kernel = zeros(((2*floor(yradius))+1, 2*(floor(xradius))+1))
for i in range(0,kernel.shape[0]):
for j in range(0,kernel.shape[1]):
icoord = abs(i - (floor(yradius)))
jcoord = abs(j - (floor(xradius)))
if ( ((icoord*icoord) / (yradius*yradius)) + ((jcoord*jcoord) / (xradius*xradius)) <= 1.0 ):
kernel[i, j] = 1.0 - ( ((icoord*icoord)/(yradius*yradius)) + ((jcoord*jcoord)/(xradius*xradius)) ) * dropoff
#Return the valid kernel array
return kernel
#------------------------------------------------------------------#
[docs]def NEP_2D(var2D_ens, threshold, kernel):
"""
A function for performing neighborhood ensemble probability (NEP)
Given a ensemble of 2D slices (var2D_ens), this function will perform the
NEP of P[var2D > threshold] and return a corresponding 2D array with the
NEP values. The NEP will use a neighborhood defined by kernel. The kernel
is a 2D array and is defined using the define_kernel function (also
contained in this library).
Args:
var2D_ens: The 2D ensemble of x-y slices to calculate NEP for [n_ens, ny, nx]
threshold: NEP will be calculated for P[var2D_ens > threshold] within the neighborhood where threshold can be any real number.
kernel: a 2D array of some size < (ny, nx) with values between 0 and 1 to serve as the kernel for the 2D convolution.
Returns:
var_NEP: a 2D [ny, nx] field containing NEP for P[var2D_ens > threshold] for the given neighborhood (defined by kernel)
"""
import numpy as np
from scipy.signal import convolve2d
neighborhood_ensemble = np.empty(var2D_ens.shape)
for ens_mem in xrange(var2D_ens.shape[0]):
neighborhood_ensemble[ens_mem] = convolve2d(var2D_ens[ens_mem] > threshold, kernel, mode='same', boundary='symm')
var_NEP = neighborhood_ensemble.sum(axis=0) / (neighborhood_ensemble.shape[0] * kernel.sum())
return var_NEP
#------------------------------------------------------------------#
[docs]def ens_prob_within_dist(var2D_ens, threshold, kernel):
"""
A function for calculating ensemble probability of the form "probability of
X within Y km of a point". Similar to, but distinct from, NEP_2D.
Given a ensemble of 2D slices (var2D_ens), this function will compute the ensemble
probability of P[var2D > threshold (anywhere within kernel)] and return a corresponding
2D array with the probability values. The kernel is a 2D array and is defined using
the define_kernel function (also contained in this library).
Args:
var2D_ens: The 2D ensemble of x-y slices to calculate NEP for [n_ens, ny, nx]
threshold: NEP will be calculated for P[var2D_ens > threshold] within the neighborhood where threshold can be any real number.
kernel: a 2D array of some size < (ny, nx) with values of 0 or 1 to serve as the kernel for the 2D convolution.
Returns:
ens_prob: a 2D [ny, nx] field containing NEP for P[var2D_ens > threshold] for the given neighborhood (defined by kernel)
"""
import numpy as np
from scipy.signal import convolve2d
neighborhood_ensemble = np.empty(var2D_ens.shape)
# for ens_mem in np.arange(1, var2D_ens.shape[0], 1):
for ens_mem in xrange(var2D_ens.shape[0]):
neighborhood_ensemble[ens_mem] = convolve2d(var2D_ens[ens_mem] > threshold, kernel, mode='same', boundary='symm')
neighborhood_ensemble = np.where(neighborhood_ensemble > 0.0, 1.0, 0.0)
ens_prob = neighborhood_ensemble.sum(axis=0) / (neighborhood_ensemble.shape[0])
return ens_prob
#--------------------------------------------------------------------#
[docs]def pm_mean_2D(var2D_ens):
"""
A function for computing probability-matched mean (PM mean)
Given an ensemble of 2D slices (var2D_ens), this function will perform the
PM mean calculation and return the PM mean as a 2D array over the same slice.
var2D_ens should have dimensions of (n_ens, ny, nx).
Args:
var2D_ens: an ensemble of 2D slices of a variable field [nens, ny, nx]
Returns:
pm_mean: the probability-matched mean of var2D_ens
"""
import numpy as np
n_ens = var2D_ens.shape[0]
ny = var2D_ens.shape[1]
nx = var2D_ens.shape[2]
ensmean = np.mean(var2D_ens, axis=0) #Calculate ensemble mean
all_values = []
for jy in range(0, ny):
for ix in range(0, nx):
for member in range(0, n_ens):
all_values.append(var2D_ens[member, jy, ix])
#Construct a tethered array to track where our ranked values are.
tethered_array = np.zeros((ny, nx))
for jy in range(0, ny):
for ix in range(0, nx):
tethered_array[jy, ix] = jy + (ny * ix)
#Construct 1-D arrays to rank values
ensmean_1d = []
tethered_1d = []
for jy in range(0, ny):
for ix in range(0, nx):
ensmean_1d.append(ensmean[jy, ix])
tethered_1d.append(tethered_array[jy, ix])
#Prepare the all_values array by sorting it from high to low:
all_values.sort()
all_values.reverse()
#Tethering array elements to their position as tuples ensures the sort is reversible.
position_tuples = zip(tethered_1d, ensmean_1d)
#The lambda performs a "mini-function" -- google "python lambda notation" for info.
sorted_tuples = sorted(position_tuples, key = lambda x: x[1])
sorted_tuples.reverse() #Put largest values first
sorted_tuples = zip(*sorted_tuples)
tethered_1d = sorted_tuples[0]
ensmean_1d = sorted_tuples[1]
pm_mean = np.zeros((ny, nx))
for point in range(0, nx*ny, 1):
#The tethered array stores information on the position of the ranked points in ensmean.
#We must first convert this data back into i- and j-locations, then assign the corresponding data.
loc_jrank = np.mod(tethered_1d[point], ny)
loc_irank = np.floor(int(tethered_1d[point] / ny))
pm_mean[loc_jrank, loc_irank] = all_values[0 + (n_ens * point)]
return pm_mean
#--------------------------------------------------------------------#