Omaha #3085 - moved common.derivparam.python, common.inventory, common.wxmath to core

Former-commit-id: 786f3c846b23ff201d32f552ed240f1e97acace4
This commit is contained in:
Steve Harris 2014-04-28 14:41:41 -05:00
parent 491b768204
commit a3618cb714
173 changed files with 0 additions and 14230 deletions

View file

@ -1,7 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<classpath>
<classpathentry kind="con" path="org.eclipse.jdt.launching.JRE_CONTAINER/org.eclipse.jdt.internal.debug.ui.launcher.StandardVMType/JavaSE-1.6"/>
<classpathentry kind="con" path="org.eclipse.pde.core.requiredPlugins"/>
<classpathentry kind="src" path="src"/>
<classpathentry kind="output" path="bin"/>
</classpath>

View file

@ -1,34 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<projectDescription>
<name>com.raytheon.uf.common.derivparam.python</name>
<comment></comment>
<projects>
</projects>
<buildSpec>
<buildCommand>
<name>org.python.pydev.PyDevBuilder</name>
<arguments>
</arguments>
</buildCommand>
<buildCommand>
<name>org.eclipse.jdt.core.javabuilder</name>
<arguments>
</arguments>
</buildCommand>
<buildCommand>
<name>org.eclipse.pde.ManifestBuilder</name>
<arguments>
</arguments>
</buildCommand>
<buildCommand>
<name>org.eclipse.pde.SchemaBuilder</name>
<arguments>
</arguments>
</buildCommand>
</buildSpec>
<natures>
<nature>org.eclipse.pde.PluginNature</nature>
<nature>org.eclipse.jdt.core.javanature</nature>
<nature>org.python.pydev.pythonNature</nature>
</natures>
</projectDescription>

View file

@ -1,7 +0,0 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<?eclipse-pydev version="1.0"?>
<pydev_project>
<pydev_property name="org.python.pydev.PYTHON_PROJECT_INTERPRETER">Default</pydev_property>
<pydev_property name="org.python.pydev.PYTHON_PROJECT_VERSION">python 2.7</pydev_property>
</pydev_project>

View file

@ -1,8 +0,0 @@
#Thu Dec 16 13:23:53 CST 2010
eclipse.preferences.version=1
org.eclipse.jdt.core.compiler.codegen.inlineJsrBytecode=enabled
org.eclipse.jdt.core.compiler.codegen.targetPlatform=1.6
org.eclipse.jdt.core.compiler.compliance=1.6
org.eclipse.jdt.core.compiler.problem.assertIdentifier=error
org.eclipse.jdt.core.compiler.problem.enumIdentifier=error
org.eclipse.jdt.core.compiler.source=1.6

View file

@ -1,19 +0,0 @@
Manifest-Version: 1.0
Bundle-ManifestVersion: 2
Bundle-Name: Derived Parameter Python Plug-in
Bundle-SymbolicName: com.raytheon.uf.common.derivparam.python
Bundle-Version: 1.14.0.qualifier
Bundle-Vendor: RAYTHEON
Bundle-RequiredExecutionEnvironment: JavaSE-1.6
Bundle-ActivationPolicy: lazy
Require-Bundle: com.raytheon.uf.common.python,
com.raytheon.uf.common.python.concurrent,
com.raytheon.uf.common.wxmath,
com.raytheon.uf.common.util,
com.raytheon.uf.common.datastorage,
com.raytheon.uf.common.derivparam;bundle-version="1.12.1174",
com.raytheon.uf.common.status;bundle-version="1.12.1174",
com.raytheon.uf.common.localization;bundle-version="1.14.0"
Export-Package: com.raytheon.uf.common.derivparam.python,
com.raytheon.uf.common.derivparam.python.function
Import-Package: com.raytheon.uf.common.inventory.tree

View file

@ -1,5 +0,0 @@
source.. = src/
output.. = bin/
bin.includes = META-INF/,\
.,\
res/

View file

@ -1,14 +0,0 @@
<beans xmlns="http://www.springframework.org/schema/beans"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://www.springframework.org/schema/beans http://www.springframework.org/schema/beans/spring-beans-3.1.xsd">
<bean id="pythonDerParFunctionAdapter"
class="com.raytheon.uf.common.derivparam.python.DerivParamPythonFunctionAdapter" />
<bean id="registeredPythonDerParFunctionAdapter"
class="com.raytheon.uf.common.derivparam.library.DerivedParameterGenerator"
factory-method="addFunctionAdapter">
<constructor-arg ref="pythonDerParFunctionAdapter" />
</bean>
</beans>

View file

@ -1,184 +0,0 @@
/**
* This software was developed and / or modified by Raytheon Company,
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
*
* U.S. EXPORT CONTROLLED TECHNICAL DATA
* This software product contains export-restricted data whose
* export/transfer/disclosure is restricted by U.S. law. Dissemination
* to non-U.S. persons whether in the United States or abroad requires
* an export license or other authorization.
*
* Contractor Name: Raytheon Company
* Contractor Address: 6825 Pine Street, Suite 340
* Mail Stop B8
* Omaha, NE 68106
* 402.291.0100
*
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
* further licensing information.
**/
package com.raytheon.uf.common.derivparam.python;
import java.io.File;
import java.io.IOException;
import java.text.SimpleDateFormat;
import java.util.Date;
import java.util.List;
import java.util.concurrent.ExecutionException;
import com.raytheon.uf.common.datastorage.records.IDataRecord;
import com.raytheon.uf.common.derivparam.DerivParamFunctionType.FunctionArgument;
import com.raytheon.uf.common.derivparam.IDerivParamFunctionAdapter;
import com.raytheon.uf.common.derivparam.library.DerivedParameterGenerator;
import com.raytheon.uf.common.localization.PathManagerFactory;
import com.raytheon.uf.common.python.concurrent.PythonJobCoordinator;
import com.raytheon.uf.common.status.IUFStatusHandler;
import com.raytheon.uf.common.status.UFStatus;
import com.raytheon.uf.common.status.UFStatus.Priority;
import com.raytheon.uf.common.util.FileUtil;
/**
* Python derived parameter adapter
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------ ---------- ----------- --------------------------
* Dec 16, 2010 mschenke Initial creation
* Jun 04, 2013 2041 bsteffen Switch derived parameters to use
* concurrent python for threading.
* Jan 30, 2014 #2725 ekladstrup Add name and extention get methods
* after removing RCP extention point
*
* </pre>
*
* @author mschenke
* @version 1.0
*/
public class DerivParamPythonFunctionAdapter implements
IDerivParamFunctionAdapter {
private static final transient IUFStatusHandler statusHandler = UFStatus
.getHandler(DerivParamPythonFunctionAdapter.class);
private static final String PYTHON = "python";
private static final String TEMPLATE_FILE = DerivedParameterGenerator.DERIV_PARAM_DIR
+ File.separator + PYTHON + File.separator + "functionTemplate.txt";
private PythonJobCoordinator<MasterDerivScript> coordinator;
/*
* (non-Javadoc)
*
* @see
* com.raytheon.uf.viz.derivparam.IDerivParamFunctionAdapter#createNewFunction
* (java.lang.String,
* com.raytheon.uf.viz.derivparam.DerivParamFunctionType.FunctionArgument[])
*/
@Override
public String createNewFunction(String functionName,
FunctionArgument[] arguments) {
File template = PathManagerFactory.getPathManager().getStaticFile(
TEMPLATE_FILE);
try {
String templateText = new String(FileUtil.file2bytes(template));
SimpleDateFormat dateFormat = new SimpleDateFormat("MMM dd, yyyy");
String date = dateFormat.format(new Date());
String argList = "";
boolean first = true;
for (FunctionArgument arg : arguments) {
if (!first) {
argList += ",";
}
first = false;
argList += arg.name;
}
return String.format(templateText, date, argList);
} catch (IOException e) {
statusHandler.handle(Priority.PROBLEM,
"Error reading in python template file");
}
return null;
}
/*
* (non-Javadoc)
*
* @see
* com.raytheon.uf.viz.derivparam.IDerivParamFunctionAdapter#getArgumentTypes
* ()
*/
@Override
public String[] getArgumentTypes() {
return new String[] { "Parameter" };
}
/*
* (non-Javadoc)
*
* @see com.raytheon.uf.viz.derivparam.IDerivParamFunctionAdapter#init()
*/
@Override
public void init() {
if (coordinator != null) {
coordinator.shutdown();
}
coordinator = PythonJobCoordinator
.newInstance(new MasterDerivScriptFactory());
}
/*
* (non-Javadoc)
*
* @see
* com.raytheon.uf.viz.derivparam.IDerivParamFunctionAdapter#executeFunction
* (java.lang.String, java.util.List)
*/
@Override
public List<IDataRecord> executeFunction(String name, List<Object> arguments)
throws ExecutionException {
try {
return coordinator.submitSyncJob(new MasterDerivScriptExecutor(
name, arguments));
} catch (InterruptedException e) {
throw new ExecutionException(e);
}
}
/*
* (non-Javadoc)
*
* @see com.raytheon.uf.viz.derivparam.IDerivParamFunctionAdapter#shutdown()
*/
@Override
public void shutdown() {
if (coordinator != null) {
coordinator.shutdown();
}
}
/*
* (non-Javadoc)
*
* @see com.raytheon.uf.viz.derivparam.IDerivParamFunctionAdapter#getName()
*/
@Override
public String getName() {
return "Python";
}
/*
* (non-Javadoc)
*
* @see
* com.raytheon.uf.viz.derivparam.IDerivParamFunctionAdapter#getExtension()
*/
@Override
public String getExtension() {
return "py";
}
}

View file

@ -1,450 +0,0 @@
/**
* This software was developed and / or modified by Raytheon Company,
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
*
* U.S. EXPORT CONTROLLED TECHNICAL DATA
* This software product contains export-restricted data whose
* export/transfer/disclosure is restricted by U.S. law. Dissemination
* to non-U.S. persons whether in the United States or abroad requires
* an export license or other authorization.
*
* Contractor Name: Raytheon Company
* Contractor Address: 6825 Pine Street, Suite 340
* Mail Stop B8
* Omaha, NE 68106
* 402.291.0100
*
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
* further licensing information.
**/
package com.raytheon.uf.common.derivparam.python;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import jep.JepException;
import com.raytheon.uf.common.datastorage.records.ByteDataRecord;
import com.raytheon.uf.common.datastorage.records.FloatDataRecord;
import com.raytheon.uf.common.datastorage.records.IDataRecord;
import com.raytheon.uf.common.datastorage.records.IntegerDataRecord;
import com.raytheon.uf.common.datastorage.records.LongDataRecord;
import com.raytheon.uf.common.datastorage.records.ShortDataRecord;
import com.raytheon.uf.common.datastorage.records.StringDataRecord;
import com.raytheon.uf.common.derivparam.library.DerivedParameterRequest;
import com.raytheon.uf.common.inventory.tree.CubeLevel;
import com.raytheon.uf.common.python.PythonInterpreter;
/**
* A script for running the master derived parameter script, which can run any
* of the derived parameter scripts
*
* <pre>
* SOFTWARE HISTORY
* Date Ticket# Engineer Description
* ------------- -------- ----------- --------------------------
* Jul 08, 2008 njensen Initial creation
* Nov 20, 2009 3387 jelkins Use derived script's variableId instead of filename
* Nov 21, 2009 3576 rjpeter Refactored to populate DerivParamDesc.
* Oct 29, 2013 2476 njensen Renamed numeric methods to numpy
* Apr 11, 2014 2947 bsteffen Allow returning NaN
*
* </pre>
*
* @author njensen
* @version 1.0
*/
public class MasterDerivScript extends PythonInterpreter {
protected static final String RESULT = "__result";
private static final String DATA_NAME = "Data";
private final Map<Object, List<String>> prevArgs = new HashMap<Object, List<String>>();
/**
* Constructor
*
* @param anIncludePath
* the python include path
* @param aClassLoader
* the java classloader
* @param preEvals
* python statements to be executed before the script of
* aFilePath is run. This should include the statement
* "derivParamDir = x" where x is the location of the derived
* parameter scripts
* @throws JepException
*/
public MasterDerivScript(String anIncludePath, ClassLoader aClassLoader,
List<String> preEvals) throws JepException {
super(anIncludePath, aClassLoader, preEvals);
}
public Object executeFunction(String name, List<Object> args)
throws JepException {
this.prevArgs.clear();
executeFunctionInternal(name, args);
jep.eval("globalsRef=globals()");
for (List<String> pArgs : prevArgs.values()) {
for (String arg : pArgs) {
jep.eval(arg + " = None");
jep.eval("del globalsRef['" + arg + "']");
}
}
this.prevArgs.clear();
return getExecutionResult();
}
private void executeFunctionInternal(String name, List<Object> args)
throws JepException {
StringBuilder functionCall = new StringBuilder();
functionCall.append(RESULT).append(" = execute(");
for (int i = 0; i < args.size(); i++) {
if (i != 0) {
functionCall.append(", ");
}
Object arg = args.get(i);
String argName = "arg" + Integer.toHexString((arg.hashCode()));
evaluateArgument(argName, arg);
functionCall.append(argName);
}
functionCall.append(")");
if (name.contains(".")) {
int lastIdx = name.lastIndexOf(".");
String functionName = name.substring(lastIdx + 1);
String path = name.substring(0, lastIdx);
jep.eval("from " + path + " import " + functionName + " as execute");
} else {
jep.eval("from " + name + " import execute");
}
jep.eval(functionCall.toString());
}
@Override
protected void evaluateArgument(String argName, Object argValue)
throws JepException {
if (prevArgs.containsKey(argValue)) {
jep.eval(argName + " = " + prevArgs.get(argValue).get(0));
} else if (argValue instanceof List) {
@SuppressWarnings({ "rawtypes" })
List valList = (List) argValue;
Object val = valList.get(0);
if (val instanceof CubeLevel) {
// process as cube
jep.eval("import numpy");
jep.eval(argName + " = []");
boolean cubeCreated = false;
@SuppressWarnings("unchecked")
List<CubeLevel<Object, Object>> levelList = valList;
StringBuilder temp = new StringBuilder();
for (int i = 0; i < levelList.size(); i++) {
CubeLevel<Object, Object> cubeLevel = levelList.get(i);
Object press = cubeLevel.getPressure();
String pressKey = argName + "_tmpPress";
if (press != null) {
pressKey += Integer.toHexString((press.hashCode()));
}
evaluateArgument(pressKey, press);
Object param = cubeLevel.getParam();
String paramKey = argName + "_tmpParam";
if (param != null) {
paramKey += Integer.toHexString((param.hashCode()));
}
evaluateArgument(paramKey, param);
if (!cubeCreated) {
jep.eval("import numpy");
jep.eval("paramShape = " + paramKey + ".shape");
jep.eval(argName
+ ".append(numpy.ndarray(["
+ valList.size()
+ ", paramShape[0], paramShape[1]], 'float32'))");
// handle pressure shaping
temp.setLength(0);
temp.append("if type(");
temp.append(pressKey);
temp.append(") == float:\n\t");
temp.append(argName);
temp.append(".append(numpy.ndarray([1,");
temp.append(valList.size());
temp.append("], 'float32'))\nelse:\n\t");
temp.append(argName);
temp.append(".append(numpy.ndarray([");
temp.append(valList.size());
temp.append(", paramShape[0], paramShape[1]], 'float32'))");
jep.eval(temp.toString());
cubeCreated = true;
}
jep.eval(argName + "[0][" + i + "] = " + paramKey);
temp.setLength(0);
temp.append("if type(");
temp.append(pressKey);
temp.append(") == float:\n\t");
temp.append(argName);
temp.append("[1][0][");
temp.append(i);
temp.append("] = ");
temp.append(pressKey);
temp.append("\nelse:\n\t");
temp.append(argName);
temp.append("[1][");
temp.append(i);
temp.append("] = ");
temp.append(pressKey);
jep.eval(temp.toString());
}
} else {
// process as list
if (valList.size() == 1) {
// treat as an unwrapped array
evaluateArgument(argName, val);
} else {
// create a list of arrays
jep.eval(argName + " = []");
for (int argIdx = 0; argIdx < valList.size(); argIdx++) {
val = valList.get(argIdx);
String argKey = argName + "_" + argIdx;
if (val != null) {
argKey += "_"
+ Integer.toHexString((val.hashCode()));
}
// setNumpy won't work with indexed objects
evaluateArgument(argKey, val);
jep.eval(argName + ".append(" + argKey + ")");
}
}
}
} else if (argValue instanceof IDataRecord[]) {
IDataRecord[] valList = (IDataRecord[]) argValue;
if (valList.length == 1) {
// treat as an unwrapped array
evaluateArgument(argName, valList[0]);
} else {
// create a list of arrays
jep.eval(argName + " = []");
for (int argIdx = 0; argIdx < valList.length; argIdx++) {
IDataRecord val = valList[argIdx];
jep.eval(argName + ".append(None)");
// setNumpy won't work with indexed objects
evaluateArgument("__tmp", val);
jep.eval(argName + "[" + argIdx + "] = __tmp");
}
jep.eval(argName + " = tuple(" + argName + ")");
}
} else if (argValue instanceof IDataRecord) {
setDataRecordArg(argName, (IDataRecord) argValue);
} else if (argValue instanceof float[]) {
float[] val = (float[]) argValue;
jep.setNumpy(argName, val, val.length, 1);
} else if (argValue instanceof int[]) {
int[] val = (int[]) argValue;
jep.setNumpy(argName, val, val.length, 1);
} else if (argValue instanceof Float) {
jep.set(argName, (argValue));
} else if (argValue instanceof DerivedParameterRequest) {
DerivedParameterRequest request = (DerivedParameterRequest) argValue;
executeFunctionInternal(request.getMethod(),
Arrays.asList(request.getArgumentRecords()));
jep.eval(argName + "=" + RESULT);
} else {
super.evaluateArgument(argName, argValue);
}
List<String> pArgs = prevArgs.get(argValue);
if (pArgs == null) {
pArgs = new ArrayList<String>();
prevArgs.put(argValue, pArgs);
}
pArgs.add(argName);
}
/**
* Retrieves the result of the method execution
*/
protected List<?> getExecutionResult() throws JepException {
// Retrieve the results and return them
List<IDataRecord> result = null;
Boolean isTuple = (Boolean) jep.getValue("isinstance(" + RESULT
+ ",tuple)");
if (isTuple) {
// figure out how long the tuple is
int lenTuple = ((Integer) jep.getValue("len(" + RESULT + ")"))
.intValue();
// create result as a list of arrays
result = new ArrayList<IDataRecord>(lenTuple);
jep.eval("__tmp = " + RESULT);
// get each array and put it in the list
jep.eval(RESULT + " = __tmp[0]");
getExecutionResult(result, null);
long[] shape = result.get(0).getSizes();
for (int tupleElem = 1; tupleElem < lenTuple; tupleElem++) {
jep.eval(RESULT + " = __tmp[" + tupleElem + "]");
getExecutionResult(result, shape);
}
} else {
result = new ArrayList<IDataRecord>(1);
getExecutionResult(result, null);
}
jep.eval("del globals()['" + RESULT + "']");
return result;
}
protected void getExecutionResult(List<IDataRecord> result, long[] shape)
throws JepException {
filterResult();
// create result as a list with a single float array
Object valObj = jep.getValue(RESULT);
if (valObj instanceof float[]) {
if (shape == null) {
shape = getResultShape();
}
result.add(new FloatDataRecord(DATA_NAME, "", (float[]) valObj,
shape.length, shape));
} else if (valObj instanceof double[]) {
if (shape == null) {
shape = getResultShape();
}
double[] dData = (double[]) valObj;
float[] fData = new float[dData.length];
for (int i = 0; i < dData.length; i++) {
fData[i] = (float) dData[i];
}
result.add(new FloatDataRecord(DATA_NAME, "", fData, shape.length,
shape));
} else if (valObj instanceof int[]) {
if (shape == null) {
shape = getResultShape();
}
result.add(new IntegerDataRecord(DATA_NAME, "", (int[]) valObj,
shape.length, shape));
} else if (valObj instanceof byte[]) {
if (shape == null) {
shape = getResultShape();
}
result.add(new ByteDataRecord(DATA_NAME, "", (byte[]) valObj,
shape.length, shape));
} else if (valObj instanceof String[]) {
if (shape == null) {
shape = getResultShape();
}
result.add(new StringDataRecord(DATA_NAME, "", (String[]) valObj,
shape.length, shape));
} else {
// wrap value in containers to meet return type requirements
float[] oneVal = new float[1];
if (!(valObj instanceof Float)) {
// try to coerce it to a float
jep.eval(RESULT + " = float(" + RESULT + ")");
valObj = jep.getValue(RESULT);
}
oneVal[0] = ((Float) valObj).floatValue();
result.add(new FloatDataRecord(DATA_NAME, "", oneVal, 1,
new long[] { 1 }));
}
}
private void filterResult() throws JepException {
// String conversion
jep.eval("import numpy");
StringBuilder script = new StringBuilder();
script.append("if isinstance(" + RESULT + ", numpy.ndarray) and "
+ RESULT + ".dtype.kind == \"S\":\n");
script.append(" " + RESULT + "=" + RESULT + ".flatten().tolist()\n");
jep.eval(script.toString());
jep.eval("numpy = None");
jep.eval("del globals()['numpy']");
}
private long[] getResultShape() throws JepException {
jep.eval("import numpy");
int nDims = (Integer) jep.getValue("len(numpy.shape(" + RESULT + "))");
long[] dims = new long[nDims];
for (int i = 0; i < nDims; i++) {
dims[i] = ((Integer) jep.getValue("numpy.shape(" + RESULT + ")["
+ i + "]")).longValue();
}
if (dims.length == 2) {
// TODO: this is ridiculous
long swap = dims[0];
dims[0] = dims[1];
dims[1] = swap;
}
jep.eval("numpy = None");
jep.eval("del globals()['numpy']");
return dims;
}
private void setDataRecordArg(String argName, IDataRecord argValue)
throws JepException {
boolean reshape = true;
long[] sizes = argValue.getSizes();
if (argValue instanceof FloatDataRecord) {
FloatDataRecord record = (FloatDataRecord) argValue;
if (sizes.length == 2) {
jep.setNumpy(argName, record.getFloatData(), (int) sizes[0],
(int) sizes[1]);
reshape = false;
} else {
evaluateArgument(argName, record.getFloatData());
}
jep.eval("import numpy");
jep.eval(argName + "[" + argName + " <= -9999] = numpy.NaN");
jep.eval(argName + "[" + argName + " >= 999999] = numpy.NaN");
jep.eval("numpy = None");
jep.eval("del globals()['numpy']");
} else if (argValue instanceof IntegerDataRecord) {
IntegerDataRecord record = (IntegerDataRecord) argValue;
if (sizes.length == 2) {
jep.setNumpy(argName, record.getIntData(), (int) sizes[0],
(int) sizes[1]);
reshape = false;
} else {
evaluateArgument(argName, record.getIntData());
}
} else if (argValue instanceof ByteDataRecord) {
ByteDataRecord record = (ByteDataRecord) argValue;
if (sizes.length == 2) {
jep.setNumpy(argName, record.getByteData(), (int) sizes[0],
(int) sizes[1]);
reshape = false;
} else {
evaluateArgument(argName, record.getByteData());
}
} else {
jep.set(argName, argValue);
jep.eval("import numpy");
if (argValue instanceof LongDataRecord) {
jep.eval(argName + " = numpy.array(" + argName
+ ".getLongData())");
} else if (argValue instanceof ShortDataRecord) {
jep.eval(argName + " = numpy.array(" + argName
+ ".getShortData())");
} else if (argValue instanceof StringDataRecord) {
jep.eval(argName + " = numpy.array(" + argName
+ ".getStringData())");
} else {
jep.eval(argName + " = numpy.array(" + argName
+ ".getDataObject())");
}
jep.eval("numpy = None");
jep.eval("del globals()['numpy']");
}
if (reshape) {
jep.set("shape", argValue.getSizes());
jep.eval(argName + " = " + argName + ".reshape(tuple(shape))");
}
}
}

View file

@ -1,65 +0,0 @@
/**
* This software was developed and / or modified by Raytheon Company,
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
*
* U.S. EXPORT CONTROLLED TECHNICAL DATA
* This software product contains export-restricted data whose
* export/transfer/disclosure is restricted by U.S. law. Dissemination
* to non-U.S. persons whether in the United States or abroad requires
* an export license or other authorization.
*
* Contractor Name: Raytheon Company
* Contractor Address: 6825 Pine Street, Suite 340
* Mail Stop B8
* Omaha, NE 68106
* 402.291.0100
*
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
* further licensing information.
**/
package com.raytheon.uf.common.derivparam.python;
import java.util.List;
import jep.JepException;
import com.raytheon.uf.common.datastorage.records.IDataRecord;
import com.raytheon.uf.common.python.concurrent.IPythonExecutor;
/**
* Executor for calling executeFunction on a MasterDerivScript
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------ ---------- ----------- --------------------------
* Jun 04, 2013 2041 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
* @version 1.0
*/
public class MasterDerivScriptExecutor implements
IPythonExecutor<MasterDerivScript, List<IDataRecord>> {
private final String name;
private final List<Object> arguments;
public MasterDerivScriptExecutor(String name, List<Object> arguments) {
this.name = name;
this.arguments = arguments;
}
@SuppressWarnings("unchecked")
@Override
public List<IDataRecord> execute(MasterDerivScript script)
throws JepException {
return (List<IDataRecord>) script.executeFunction(name, arguments);
}
}

View file

@ -1,127 +0,0 @@
/**
* This software was developed and / or modified by Raytheon Company,
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
*
* U.S. EXPORT CONTROLLED TECHNICAL DATA
* This software product contains export-restricted data whose
* export/transfer/disclosure is restricted by U.S. law. Dissemination
* to non-U.S. persons whether in the United States or abroad requires
* an export license or other authorization.
*
* Contractor Name: Raytheon Company
* Contractor Address: 6825 Pine Street, Suite 340
* Mail Stop B8
* Omaha, NE 68106
* 402.291.0100
*
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
* further licensing information.
**/
package com.raytheon.uf.common.derivparam.python;
import java.io.File;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import jep.JepException;
import com.raytheon.uf.common.derivparam.library.DerivedParameterGenerator;
import com.raytheon.uf.common.localization.IPathManager;
import com.raytheon.uf.common.localization.LocalizationContext.LocalizationType;
import com.raytheon.uf.common.localization.LocalizationFile;
import com.raytheon.uf.common.localization.LocalizationUtil;
import com.raytheon.uf.common.localization.PathManagerFactory;
import com.raytheon.uf.common.python.PyUtil;
import com.raytheon.uf.common.python.concurrent.AbstractPythonScriptFactory;
/**
* Factory for creating and initializing MasterDerivScript.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------ ---------- ----------- --------------------------
* Jun 04, 2013 2041 bsteffen Initial creation
* Aug 26, 2013 2289 bsteffen Make number of deriv param threads
* configurable.
*
* </pre>
*
* @author bsteffen
* @version 1.0
*/
public class MasterDerivScriptFactory extends
AbstractPythonScriptFactory<MasterDerivScript> {
private static final int DEFAULT_MAX_THREADS = 3;
public static final String NAME = "DerivedParameterPython";
private static final String INTERFACE_SCRIPT = DerivedParameterGenerator.DERIV_PARAM_DIR
+ File.separator
+ "python"
+ File.separator
+ "DerivParamImporter.py";
public MasterDerivScriptFactory() {
super(NAME, getMaxThreadsProperty());
}
@Override
public MasterDerivScript createPythonScript() throws JepException {
IPathManager pm = PathManagerFactory.getPathManager();
File script = pm.getStaticFile(INTERFACE_SCRIPT);
// Get list of all files for search hierarch of CAVE_STATIC
LocalizationFile[] derivParamFiles = pm.listFiles(
pm.getLocalSearchHierarchy(LocalizationType.COMMON_STATIC),
DerivedParameterGenerator.DERIV_PARAM_DIR, null, false, false);
List<String> functionDirs = new ArrayList<String>(
derivParamFiles.length);
functionDirs.add(script.getParent());
Arrays.sort(derivParamFiles);
for (LocalizationFile file : derivParamFiles) {
if (file.isDirectory()
&& DerivedParameterGenerator.FUNCTIONS
.equals(LocalizationUtil.extractName(file.getName()))) {
// If it is a derived parameters functions directory, add to
// search list
functionDirs.add(file.getFile().getAbsolutePath());
}
}
// Create path from function dir list
String PATH = PyUtil.buildJepIncludePath(functionDirs
.toArray(new String[functionDirs.size()]));
List<String> preEvals = new ArrayList<String>(2);
preEvals.add("import DerivParamImporter");
StringBuilder cmd = new StringBuilder(200);
cmd.append("sys.meta_path.append(DerivParamImporter.DerivParamImporter(");
// Pass in directories to search based on function directories
int size = functionDirs.size() - 1;
for (int i = size; i > 0; --i) {
if (i < size) {
cmd.append(", ");
}
cmd.append("'").append(functionDirs.get(i)).append("'");
}
cmd.append("))");
preEvals.add(cmd.toString());
return new MasterDerivScript(PATH,
MasterDerivScript.class.getClassLoader(), preEvals);
}
private static int getMaxThreadsProperty() {
return Integer.getInteger(
"com.raytheon.uf.viz.derivparam.python.threads",
DEFAULT_MAX_THREADS);
}
}

View file

@ -1,94 +0,0 @@
/**
* This software was developed and / or modified by Raytheon Company,
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
*
* U.S. EXPORT CONTROLLED TECHNICAL DATA
* This software product contains export-restricted data whose
* export/transfer/disclosure is restricted by U.S. law. Dissemination
* to non-U.S. persons whether in the United States or abroad requires
* an export license or other authorization.
*
* Contractor Name: Raytheon Company
* Contractor Address: 6825 Pine Street, Suite 340
* Mail Stop B8
* Omaha, NE 68106
* 402.291.0100
*
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
* further licensing information.
**/
package com.raytheon.uf.common.derivparam.python.function;
import jep.INumpyable;
/**
* Calls {@link com.raytheon.uf.common.wxmath.CapeFunc} and transforms the
* output into an INumpyable.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------ ---------- ----------- --------------------------
* Aug 13, 2013 njensen Initial creation
*
* </pre>
*
* @author njensen
* @version 1.0
*/
public class CapeFuncPythonAdapter {
public static class CapeCinPair implements INumpyable {
private final int nx;
private final int ny;
private final float[] cape;
private final float[] cin;
public CapeCinPair(int nx, int ny, float[] cape, float[] cin) {
this.nx = nx;
this.ny = ny;
this.cape = cape;
this.cin = cin;
}
@Override
public Object[] getNumpy() {
return new Object[] { cape, cin };
}
@Override
public int getNumpyX() {
return nx;
}
@Override
public int getNumpyY() {
return ny;
}
}
public static INumpyable capeFunc(float usetv, float[] p_dat,
float[] tve_dat, float[] p0, float[] th0, float[] sh0, int nx,
int ny, int nz) {
float[][] result = com.raytheon.uf.common.wxmath.CapeFunc.capeFunc(
usetv, p_dat, tve_dat, p0, th0, sh0, nx, ny, nz);
return new CapeCinPair(ny, nx, result[0], result[1]);
}
public static INumpyable capeFuncTop(float usetv, float[] p_dat,
float[] tve_dat, float[] p0, float[] th0, float[] sh0,
float[] ptop, int nx, int ny, int nz) {
float[][] result = com.raytheon.uf.common.wxmath.CapeFunc.capeFuncTop(
usetv, p_dat, tve_dat, p0, th0, sh0, ptop, nx, ny, nz);
return new CapeCinPair(ny, nx, result[0], result[1]);
}
}

View file

@ -1,55 +0,0 @@
/**
* This software was developed and / or modified by Raytheon Company,
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
*
* U.S. EXPORT CONTROLLED TECHNICAL DATA
* This software product contains export-restricted data whose
* export/transfer/disclosure is restricted by U.S. law. Dissemination
* to non-U.S. persons whether in the United States or abroad requires
* an export license or other authorization.
*
* Contractor Name: Raytheon Company
* Contractor Address: 6825 Pine Street, Suite 340
* Mail Stop B8
* Omaha, NE 68106
* 402.291.0100
*
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
* further licensing information.
**/
package com.raytheon.uf.common.derivparam.python.function;
import jep.INumpyable;
import com.raytheon.uf.common.python.PythonNumpyFloatArray;
/**
* Calls {@link com.raytheon.uf.common.wxmath.DCapeFunc} and transforms the
* output into an INumpyable for python.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------ ---------- ----------- --------------------------
* Aug 13, 2013 njensen Initial creation
*
* </pre>
*
* @author njensen
* @version 1.0
*/
public class DCapeFuncPythonAdapter {
public static INumpyable dcapeFunc(float usetv, float[] p_dat,
float[] t_dat, float[] td_dat, float[] p0, float[] th0,
float[] sh0, int nx, int ny, int nz, float max_evap, float max_rh) {
float[] result = com.raytheon.uf.common.wxmath.DCapeFunc.dcapeFunc(
usetv, p_dat, t_dat, td_dat, p0, th0, sh0, nx, ny, nz,
max_evap, max_rh);
return new PythonNumpyFloatArray(result, ny, nx);
}
}

View file

@ -1,53 +0,0 @@
/**
* This software was developed and / or modified by Raytheon Company,
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
*
* U.S. EXPORT CONTROLLED TECHNICAL DATA
* This software product contains export-restricted data whose
* export/transfer/disclosure is restricted by U.S. law. Dissemination
* to non-U.S. persons whether in the United States or abroad requires
* an export license or other authorization.
*
* Contractor Name: Raytheon Company
* Contractor Address: 6825 Pine Street, Suite 340
* Mail Stop B8
* Omaha, NE 68106
* 402.291.0100
*
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
* further licensing information.
**/
package com.raytheon.uf.common.derivparam.python.function;
import jep.INumpyable;
import com.raytheon.uf.common.python.PythonNumpyFloatArray;
/**
* Calls {@link com.raytheon.uf.common.wxmath.DistFilter} and transforms the
* output into an INumpyable.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------ ---------- ----------- --------------------------
* Aug 13, 2013 njensen Initial creation
*
* </pre>
*
* @author njensen
* @version 1.0
*/
public class DistFilterPythonAdapter {
public static INumpyable filter(float[] input, float npts, int nx, int ny,
int times) {
float[] result = com.raytheon.uf.common.wxmath.DistFilter.filter(input,
npts, nx, ny, times);
return new PythonNumpyFloatArray(result, nx, ny);
}
}

View file

@ -1,47 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import add,array,zeros_like
import Vector
def execute(*args):
""" Perform scalar or vector addition
"""
if len(args) == 1 and isinstance(args[0], list):
return execute(*args[0])
elif isinstance(args[0], tuple):
return vectorAddition(args)
else:
return scalarAddition(args)
def scalarAddition(args):
return reduce(add, args)
def vectorAddition(args):
uResult = zeros_like(args[0][0])
vResult = zeros_like(args[0][0])
for u, v in args:
uResult += u
vResult += v
return Vector.componentsTo(uResult, vResult)

View file

@ -1,59 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file AdiabaticTemp.py
from numpy import exp
# Constants from the old C code.
const1 = 26.660820
const2 = 0.0091379024
const3 = 6106.396
const4 = 0.622
const5 = 2740.0
##
# The method that does the actual calculation.
# Adapted from the old adiabatic_te.c function.
# This method skips conversion of temp and pressure to NaN arrays
# and back again. It can be called directly if temp and pressure are
# guaranteed to be valid, or if they are already NaN arrays.
#
# @param temp: Temperature (K)
# @param pressure: Pressure (mb)
# @return: Adiabatic temperature (K)
# @rtype: numpy masked array
def calculate(temp, pressure):
ee = exp(const1 - const2*temp - const3/temp)
pme = pressure-ee
ee *= const4
ee = ee / pme
result = temp * exp(const5*ee/temp)
return result
##
# Calculate adaibatic temperature.
#
# @param temp: Temperature (K)
# @param pressure: Pressure (mb)
# @return: Adiabatic temperature
def execute(temp, pressure):
result = calculate(temp, pressure)
return result

View file

@ -1,86 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import empty, shape, NaN
##
# Calculate the advection of quantity by Vector.
#
# @param Vector: A 2-tuple (U,V) of vector components. U and V must have the same shape,
# be rank 2 or more, and have a shape of 3 or more in the first two dimensions.
# @param quantity: A concentration array. This must be an array with the same dimensions as U
# and V.
# @param dx: The spacing between adjacent data points in the X direction. This can be an array
# with the same dimensions as U or a scalar.
# @param dy: The spacing between adjacent data points in the Y direction. This can be an array
# with the same dimensions as U or a scalar.
# @return: The advection array.
# @rtype: An array of at least rank 2 with a shape of 3 or more in the first two dimensions.
# The outer edges of the returned cannot be calculated and are set to NaN.
def execute(U, V, quantity, dx, dy):
"""Calculate the advection of Vector.
Parameters:
Vector - tuple(U,V) of arrays, at least 3x3
quantity - quantity to advect: an array the same shape as U
dx - X dimension spacing: an array the same shape as U or a scalar
dy - Y dimension spacing: an array the same shape as U or a scalar
Returns:
Advection: An array the same shape as U
"""
result = empty(shape(U), dtype=U.dtype)
result[0,:] = NaN
result[-1,:] = NaN
result[1:-1,0] = NaN
result[1:-1,-1] = NaN
# If dx and dy are matrices, we never use the outer edge,
# so strip it off so we don't have to use slice notation in the math.
# If they're actually scalars or 1-element matrices, we can't
# slice them so don't try.
shapedx = shape(dx)
if len(shapedx) < sum(shapedx):
dx = dx[1:-1,1:-1]
shapedy = shape(dy)
if len(shapedy) < sum(shapedy):
dy = dy[1:-1, 1:-1]
# create a partial answer from horizontal neighbors, U, and dx
ans = quantity[1:-1,0:-2] - quantity[1:-1,2:]
# U data is coming in negated
ans *= -U[1:-1,1:-1]
ans /= dx
# create another partial answer from vertical neighbors, V, and dy
term = quantity[0:-2,1:-1] - quantity[2:,1:-1]
term *= V[1:-1,1:-1]
term /= dy
# average the two partial results
ans += term
ans /= 2
# put ans in the middle block of result
# Answer is reversed
result[1:-1,1:-1] = -ans
return result

View file

@ -1,33 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
# ----------------------------------------------------------------
#
#
# ----------------------------------------------------------------
import numpy
const1 = 33.86389
const2 = 100
def execute(altimeter, accum_altimeter24):
ALT0 = numpy.where(altimeter > 0 , altimeter, numpy.NaN)
ALT24 = numpy.where(accum_altimeter24 > 0 , accum_altimeter24, numpy.NaN)
return (ALT0 - ALT24) * const1 * const2

View file

@ -1,37 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import power
##
# Routine to pressure from elevation and altimeter setting.
#
# @param alt: Altimeter setting (X)
# @param z: Elevation in meters.
# @return: Pressure (X)
def execute(alt, z):
T0 = 288.0
gamma = 0.0065
g_Rgamma = 5.2532
result = (T0-gamma*z)/T0
result = power(result, g_Rgamma)
result = alt*result
return result

View file

@ -1,48 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from Add import execute as Add
from Divide import execute as Divide
def execute(*args):
""" Calculate the arithmetic average of any number of arguments.
"""
if len(args) == 1 and isinstance(args[0], list):
return execute(*args[0])
return Divide(Add(*args),float(len(args)))
def test():
from numpy import array
if not( all(execute(array([1.,2.]),array([3.,4.])) == array([2,3]))):
raise Exception
from Vector import execute as Vector
Vector1 = Vector(array([1.,2.]),array([0.,270.]),True)
Vector2 = Vector(array([3.,4.]),array([180.,90.]),True)
(mag, dir, u, v) = execute(Vector1,Vector2)
if not( all(mag == array([1., 1.])) and all(dir.round(0) == array([180., 90.]))):
raise Exception
print "Average Test Complete"

View file

@ -1,23 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
import numpy
def execute(POP_hour_bestCatStation):
return (POP_hour_bestCatStation <= 1).astype('float32')

View file

@ -1,98 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from com.raytheon.uf.viz.derivparam.python.function import CapeFuncPythonAdapter as CapeFunc
from numpy import zeros
def execute(*args):
return __execute(*args)[0]
def __execute(*args):
""" Calculate Convective Available Potential Energy.
@param temperature
if there are five or six input arguments the leading argument is
3D (virtual) temperature.
Otherwise, the leading inputs are a list of pressure values and a
list of corresponding (virtual) temperature values.
@param pressure
@param potentialTemperature
@param specificHumidity
of the starting parcel
@param useVirtualTemp
if true (1) use virtual temperatures, otherwise (0) use plain
temperatures
@param upperTerminationPressure
optional argument designating the upper pressure at which the CAPE
computation terminates
@return a scalar
Temperature must be in Kelvins
"""
argLength = len(args)
upperTerminationPressure = None
useVirtualTemp = args[argLength - 1]
specificHumidity = args[argLength - 2]
potentialTemperature = args[argLength - 3]
pressure = args[argLength - 4]
temperatureListLengthOffset = 4
if type(useVirtualTemp) != float:
# an upperTerminationPressure has been supplied
upperTerminationPressure = useVirtualTemp
useVirtualTemp = args[argLength - 2]
specificHumidity = args[argLength - 3]
potentialTemperature = args[argLength - 4]
pressure = args[argLength - 5]
temperatureListLengthOffset = 5
pressureValues = []
temperatureValues = []
if type(args[0]) == list:
# we have 3D temperature data
temperatureValues = args[0][0]
pressureValues = zeros(temperatureValues.shape, temperatureValues.dtype)
else:
# we have a list of pressures followed by temperatures
temperatureListLength = (argLength - temperatureListLengthOffset) / 2
for i in range(0, temperatureListLength):
pressureValues.append(args[i])
temperatureValues.append(args[i + temperatureListLength])
# expand each pressure value to a full grid size
for pressureLevel in range(0, pressureValues.shape[0]):
gridLevelValues = pressureValues[pressureLevel]
gridLevelValues[:] = args[0][1][0][pressureLevel]
pressureValue = pressure
pressure = zeros(temperatureValues.shape, temperatureValues.dtype)
pressure[:] = pressureValue
if upperTerminationPressure is None:
threeDshape = pressureValues.shape
return CapeFunc.capeFunc(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, threeDshape[1], threeDshape[2], threeDshape[0]).__numpy__
else:
threeDshape = pressureValues.shape
return CapeFunc.capeFuncTop(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, upperTerminationPressure, threeDshape[1], threeDshape[2], threeDshape[0]).__numpy__

View file

@ -1,24 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from Cape import __execute
def execute(*args):
return __execute(*args)[1]

View file

@ -1,90 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import zeros, float32, NaN, isnan
##
# Sample a column of profiler data in a way which matches A1
# The most notable feature is that A1 stores all profiler data at evenly spaced level
# where as A2 stores only the levels which have data. This function essentially
# reinserts the missing level and then selects the closest level
#
def execute(paramArray, vertArray, numLevels, vertPoints, sfcParams=None, sfcVerts=None):
ret = zeros(numLevels.shape, 'float32')
ret.fill(NaN)
for i in range(len(vertArray)):
count = numLevels[i]
verts = vertArray[i][:count]
params = paramArray[i][:count]
if isinstance(vertPoints, float):
vertPoint = vertPoints
elif isinstance(vertPoints, float32):
vertPoint = vertPoints
else:
vertPoint = vertPoints[i]
if sfcParams == None:
sfcParam = NaN
elif isinstance(sfcParams, float):
sfcParam = sfcParams
elif isinstance(sfcParams, float32):
sfcParam = sfcParams
else:
sfcParam = sfcParams[i]
if sfcVerts == None:
sfcVert = NaN
elif isinstance(sfcVerts, float):
sfcVert = sfcVerts
elif isinstance(sfcVerts, float32):
sfcVert = sfcVerts
else:
sfcVert = sfcVerts[i]
# expand the input vertices to include surface and every 250 meters
inVert = [sfcVert]
inData = [sfcParam]
for j in range(0,len(verts)):
# Profiler data should be spaced every 250 from 500 to 9250, and then every 1000 above that
heightAGL = verts[j] - sfcVert
if heightAGL < 0:
continue
if heightAGL > 9250:
index = 36 + int((verts[j] - sfcVert - 9250)/1000)
else:
index = int((verts[j] - sfcVert - 250)/250)
while len(inVert) <= index:
inData.append(NaN)
if len(inVert) > 36:
inVert.append(9250 + sfcVert + 1000*(len(inVert) - 36))
else:
inVert.append(250 + sfcVert + 250*len(inVert))
inVert[index] = verts[j]
inData[index] = params[j]
# add a NaN on the top to prevent interpolating over.
inData.append(NaN)
if len(inData) > 36:
inVert.append(9250 + sfcVert + 1000*(len(inVert) - 36))
else:
inVert.append(250 + sfcVert + 250*len(inVert))
#grab the closest level
bestDist = 999999
for j in range(0,len(inVert)):
dist = abs(inVert[j]-vertPoint)
if dist<bestDist:
bestDist = dist
ret[i] = inData[j];
return ret

View file

@ -1,104 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file CompBy.py
from numpy import cos, sin, sqrt
import Vector
##
# This function is based on the comp_by.f Fortran function by J. Ramer, Jul 1 2003.
# It returns the component of Vecs in the direction of ThGrd, normalized by the
# magnitude of ThGrd, as a two-component vector. ThGrd is rotated by Angle before
# the component calculation.
#
# comp_by.f performs different calculations based on whether Angle is supplied at all,
# whether it is an integer, and the value of abs(Angle)/1000. This function corresponds
# to the calculations performed in comp_by.f when Angle is an integer and abs(Angle)/1000
# is equal to 2. Unless Angle is 0. ( see comp_by.f )
#
# @param Vecs: Data vector as a tuple of arrays
# @param ThGrd: Potential temperature gradient vector as a tuple of arrays
# @param Angle: Degrees through which to rotate first vector before
# dotting it with the second vector. May be a scalar or array.
# Any 1000s component is ignored for compatibility with old code;
# otherwise, valid values fall in the range [-180..180]. Angles
# outside the allowed range are treated as an angle of 0.
#
# @return: a tuple of 2 arrays representing the vector components of Vecs in the
# direction of ThGrd, normalized by the magnitude of ThGrd.
#
def execute(Vecs, ThGrd, Angle=0):
""
# pi/180 for degrees to radians conversion
dgtord = 0.01745329252
Vecs_U, Vecs_V = Vecs
ThGrd_U, ThGrd_V = ThGrd
if Angle != 0:
# Strip off any 1000s from angle.
# The second part makes sure -2090 ends up as -90 rather than 910.
# This works with scalars or arrays, where "if angle<0:" wouldn't.
Angle = Angle % 1000 + (Angle < 0) * -1000
# Angle has to be in the [-180..180] range.
# If not, convert it to zero.
Angle = Angle - (Angle > 180) * Angle
Angle = Angle - (Angle < -180) * Angle
# convert angle in degrees to U and V vector components
rotate_u = cos( dgtord * Angle )
rotate_v = sin( dgtord * Angle )
# find the rotated vector
rotated_u = ThGrd_U * rotate_u + ThGrd_V * rotate_v
rotated_v = ThGrd_V * rotate_u - ThGrd_U * rotate_v
mag = rotated_u * rotated_u + rotated_v * rotated_v
# we have to divide by the magnitude, so mask away any zero values
# find the components
zeroMag = (mag == 0)
mag[zeroMag] = 1
mag = (rotated_u * Vecs_U + rotated_v * Vecs_V)/mag
mag[zeroMag] = 0
comp = rotated_u * mag
comp2 = rotated_v * mag
return Vector.componentsTo(comp, comp2)
else:
# initialize to u and v from ThGrd
rotated_u = ThGrd_U
rotated_v = ThGrd_V
#calculate the magnitude
mag = rotated_u * rotated_u + rotated_v * rotated_v
mag = sqrt(mag)
#magical math that I don't know what it does
zeroMag = (mag == 0)
mag[zeroMag] = 1
mag = (rotated_u * Vecs_U + rotated_v * Vecs_V)/mag
mag[zeroMag] = 0
return mag

View file

@ -1,54 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import zeros
from numpy import log
from numpy import sqrt
CONST_A = 0.0091379024
CONST_B = 6106.396
CONST_C = 223.1986
CONST_D = 0.0182758048
CONST_E = 0.37329638
CONST_F = 41.178204
CONST_G = 0.0015945203
CONST_H = 3.498257
##
# Calculate condensation pressure from pressure, temperature, and relative
# humidity
# @attention: Result may contain NaN
#
# @param p: Pressure in millibars
# @type p: scalar or numpy array
# @param t: Temperature in degrees Kelvin
# @type t: scalar or numpy array
# @param rh: Relative humidity from 0.0 to 100.0
# @type rh: scalar or numpy array
#
# @return: Condensation pressure in millibars
# @rtype: numpy array of float
def execute(p, t, rh):
rhqc = rh.clip(1.0,100.0)
b = CONST_A * t + CONST_B / t - log(rhqc/100.0)
tdp = (b - sqrt(b**2.0 - CONST_C)) / CONST_D
tcp = tdp - (t-tdp) * (-CONST_E + CONST_F / t + CONST_G * tdp)
q = p * (tcp/t) ** CONST_H
return q

View file

@ -1,66 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# SOFTWARE HISTORY
#
# Date Ticket# Engineer Description
# ------------ ---------- ----------- --------------------------
# Feb 17, 2010 4502 jelkins Initial Creation.
# Jun 05, 2013 2043 bsteffen Ported from meteolib C
from numpy import zeros
from com.raytheon.uf.viz.derivparam.python.function import DCapeFuncPythonAdapter as DCapeFunc
def execute(threeDtemperature, threeDdewpoint, pressure, potentialTemperature, specificHumidity,maxEvaporation,maxRelativeHumidity,useVirtualTemp):
""" Calculate Downdraft Convective Available Potential Energy
@param threeDtemperature
3D array of temperatures
@param threeDdewpoint
3D array of dewpoints
@param pressure
@param potentialTemperature
@param specificHumidity
@param maxEvaporation
maximum amount of liquid water available to evaporate into the parcel
as it descends, in g/g
@param maxRelativeHumidity
the desired maximum RH of the descending parcel as it reaches the surface
@param useVirtualTemp
use virtual (1) or plain (0) temperatures
@return
"""
threeDpressure = zeros(threeDtemperature[0].shape, threeDtemperature[0].dtype)
# fill in the 3D pressure
for pressureLevel in range(0, threeDpressure.shape[0]):
pressureGrid = threeDpressure[pressureLevel]
pressureGrid[:] = threeDtemperature[1][0][pressureLevel]
# expand the single pressure value into a complete grid
pressureValue = pressure
pressure = zeros(potentialTemperature.shape, potentialTemperature.dtype)
pressure[:] = pressureValue
threeDshape = threeDpressure.shape
return DCapeFunc.dcapeFunc(useVirtualTemp, threeDpressure, threeDtemperature[0], threeDdewpoint[0], pressure, potentialTemperature, specificHumidity, threeDshape[1], threeDshape[2], threeDshape[0], maxEvaporation, maxRelativeHumidity).__numpy__[0]

View file

@ -1,74 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import shape, empty, hypot, NaN
##
# Calculate the deformation of Vector.
#
# This function is currently only used by Tdef. It was put in the
# functions module to keep it in the same module as other functions
# derived from g2gkinematics.f.
#
# @param U: U and V must be the same size,
# at least 2d, and at least 3x3 in the first two dimensions.
# @param V: U and V must be the same size,
# at least 2d, and at least 3x3 in the first two dimensions.
# @param dx: The spacing between data points in the X direction (array or scalar)
# @param dy: The spacing between data points in the Y direction (array or scalar)
# @return: The deformation array of Vector.
def execute(U, V, dx, dy):
# create an empty array the same dimensions as U
result = empty(U.shape, U.dtype)
# flag edges as invalid
result[0,:] = NaN
result[-1,:] = NaN
result[1:-1,0] = NaN
result[1:-1,-1] = NaN
# strip off edges of dx if it's not a scalar
shapedx = shape(dx)
if len(shapedx) < sum(shapedx):
dx = dx[1:-1,1:-1]
# strip off edges of dy if it's not a scalar
shapedy = shape(dy)
if len(shapedy) < sum(shapedy):
dy = dy[1:-1,1:-1]
# Calculate deformation vector components.
qqq = 0.5/dx
www = 0.5/dy
dst = (U[1:-1,0:-2] - U[1:-1,2:]) * qqq
dst += (V[0:-2,1:-1] - V[2:,1:-1]) * www
dsh = (U[0:-2,1:-1] - U[2:,1:-1]) * qqq
dsh += (V[1:-1,2:] - V[1:-1,0:-2]) * www
# find vector magnitude of dst and dsh
ans = hypot(dst, dsh) * 100000.0
# put ans in the valid block of result
result[1:-1,1:-1] = ans
return result

View file

@ -1,83 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import sqrt
from numpy import shape
from numpy import empty
from numpy import hypot
from numpy import where
from numpy import NaN
import Vector as Vector
##
# Calculate the deformation components of Vector.
#
# @param Vector: A tuple(U,V) of numpy arrays. U and V must be the same size,
# at least 2d, and at least 3x3 in the first two dimensions.
# @param dx: The spacing between data points in the X direction (array or scalar)
# @param dy: The spacing between data points in the Y direction (array or scalar)
# @return: The deformation array of Vector.
def execute(U, V, dx, dy):
# create an empty array the same dimensions as U
resultx = empty(U.shape, U.dtype)
resulty = empty(U.shape, U.dtype)
# flag edges as invalid
resultx[0,:] = NaN
resultx[-1,:] = NaN
resultx[1:-1,0] = NaN
resultx[1:-1,-1] = NaN
resulty[0,:] = NaN
resulty[-1,:] = NaN
resulty[1:-1,0] = NaN
resulty[1:-1,-1] = NaN
# strip off edges of dx if it's not a scalar
shapedx = shape(dx)
if len(shapedx) < sum(shapedx):
dx = dx[1:-1,1:-1]
# strip off edges of dy if it's not a scalar
shapedy = shape(dy)
if len(shapedy) < sum(shapedy):
dy = dy[1:-1,1:-1]
# Calculate deformation vector components.
qqq = 0.5/dx
www = 0.5/dy
dst = (U[1:-1,0:-2] - U[1:-1,2:]) * qqq
dst += (V[0:-2,1:-1] - V[2:,1:-1]) * www
dsh = (U[0:-2,1:-1] - U[2:,1:-1]) * qqq
dsh += (V[1:-1,2:] - V[1:-1,0:-2]) * www
qqq = dst*dst + dsh*dsh
www = sqrt(qqq)
ans = sqrt((qqq+www*dst)/2)
Q = where(ans != 0.0, www*dsh/(2*ans), www)
# put ans in the valid block of result
resultx[1:-1,1:-1] = Q
resulty[1:-1,1:-1] = ans
return Vector.execute(resultx,resulty)

View file

@ -1,57 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file Derivative.py
from numpy.ma.core import masked_greater
from numpy.ma.core import masked_values
from numpy.ma.core import filled
import Vector
##
# Calculate the derivative of A with respect to B.
#
# @param A0: Initial value of A
# @type A0: numpy array
# @param A1: Final value of A
# @type A0: numpy array
# @param B0: Initial value of B
# @type A0: numpy array
# @param B1: Final value of B
# @type A0: numpy array
# @return: The derivative of A with respect to B
# @rtype: numpy array
def execute(A0, A1, B0, B1):
"Calculate the derivative of A with respect to B."
if isinstance(A0, tuple):
# Do derivitive of components of the vector
u0,v0 = A0
u1,v1 = A1
uDeriv = execute(u0, u1, B0, B1)
vDeriv = execute(v0, v1, B0, B1)
return Vector.componentsTo(uDeriv, vDeriv)
Adiff = A1-A0
Bdiff = B1-B0
result = Adiff / Bdiff
return result

View file

@ -1,91 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file DgeoComps.py
from numpy import isscalar
from numpy import copy as Copy
from numpy import NaN
from numpy import ones
from numpy import shape
##
# Compute the components of geostrophic wind for this level.
# @attention: result may contain NaN.
#
# @param height: Height (m)
# @type height: 2D numpy array, at least 3x3
# @param dx: Spacing between data points in the X direction
# @type dx: scalar or 2D numpy array at least 3x3
# @param dy: Spacing between data points in the Y direction
# @type dy: scalar or 2D numpy array at least 3x3
# @param coriolis: Coriolis force (kg-m/s^2)
# @type coriolis: scalar or 2D numpy array at least 3x3
# @return: d/dx of u component of geostrophic wind
# d/dy of u component of geostrophic wind
# d/dx of v component of geostrophic wind
# d/dy of v component of geostrophic wind
# @rtype: tuple(dugdx, dugdy, dvgdx, dvgdy)
def execute(height, dx, dy, coriolis):
"Compute the components of geostrophic wind for this level."
# Generate output arrays as arrays of NaN,
# the same size and dtype as Height.
# dugdx = masked_where(True, Height) -- not needed
dugdy = ones(shape(height), dtype=height.dtype) * NaN
dvgdx = Copy(dugdy)
dvgdy = Copy(dugdy)
gravAcc = 9.806
# use cropped versions of dx, dy, and coriolis in calculations
if not isscalar(dx):
dx = dx[1:-1,1:-1]
if not isscalar(dy):
dy = dy[1:-1,1:-1]
if not isscalar(coriolis):
coriolis = coriolis[1:-1,1:-1]
# intermediate values that are used repeatedly
qqq = gravAcc / coriolis
www = height[1:-1,1:-1] * 2.0
# make calculations in arrays with outer edge cropped
# then paste into full-sized array
croppedDvgdy = height[2:,2:] + height[0:-2,0:-2]
croppedDvgdy -= height[2:,0:-2] + height[0:-2,2:]
croppedDvgdy *= qqq
croppedDvgdy /= 4 * dx * dy
dvgdy[1:-1,1:-1] = croppedDvgdy
dugdx = -dvgdy#*0.0
croppedDugdy = www - height[2:,1:-1] - height[0:-2,1:-1]
croppedDugdy *= qqq
croppedDugdy /= dy * dy
dugdy[1:-1,1:-1] = croppedDugdy
croppedDvgdx = height[1:-1,2:] + height[1:-1,0:-2] - www
croppedDvgdx *= qqq
croppedDvgdx /= dx * dx
dvgdx[1:-1,1:-1] = croppedDvgdx
return (dugdx, dugdy, dvgdx, dvgdy)

View file

@ -1,38 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from Add import execute as Add
from numpy import multiply, copy
def execute(*args):
""" Perform scalar or vector subtraction.
"""
diffArgs = list(args)
if type(diffArgs[0]) == tuple:
for i in range(1, len(diffArgs)):
diffArgs[i] = (-diffArgs[i][0], -diffArgs[i][1])
return apply(Add, diffArgs)
else:
result = 0
result += diffArgs[0]
for i in range(1, len(diffArgs)):
result -= diffArgs[i]
return result

View file

@ -1,44 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
# ----------------------------------------------------------------
#
#
# ----------------------------------------------------------------
import numpy
import U
import V
from numpy.ma.core import filled
from numpy import arctan2
def execute(windSpeed, windDir, accum_windSpeed24, accum_windDir24):
U0 = U.calculate(windSpeed, windDir)
V0 = V.calculate(windSpeed, windDir)
U24 = U.calculate(accum_windSpeed24, accum_windDir24)
V24 = V.calculate(accum_windSpeed24, accum_windDir24)
DU = U0 - U24
DV = V0 - V24
WD = 57.2957 * arctan2(-DU,-DV)
WD = numpy.where(WD < 0, WD + 360, WD)
return WD

View file

@ -1,35 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# Calculate Direction
# ----------------------------------------------------------------
from numpy import arctan2
from numpy import any
def execute(uW, vW, copy=False):
# 180.0 / pi = 57.2957
dir = 57.2957 * arctan2(-uW,-vW)
dirN = dir < 0
if any(dirN):
dir[dirN] += 360.0
return dir

View file

@ -1,102 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import empty, isscalar, shape, NaN
##
# Calculate the divergence of Vector.
# Because of the way in which this is calculated, the outermost
# row and column of the output array is always invalid.
#
# @change: quan parm added 2008/06/23 wldougher
#
# @param Vector: 2-tuple (U,V) of arrays of values
# @param dx: 2D array of spacing between grid points
# @param dy: 2D array of spacing between grid points
# @param quan: Quantity parameter (optional; defaults to constant 1)
# @return: divergence
def execute(vec, dx, dy, quan=1.0):
"""Calculate the divergence of Vector from dx and dy.\n\
for each output cell, \n
u = Vector[0] \n
v = Vector[1] \n
diff_u = u[i+1,j]*q[i+1,j] - u[i-1,j]*q[i-1,j] \n
diff_v = u[i,j+1]*q[i,j+1] - u[i,j-1]*q[i,j-1] \n
dudx = diff_u/dx[i,j] \n
dvdy = diff_v/dy[i,j] \n
diverg[i,j] = (dudx + dvdy)/2 \n
"""
vec_U,vec_V = vec
vshape = shape(vec_U)
if len(vshape) < 2:
raise TypeError, "Divergence: Vector must be at least 2-D."
if vshape[0] < 3:
raise ValueError, "Divergence: Vector's first dimension must be at least 3."
if vshape[1] < 3:
raise ValueError, "Divergence: Vector's second dimension must be at least 3."
rslt = empty(vshape, dtype=vec_U.dtype)
rslt[0,:] = NaN
rslt[-1,:] = NaN
rslt[1:-1, 0] = NaN
rslt[1:-1, -1] = NaN
if shape(quan) == ():
quan_rl = quan
quan_tb = quan
else:
quan_rl = quan[1:-1]
quan_tb = quan[:,1:-1]
QU = quan_rl * vec_U[1:-1,:]
QV = quan_tb * vec_V[:,1:-1]
diff_U = QU[:,2:] - QU[:,0:-2]
diff_V = QV[0:-2,:] - QV[2:,:]
# if dx is a constant or one-element array, just divide by it.
# otherwise, divide each cell in diff_U by the corresponding cell in dx.
dxshape = shape(dx)
ldxs = len(dxshape)
if ldxs == 0 or sum(dxshape)==ldxs:
dudx = diff_U/dx
elif dxshape==vshape:
dudx = diff_U/dx[1:-1,1:-1]
else:
raise TypeError, "Divergence: dx must be a scalar or the same shape as Vector."
# if dy is a constant or one-element array, just divide by it.
# otherwise, divide each cell in diff_V by the corresponding cell in dy.
dyshape = shape(dy)
ldys = len(dyshape)
if ldys == 0 or sum(dyshape)==ldys:
dvdy = diff_V/dy
elif dyshape==vshape:
dvdy = diff_V/dy[1:-1,1:-1]
else:
raise TypeError, "Divergence: dy must be a scalar or the same shape as Vector."
rslt[1:-1,1:-1] = (dudx + dvdy)/2
return rslt

View file

@ -1,55 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from Multiply import execute as Multiply
from numpy import NaN, where
#
def execute(*args):
""" Divide a scalar by a scalar or a vector by a scalar.
"""
divArgs = list(args)
for i in range(1,len(divArgs)):
divArgs[i] = where(divArgs[i] == 0, NaN, 1/divArgs[i] )
return apply(Multiply,divArgs)
def test():
from numpy import array
if not( all(execute(array([2.,4.]),array([2.,2.])) == array([1.,2.]))):
raise Exception
from Vector import execute as Vector
vector = Vector(array([2.,4.]),array([0.,270.]),True)
(mag, dir, u, v) = execute(vector,2.)
if not( all(mag == array([1.,2.])) and all(dir == array([0.,270.]))):
raise Exception
if not( all(execute(array([2.,4.]),array([2.,0.])) == array([1.,NaN]))):
raise Exception
print "Divide Test Complete"

View file

@ -1,32 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
# ----------------------------------------------------------------
#
#
# ----------------------------------------------------------------
import numpy
def execute(dewpoint, dpFromTenths, accum_dewpoint24, accum_dpFromTenths24):
Td24 = numpy.where(numpy.isnan(dpFromTenths), dewpoint, dpFromTenths)
Td0 = (dewpoint * 1.8) + 32
Td24 = numpy.where(numpy.isnan(accum_dpFromTenths24), accum_dewpoint24, accum_dpFromTenths24)
Td24 = (accum_dpFromTenths24 * 1.8) + 32
return Td0 - Td24

View file

@ -1,85 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# Calculate Dewpoint Temperature (K) from Temperature (K) and Relative
# Humdity(0 to 100).
# ----------------------------------------------------------------
import numpy
from numpy import clip
from numpy import log
from numpy import sqrt
from numpy import zeros
from numpy import where
from numpy import equal
from numpy import isnan, NaN
from unit import celciusToKelvin
# build an array that contains the higher resolution dpFromTenths where possible, otherwise use the lower resolution dewpoint
def combineTemps(dewpoint,dpFromTenths):
return numpy.where(numpy.isnan(dpFromTenths), dewpoint, dpFromTenths)
##
# Calculate Dewpoint Temperature (K) from Temperature (K) and Relative
# Humdity(0 to 100).
# This function can operate on numpy arrays of the appropriate values.
#
# @param T: Temperature in degrees K
# @param RH: Relative humidity from 0 to 100
# @return: Dewpoint temperature in degrees K
# @rtype: numpy array of Python floats or Python float
def execute1(T,RH):
"Calculate dewpoint temperature(K) from temperature(K) and relative \
humidity(0 to 100)."
rhqc=clip(RH,1.0,100.0)
b=0.0091379024*T
b += 6106.396/T
b -= log(rhqc/100)
val = b*b
val -= 223.1986
val = sqrt(val)
DpT = b-val
DpT /= 0.0182758048
return DpT
def execute3(P,T,SHx):
eee=P*SHx/(622.0+0.378*SHx)
b=26.66082-log(eee)
result = (b-sqrt(b*b-223.1986))/0.0182758048
result[isnan(T)] = NaN
result[eee>980.5386] = NaN
result[eee<3.777647E-05] = NaN
return result
# @param dewpoint: dewpoint in degrees C
# @param dpFromTenths: dewpoint from thenths in degrees C
# @return: Dewpoint temperature in degrees K
# @rtype: numpy array of Python floats or Python float
def execute4(dewpoint,dpFromTenths):
return celciusToKelvin(combineTemps(dewpoint, dpFromTenths))
# @param dewpoint: dewpoint in degrees K
# @param temperature: Temperature in degrees K
# @param relHumidity: Relative humidity from 0 to 100
# @return: Dewpoint temperature in degrees K
# @rtype: numpy array of Python floats or Python float
def execute5(dewpoint,temperature,relHumidity):
TdRH = execute1(temperature, relHumidity)
return combineTemps(TdRH,dewpoint)

View file

@ -1,106 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import ndarray, float32, NaN
from numpy import sin, isnan, invert
from com.raytheon.uf.viz.derivparam.python.function import DistFilterPythonAdapter as DistFilter
MAX_WAVE_NUMBER = 15
##
# Apply a spatial filter to the first argument. The second argument is a constant which determines the distance over which to apply the filter, in km if positive, in number of grid points if negative.
#
# @param grid: the grid to filter
# @param distance: the distance over which to filter
# @param dx: Grid spacing in the x-direction (m)
# @param dy: Grid spacing in the y-direction (m)
# @param times: Number of times to filter
# @return: A filtered grid
def execute(input, dist, dx, dy, times=1):
if (dist<0):
npts = -dist;
else:
npts = dist * (1000/dx[input.shape[0]/2,input.shape[1]/2]);
return executeJava(input, npts, times)
#return executePython(input, npts, times)
def executeJava(input, npts, times):
output = DistFilter.filter(input, npts, input.shape[1], input.shape[0], times).__numpy__[0]
output[output==1e37] = NaN
return output
def executePython(input, npts, times=1):
n = int(npts + 0.99)
if (n<2):
n = 2;
elif (n>MAX_WAVE_NUMBER):
n = MAX_WAVE_NUMBER
d = (n+1)/2;
dd = d+d;
m = dd+1;
#Calculate wave table
waveno = 3.14159265/npts
jweights = ndarray([m,m], float32)
iweights = ndarray([m,m], float32)
for c in range(-d, d+1):
if (c != 0):
val = sin(waveno*c)/(waveno*c)
iweights[c+d,:] = val
jweights[:,c+d] = val
else:
iweights[c+d,:] = 1
jweights[:,c+d] = 1
weights = iweights*jweights
weights = weights/weights.sum()
output = input;
for i in range(times):
input = output;
output = ndarray(input.shape, float32)
for i in range(input.shape[0]):
for j in range(input.shape[1]):
wi1 = 0
wi2 = m
wj1 = 0
wj2 = m
i1 = i - d
i2 = i + d + 1
j1 = j - d
j2 = j + d + 1
if (i1 < 0):
wi1 = -i1
i1 = 0
if (i2 > input.shape[0]):
wi2 = - (i2 - input.shape[0])
i2 = input.shape[0]
if (j1 < 0):
wj1 = -j1
j1 = 0
if (j2 > input.shape[1]):
wj2 = - (j2 - input.shape[1])
j2 = input.shape[1]
notnanmask = invert(isnan(input[i1:i2,j1:j2]))
tot = weights[wi1:wi2,wj1:wj2][notnanmask].sum()
if tot >= 0.95:
output[i,j] = (input[i1:i2,j1:j2][notnanmask]*weights[wi1:wi2,wj1:wj2][notnanmask]).sum()/tot
else:
output[i,j] = NaN
return output

View file

@ -1,75 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
def execute(*args):
"""Test all derived functions for proper operation
"""
import Vector
reload(Vector)
Vector.test()
import Add
reload(Add)
Add.test()
import Multiply
reload(Multiply)
Multiply.test()
import Divide
reload(Divide)
Divide.test()
import Difference
reload(Difference)
Difference.test()
import Poly
reload(Poly)
Poly.test()
import Average
reload(Average)
Average.test()
import LinTrans
reload(LinTrans)
LinTrans.test()
import Test
reload(Test)
Test.test()
import Rotate
reload(Rotate)
Rotate.test()
import Magnitude
reload(Magnitude)
Magnitude.test()
import Cape
reload(Cape)
Cape.test()
print "-"*60
print "Function Test Complete"

View file

@ -1,43 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
import meteolib
from numpy import concatenate
from numpy import where
from numpy import equal
from numpy import log
from numpy import greater
from numpy import save
import P
def execute1(levels, staElev):
staElev = staElev.reshape(-1, 1);
hft2m = 30.48
return concatenate((staElev, levels*30.48), 1)
def execute2(pres):
return meteolib.ptozsa(pres)*4
def execute3(prCloudStation,lowCldStation,midCldStation,hiCldStation):
prCloudClg = P.execute6(prCloudStation,lowCldStation,midCldStation,hiCldStation)
return meteolib.ptozsa(prCloudClg)

View file

@ -1,82 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file Geowind.py
from numpy import isscalar,NaN
import Vector
g = 9.806 # Gravitational acceleration (m/s^2)
##
# Calculate geo wind from height, dx, dy, and coriolis parameter.
#
# @param Height: Height (m)
# @type Height: numpy array
# @param dx: Spacing between data points in X direction (m).
# @type dx: numpy array or scalar
# @param dy: Spacing between data points in Y direction (m).
# @type dy: numpy array or scalar
# @param coriolis: coriolis effect
# @type coriolis: numpy array or scalar
# @return: geological wind
# @rtype: tuple(mag,dir,U,V) of numpy arrays of float
def execute(Height, dx, dy, coriolis):
""
# assume dx, dy, and coriolis are OK
# Because we're using the adjacent points to calculate the result, we can't
# find values for points along the edges. Since we never use any points of
# dx, dy, and coriolis except the middle block, redefine them as the slice
# we actually use. This also allows us to deal with scalars.
if not isscalar(dx):
dx = dx[1:-1, 1:-1]
if not isscalar(dy):
dy = dy[1:-1, 1:-1]
if not isscalar(coriolis):
coriolis = coriolis[1:-1, 1:-1]
# Create the cropped answer arrays.
ans_V = Height[1:-1,2:] - Height[1:-1, 0:-2]
ans_V *= g
ans_V /= 2 * dy * coriolis
ans_U = Height[2:, 1:-1] - Height[0:-2, 1:-1]
ans_U *= g
ans_U /= 2 * dx * coriolis
# Create full-sized result arrays with all cells masked.
result_U = Height + NaN
result_V = Height + NaN
# Paste the cropped arrays into the valid portion of the answer arrays.
result_U[1:-1, 1:-1] = ans_U
result_V[1:-1, 1:-1] = ans_V
# Any masked cells become NaN
# This includes the outer edges and any cells that used a masked
# value from Height.
result_U = result_U
result_V = result_V
result = Vector.componentsTo(result_U, result_V)
return result

View file

@ -1,70 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file Gradient.py
from numpy import empty, shape, NaN
import Vector
##
# Calculate the X and Y gradient arrays of a 2+D array of scalars.
#
# @param scalar: Array of scalar values
# @param dx: Spacing between data points in X direction.
# @param dy: Spacing between data points in Y direction.
# @return: gradient of scalar
# @rtype: tuple(u, v)
def execute(scalar, dx, dy):
"Calculate the 2D gradient arrays of a 2+D scalar array."
result_u = empty(shape(scalar), dtype=scalar.dtype)
result_v = empty(shape(scalar), dtype=scalar.dtype)
# Left/rt edges of result_u can't be calculated.
result_u[:,0] = NaN
result_u[:,-1] = NaN
# Top/bot edges of result_v can't be calculated.
result_v[0,:] = NaN
result_v[-1,:] = NaN
# If dx and dy are arrays, remove extra cells.
shapedx = shape(dx)
if len(shapedx) < sum(shapedx):
dx = dx[:,1:-1]
shapedy = shape(dy)
if len(shapedy) < sum(shapedy):
dy = dy[1:-1,:]
# assume dx and dy are never zero
# dx = masked_values(dx, 0.0, copy=False)
# dy = masked_values(dy, 0.0, copy=False)
# calculate d(scalar)/dx
ans_u = scalar[:,2:] - scalar[:,0:-2]
ans_u = ans_u/(2 * dx)
# calculate d(scalar)/dy
ans_v = scalar[0:-2,:] - scalar[2:,:]
ans_v = ans_v/(2 * dy)
result_u[:,1:-1] = ans_u
result_v[1:-1,:] = ans_v
return Vector.componentsTo(result_u, result_v)

View file

@ -1,72 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# Returns heat index or wind chill
# History
# 12/02/2013 DR 14455 Qinglu Lin Changed 1.85200 to 3.6 at
# wSpd_kmh = wSpd * 1.85200
# ----------------------------------------------------------------
import numpy
import T
import DpT
import HeatIndex
import WindChill
from unit import knotToMetersPS
from unit import celciusToKelvin
# @param temperature: Temperature in degrees C
# @param tempFromTenths: tempFromTenths in degrees C
# @param dewpoint: dewpoint in degrees C
# @param tdpFromTenths: dpFromTenths in degrees C
# @param rwindSpeed: Wind Speed in knots
# @return: Heat Index or Wind Chill in degrees K
# @rtype: numpy array of Python floats or Python float
def execute1(temperature,tempFromTenths,dewpoint,dpFromTenths,windSpeed):
TK = T.execute1(temperature,tempFromTenths) #Outputs Kelvin
DpTK = DpT.execute4(dewpoint,dpFromTenths) #Outputs Kelvin
wSpd = knotToMetersPS(windSpeed)
return execute3(TK,DpTK,wSpd)
# @param temperature: Temperature in degrees K
# @param dewpoint: dewpoint in degrees K
# @param relHumidity: Relative humidity from 0 to 100
# @param windSpeed: Wind Speed in meter per second
# @return: Heat Index or Wind Chill in degrees K
# @rtype: numpy array of Python floats or Python float
def execute2(temperature,dewpoint,relHumidity,windSpeed):
DpTK = DpT.execute5(dewpoint,temperature,relHumidity) #Outputs Kelvin
return execute3(temperature,DpTK,windSpeed)
# @param T: Temperature in degrees K
# @param DpT: dewpoint in degrees K
# @param wSpd: Wind Speed in meter per second
# @return: Heat Index or Wind Chill in degrees K
# @rtype: numpy array of Python floats or Python float
def execute3(T,DpT,wSpd):
TC = T - 273.15 #convert from K to C
DpTC = DpT - 273.15 #convert to from K to C
wSpd_kmh = wSpd * 3.6 #convert from m/s to km/h
Hi = HeatIndex.calculate(TC,DpTC) #Outputs Celsius
Wc = WindChill.calculate(TC,wSpd_kmh) #Outputs Celsius
HiK = numpy.where(Hi != -9999.0, celciusToKelvin(Hi),-9999.0)
WcK = numpy.where(Wc != -9999.0, celciusToKelvin(Wc),-9999.0)
return numpy.where(HiK != -9999.0, HiK, numpy.where(WcK != -9999.0, WcK, -9999.0))

View file

@ -1,55 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# Calculates the heatIndex from temperature(C) and dewpoint temperatures(C)
# ----------------------------------------------------------------
import numpy
from numpy import power
const1 = 0.0091379024
const2 = 6106.396
def calculate(T, DpT):
# T and DpT have to be in C for the below calculation
badValue = numpy.where(T < 26.5,1,numpy.where(DpT > T,1,0))
TK = T + 273.15
DpTK = DpT + 273.15
# Legacy src/meteoLib/calcrh.f
exponent = numpy.where(T < 80.0, const1 * (TK - DpTK) + const2/TK - const2/DpTK,const1 * (T - DpT) + const2/T - const2/DpT)
RH = 100 * numpy.exp(exponent)
T1 = (T * 1.8) + 32
T1_sq = T1 * T1
RH_sq = RH * RH
# the Lans Rothfusz formula
H1 = -42.379 + (2.04901523 * T1) + (10.14333127 * RH)
H2 = (-0.22475541 * T1 * RH) - (0.00683783 * T1_sq) - (0.05481717 * RH_sq)
H3 = (0.00122874 * T1_sq * RH) + (0.00085282 * T1 * RH_sq)
H4 = (-0.00000199 * T1_sq * RH_sq)
Hi = H1 + H2 + H3 + H4
HiC = (Hi - 32) / 1.8 #convert F to C
#Returned in Celsius
return numpy.where(badValue == 1, -9999.0, HiC)

View file

@ -1,36 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
def execute(uStk, vStk, RM5):
umot,vmot = RM5
u1 = uStk[0]
v1 = vStk[0]
u2 = uStk[-1]
v2 = vStk[-1]
# First do our motion, lower bulk shear computation.
hptr = (v2-v1)*umot+(u1-u2)*vmot
for i in range(1, len(uStk)):
u1 = uStk[i-1]
v1 = vStk[i-1]
u2 = uStk[i]
v2 = vStk[i]
hptr += u2*v1-u1*v2
return hptr

View file

@ -1,34 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import power, piecewise
T0 = 288.0
gamma = 0.0065
p0 = 1013.2
c1 = 5.256
c2 = 14600
z11 = 11000
p11 = 226.0971
# Routine to calculate pressure from height based on a standard
# atmosphere
def execute(z):
return piecewise(z, [z < z11, z >= z11], [lambda z: p0*power((T0-gamma*z)/T0,c1), lambda z: p11*power(10, (z11-z)/c2)])

View file

@ -1,26 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import where, isnan
def execute(precip, accumNow, accumBefore):
if(len(precip.shape) == 2):
precip = precip.sum(1)
return where(isnan(precip), accumNow-accumBefore, precip)

View file

@ -1,87 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from scipy.interpolate import Rbf
from numpy import zeros, float32, NaN, isnan, array
##
# Designed to replace interp_up and interp_down in design files for point data
#
def execute(paramArray, vertArray, numLevels, vertPoints, maxGap=None):
ret = zeros(numLevels.shape, float32)
ret.fill(NaN)
for i in range(len(vertArray)):
verts = vertArray[i]
if type(verts) == float32:
verts = [verts]
params = paramArray[i]
# clone verts and params before modifying
verts = array(verts)
params = array(params)
gi = 0
for ci in range(len(verts)):
if isnan(verts[ci]) or isnan(params[ci]):
continue
verts[gi] = verts[ci]
params[gi] = params[ci]
gi += 1
if(gi == 0):
ret[i] = NaN
continue
verts = verts[:gi]
params = params[:gi]
if isinstance(vertPoints, float):
vertPoint = vertPoints
elif isinstance(vertPoints, float32):
vertPoint = vertPoints
else:
vertPoint = vertPoints[i]
if (maxGap != None):
below = None
above = None
gap = None
match = None
for point in verts:
if (point > vertPoint and (above == None or above > point)):
above = point
if (point < vertPoint and (below == None or below < point)):
below = point
if (point == vertPoint):
match = point
if(above == None and below == None):
continue
elif(above == None):
gap = (vertPoint - below)*2
elif(below == None):
gap = (above - vertPoint)*2
else:
gap = above - below
if((gap == None or gap > maxGap) and match == None):
ret[i] = NaN
continue
# If withion gap and only 1 value just use it
if(len(verts) == 1):
ret[i] = params[0]
continue
if(len(verts) == 1):
ret[i] = NaN
continue
rbf = Rbf(verts, params, function="linear")
ret[i] = rbf(vertPoint)
return ret

View file

@ -1,54 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file IsenStability.py
from numpy import NaN
from numpy import isnan
from numpy import maximum
from numpy import shape
##
# Calculate the isentropic stability through a layer.
#
# @param p_up: Pressure on upper isentrope (mb)
# @param p_lo: Pressure on lower isentrope (mb)
# @param o_up: Upper isentrope (K)
# @param o_lo: This (lower) isentrope (K)
# @return: Isentropic stability array through layer
# @rtype: numpy array
def execute(p_up, p_lo, o_up, o_lo):
"Calculate the isentropic stability through a layer."
# protect against divide-by-zero errors
o_diff = o_up - o_lo
dth = 1.0 / o_diff
# calculate stability [maximum() is elementwise]
p_diff = p_lo-p_up
stab = maximum(p_diff, 10.0) * dth
# maximum(NaN,X) is X. Restore the NaNs.
if shape(p_diff) == ():
if isnan(p_diff):
stab = NaN
else:
stab[isnan(p_diff)] = NaN
return stab

View file

@ -1,74 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file Laplacian.py
from numpy import empty
from numpy import shape
from numpy import NaN
##
# Function to calculate the Laplacian of qan.
#
# @param qan: The quantity array. This must be at rank 2, with at least 3
# cells in the first two dimensions.
# @param dx: The spacing between data points in the X direction.
# @param dy: The spacing between data points in the Y direction.
def execute(qan, dx, dy):
result = empty(shape(qan), dtype=qan.dtype)
# outer edges can't be calculated: fill with "invalid" value
result[0,:] = NaN
result[-1,:] = NaN
result[1:-1,0] = NaN
result[1:-1,-1] = NaN
# If dx and dy are matrices, we never use the outer edge,
# so strip it off so we don't have to use slice notation in the math.
# If they're actually scalars or 1-element matrices, we can't
# slice them so don't try.
shapedx = shape(dx)
if len(shapedx) < sum(shapedx):
dx = dx[1:-1,1:-1]
shapedy = shape(dy)
if len(shapedy) < sum(shapedy):
dy = dy[1:-1, 1:-1]
# Calculate the Q[i,j] + Q[i,j] array.
twoq = qan[1:-1,1:-1] * 2
# Calculate the X part of the answer
ans = qan[1:-1,0:-2] + qan[1:-1,2:]
ans -= twoq
ans /= dx * dx
# Calculate the Y part of the answer
term = qan[0:-2,1:-1] + qan[2:,1:-1]
term -= twoq
term /= dy * dy
# Combine pieces of the answer
ans += term
# fill middle block of result with ans
result[1:-1,1:-1] = ans
return result

View file

@ -1,51 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
## @file gamma.py
from numpy.ma.core import masked_values
from numpy.ma.core import log
##
# Calculate lapse rate from temperature and pressure pairs.
# This method is derived from the lapserate.f function, which
# used an extra integer parameter to determine which of 5 functions
# were to be used. However, only the "4" (pressure) function was
# being used.
#
# @param Tlo: Temperature at lower level (K)
# @param Plo: Pressure at lower level (mb)
# @param Thi: Temperature at upper level (K)
# @param Phi: Pressure at upper level (mb)
# @return: Lapse rate (C/km)
def execute(Tlo, Plo, Thi, Phi):
C = 0.034167
Tratio = Thi/Tlo
Pratio = Phi/Plo
logTratio = log(Tratio)
logPratio = log(Pratio)
# Guard against divide-by-zero errors, again
logPratio = masked_values(logPratio, 0, copy=False)
result = C * logTratio / logPratio
return result

View file

@ -1,118 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import copy
from numpy import exp
from numpy import isnan
from numpy import log
from numpy import maximum
from numpy import minimum
from numpy import shape
from numpy import sqrt
from numpy import NaN
from numpy import ndarray
import AdiabaticTemp
import TempOfTe
lowKTemp = 193
# Constants used in calculations.
# Wish the Fortran documented what these are.
const1 = 22.05565
const2 = 0.0091379024
const3 = 6106.396
const4 = 26.66082
const5 = 223.1986
const6 = 0.0182758048
const7 = -0.37329638
const8 = 41.178204
const9 = 0.0015945203
const10 = 3.498257
const11 = 0.286
import numpy as np
##
# Calculate lifted index.
#
# @param P: Pressure at current level (mb) [array]
# @param T: Temperature at current level (K) [array]
# @param RH: Relative humidity at current level (0-100) [array]
# @param T_500MB: Temperature at 500MB (K) [array]
# @param P_500MB: Pressure at 500MB (mb) [usually scalar==500]
# @return: Calculated lifted index array
def execute(P, T, RH, T_500MB, P_500MB=500):
"Calculate lifted index from P,T,RH,T at 500MB, P at 500MB"
# Mask any invalid inputs (even in scalars)
T = np.copy(T);
T[T < lowKTemp] = NaN
rhqc = minimum(100.0, maximum(1.0, RH))
# NaN cells in RH will be 1.0 in rhqc. Make them NaN again.
rhqc[isnan(RH)] = NaN
eee = const1 - const2 * T
# T is in degrees K, so it shouldn't ever be zero.
eee -= const3 / T
eee = exp(eee) * rhqc
b = const4 - log(eee)
tdp = b-sqrt(b*b-const5)
tdp /= const6
tc = tdp - (T-tdp)*(const7+const8/T+const9*tdp)
pc = P * (tc/T)**const10
# We need to process cells where pc <= P_500MB differently
# from the ones where it is > P_500MB. However, we don't want
# to do anything with NaN cells. Create two boolean masks for
# array indexing. They are logical inverses of one another,
# except that both contain False where cells in pc or P_500MB
# are NaN.
le_mask = pc <= P_500MB
gt_mask = pc > P_500MB
if shape(P_500MB) == shape(le_mask):
P5_le = P_500MB[le_mask]
P5_gt = P_500MB[gt_mask]
else: # assume it was a scalar
P5_le = P_500MB
P5_gt = P_500MB
# Create an output array with all cells NaN.
result = copy(T)
result[:] = NaN
intm_result = (P5_le/P)**const11
# make sure to properly mask an intm_result array
if type(intm_result) == ndarray:
intm_result = intm_result[le_mask]
intm_result *= -T[le_mask]
result[le_mask] = T_500MB[le_mask] + intm_result
result = result
tc_ge = AdiabaticTemp.calculate(tc[gt_mask], pc[gt_mask]) * (P5_gt/pc[gt_mask])**const11
result[gt_mask] = T_500MB[gt_mask] - TempOfTe.execute(tc_ge, P5_gt)
result[P < P_500MB] = NaN
return result

View file

@ -1,57 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from Add import execute as Add
from Multiply import execute as Multiply
from numpy import zeros_like, ndarray
def execute(*args):
"""Perform a linear transform
@return: result = arg1*arg2 + arg3*arg4 + arg5...
vecresult = vec1*sca1 + vec2*sca2 + vec3...
"""
result = None
if type(args[0]) == tuple:
zeroArray = zeros_like(args[0][0])
result = (zeroArray, zeroArray.copy())
else:
for arg in args:
if (type(arg) == ndarray):
result = zeros_like(arg)
break
termLength = 2
terms = len(args) / termLength
if len(args) % termLength != 0:
terms = terms + 1
for i in range(terms):
coefficient = args[i * termLength]
variable = 1 if i * termLength + 1 >= len(args) else args[i * termLength + 1]
result = Add(result, Multiply(coefficient, variable))
return result

View file

@ -1,75 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import zeros, float32, NaN, isnan
##
# Designed to replace interp_up and interp_down in design files for point data
#
def execute(paramArray, vertArray, numLevels, vertPoints, maxGap=None):
#print "paramArray =", paramArray
#print "vertArray =", vertArray
#print "numLevels =", numLevels
#print "vertPoints =", vertPoints
#print "maxGap =", maxGap
ret = zeros(numLevels.shape, 'float32')
ret.fill(NaN)
for i in range(len(vertArray)):
count = numLevels[i]
verts = vertArray[i][:count]
params = paramArray[i][:count]
if isinstance(vertPoints, float):
vertPoint = vertPoints
elif isinstance(vertPoints, float32):
vertPoint = vertPoints
else:
vertPoint = vertPoints[i]
below = None
above = None
for j in range(len(verts)):
if isnan(verts[j]) or isnan(params[j]):
continue
if (verts[j] > vertPoint and (above == None or verts[above] > verts[j])):
above = j
if (verts[j] < vertPoint and (below == None or verts[below] < verts[j])):
below = j
if(above == None and below == None):
ret[i] = NaN
continue
if (maxGap != None):
gap = None
if(above == None):
gap = (vertPoint - verts[below])*2
elif(below == None):
gap = (verts[above] - vertPoint)*2
else:
gap = verts[above] - verts[below]
if(gap == None or gap > maxGap):
ret[i] = NaN
continue
if(above == None):
ret[i] = params[below]
elif(below == None):
ret[i] = params[above]
else:
wgt = (vertPoint-verts[below])/(verts[above] - verts[below])
ret[i] = params[below] + wgt*(params[above] - params[below])
#print "ret = ", ret
#print "----------------------------"
return ret

View file

@ -1,47 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
#
#
# SOFTWARE HISTORY
#
# Date Ticket# Engineer Description
# ------------ ---------- ----------- --------------------------
# 01/19/12 DR14296 M. Huang Flipped result by applying -1
#
#
from numpy import power
from LvlQvec import calculate as lvlQvec
def execute(GHxSM, TxSM, P, dx, dy, coriolis):
slqx, slqy, dtdx, dtdy = lvlQvec(GHxSM, TxSM, P, dx, dy, coriolis)
# Compute the temperature to potential temperature ratio for this level.
t2th = power(1000/P, 0.286)
t2th *= 2
# Now calculate the QG frontogensis function for this level.
# Fgen = 2(Qx * d(theta)/dx + Qy * d(theta)/dy)
result = slqx * dtdx
result -= slqy * dtdy
result *= t2th
return result

View file

@ -1,104 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import log, exp
import PartialDerivative as Partial
import DgeoComps
import Vector
##
# Find Q vectors from height, temp, and pressure.
#
# @param height: Height (m)
# @type height: 2D numpy array
# @param temp: Temperature (K)
# @type temp: 2D numpy array
# @param pressure: Pressure (mb)
# @type pressure: scalar or 2D numpy array
# @param dx: Spacing between data points in the X direction
# @type dx: scalar or 2D numpy array
# @param dy: Spacing between data points in the Y direction
# @type dy: scalar or 2D numpy array
# @param coriolis: Coriolis force (kg-m/s^2)
# @type coriolis: 2D numpy array
# @return: Q vectors and dtemp/dx and dtemp/dy
# @rtype: tuple(U,V, dtdx, dtdy) of 2D masked numpy arrays
def execute(GHxSM, TxSM, P, dx, dy, coriolis):
result_u, result_v, dtdx, dtdy = calculate(GHxSM, TxSM, P, dx, dy, coriolis)
# convert the results we want to unmasked arrays
return Vector.componentsTo(result_u, result_v)
##
# Find Q vectors and dtemp/dx and dtemp/dy from height, temp, and pressure.
# This is an adaptation of slqvect.f. That function had
# multiple calls to a smoothing function, but I didn't find
# anything to set the number of passes to anything but 0, so
# the smoothing loops have been omitted; comments in the code
# mark where they would have been called.
#
# In slqvect.f, comments described dtdx and dtdy as work arrays, but
# slfront.f was treating them as output arrays since they were filled
# with the values slfront.f needed during slqvect.f's processing. This
# method was created to return all 4 arrays, so fGen.py wouldn't have to
# re-generate dtdx and dtdy.
#
# @attention: Returned vectors may contain NaN.
#
# @param height: Height (m)
# @type height: 2D numpy array
# @param temp: Temperature (K)
# @type temp: 2D numpy array
# @param pressure: Pressure (mb)
# @type pressure: scalar or 2D numpy array
# @param dx: Spacing between data points in the X direction
# @type dx: scalar or 2D numpy array
# @param dy: Spacing between data points in the Y direction
# @type dy: scalar or 2D numpy array
# @param coriolis: Coriolis force (kg-m/s^2)
# @type coriolis: 2D numpy array
# @return: Q vectors and dtemp/dx and dtemp/dy
# @rtype: tuple(U,V, dtdx, dtdy) of 2D masked numpy arrays
def calculate(height, temp, pressure, dx, dy, coriolis):
"Find Q vectors from height, temp, and pressure."
# assume dx, dy, and coriolis don't need masked
# Compute the temperature to potential temperature ratio for this level.
tptRatio = (1000/pressure) ** 0.286
#******* Code to smooth height omitted *******
# Compute the geostrophic wind components for this level.
dugdx, dugdy, dvgdx, dvgdy = DgeoComps.execute(height, dx, dy, coriolis )
#******* Code to smooth temp omitted *******
# Compute the derivative of temp with respect to X and Y
dtdx, dtdy = Partial.calculate(temp, dx, dy)
# Calculate the X and Y components of Q.
result_x = dugdx * dtdx
result_x += dvgdx * dtdy
result_x *= tptRatio
result_y = dugdy * dtdx
result_y += dvgdy * dtdy
result_y *= -tptRatio
return (result_x, result_y, dtdx, dtdy)

View file

@ -1,30 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from LyrQvec import calculate as lyrQvec
def execute(height_up, height_lo, pressure_up, pressure_lo, dx, dy, coriolis):
qx, qy, dtdx, dtdy = lyrQvec(height_up, height_lo, pressure_up, pressure_lo, dx, dy, coriolis)
fgen = qx * dtdx
fgen -= qy * dtdy
fgen *= 2
return fgen

View file

@ -1,125 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import log, exp
import PartialDerivative as Partial
import DgeoComps as DgeoComps
import Vector
##
# Find Q vectors from upper and lower height and pressure.
#
# @param height_up: Height at each grid point for the top of the layer (m)
# @type height_up: 2D numpy array
# @param height_lo: Height at each grid point for the bottom of the layer (m)
# @type height_lo: 2D numpy array
# @param pressure_up: Pressure level corresponding to height_up (mb)
# @type pressure_up: scalar or 2D numpy array
# @param pressure_lo: Pressure level corresponding to height_lo (mb)
# @type pressure_lo: scalar or 2D numpy array
# @param dx: Spacing between data points in the X direction
# @type dx: scalar or 2D numpy array
# @param dy: Spacing between data points in the Y direction
# @type dy: scalar or 2D numpy array
# @param coriolis: Coriolis force (kg-m/s^2)
# @type coriolis: 2D numpy array
# @return: Q vectors
# @rtype: tuple(U,V) of 2D numpy arrays
def execute(height_up, height_lo, pressure_up, pressure_lo, dx, dy, coriolis):
qx, qy, dtdx, dtdy = calculate(height_up, height_lo, pressure_up, pressure_lo, dx, dy, coriolis)
# unmask the arrays we're interested in
return Vector.componentsTo(qx, qy)
##
# Find Q vectors and dtdx and dtdy from upper and lower height and pressure.
#
# In qvector.f, comments described dtdx and dtdy as work arrays, but
# frontogen.f was treating them as output arrays since they were filled
# with the values frontogen.f needed during qvector.f's processing. This
# method was created to return all 4 arrays, so fGen.py wouldn't have to
# re-generate dtdx and dtdy.
#
# @attention: Returned vectors may contain NaN.
#
#
# @param height_up: Height at each grid point for the top of the layer (m)
# @type height_up: 2D numpy array
# @param height_lo: Height at each grid point for the bottom of the layer (m)
# @type height_lo: 2D numpy array
# @param pressure_up: Pressure level corresponding to height_up (mb)
# @type pressure_up: scalar or 2D numpy array
# @param pressure_lo: Pressure level corresponding to height_lo (mb)
# @type pressure_lo: scalar or 2D numpy array
# @param dx: Spacing between data points in the X direction
# @type dx: scalar or 2D numpy array
# @param dy: Spacing between data points in the Y direction
# @type dy: scalar or 2D numpy array
# @param coriolis: Coriolis force (kg-m/s^2)
# @type coriolis: 2D numpy array
# @return: Q vectors
# @rtype: tuple(U,V) of 2D numpy arrays
def calculate(height_up, height_lo, pressure_up, pressure_lo, dx, dy, coriolis):
"Find Q vectors from upper and lower height and pressure."
# assume dx, dy, and coriolis don't need masked
# Acceleration due to gravity (m/s**2)
gravAcc = 9.806
# magic numbers used in calculations.
# Taken from original Fortran code; not sure what they are.
magic_exp = 0.286
magic_fac = 287
height_mid = height_up + height_lo
height_mid /= 2
#******* Code to smooth height omitted *******
# calculate components of geostrophic wind
dugdx, dugdy, dvgdx, dvgdy = DgeoComps.execute(height_mid, dx, dy, coriolis)
# get mean pressure
# copied from Fortran; wouldn't sqrt(p_up * p_lo) be better?
pavg = log(pressure_up) + log(pressure_lo)
pavg /= 2
pavg = exp( pavg )
# derive cnvFac to convert potential temp to thickness
pdiff = pressure_up - pressure_lo
gravAcc = 9.806
denom = (pavg/1000) ** magic_exp
denom *= pdiff * magic_fac
cnvFac = -gravAcc * pavg / denom
temp = cnvFac * (height_up - height_lo)
dtdx, dtdy = Partial.execute(temp, dx, dy)
qx = -dugdx * dtdx
qx -= dvgdx * dtdy
qy = -dugdy * dtdx
qy -= dvgdy * dtdy
return (-qx, qy, dtdx, dtdy)

View file

@ -1,41 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
#Returns the sum of elevation and the 1st skyLayerBase
#
# ----------------------------------------------------------------
import numpy
##
# Returns the sum of elevation and the 1st skyLayerBase
#
# @param skyLayerBase, : skyLayerBase,
# @param elevation : elevation
# @return:
# @rtype:
#
def execute(skyLayerBase,elevation, index):
#convert elevation (m) to elevation (ft)
elevation_ft = (elevation * 3.2808399)
return numpy.where(skyLayerBase[:,index] != -9999.0, skyLayerBase[:,index] + elevation_ft, -9999.0)

View file

@ -1,34 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# Calculate magnitude
# ----------------------------------------------------------------
from numpy import hypot
def execute(*args):
""" Return the magnitude of a given vector or vector components
"""
if type(args[0]) == tuple:
return hypot(args[0][0], args[0][1])
else:
return hypot(args[0], args[1])

View file

@ -1,38 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
import numpy
from numpy import zeros
def execute(values, *mappings):
tmp = zeros(values.shape, 'float32')
size = len(mappings)
i = 0
while i < size:
checkFor = mappings[i]
i += 1
setTo = mappings[i]
i += 1
tmp[values == checkFor] = setTo
return tmp

View file

@ -1,67 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# SOFTWARE HISTORY
#
# Date Ticket# Engineer Description
# ------------ ---------- ----------- --------------------------
# 18 Feb 2010 #4502 jelkins Initial Creation.
from numpy import array, nanmax
def execute(*args):
""" Return the max value at each grid point
@param args
a list of grids or a single 3D grid. If a single 3D variable is passed in,
will compute without considering the vertical coordinate information.
"""
if len(args) == 1 and isinstance(args[0], list):
if len(args[0]) == 2 and len(args[0][0].shape) == 3:
# we have been passed 3D data
return nanmax(args[0][0], 0)
else:
return execute(*args[0])
else:
return nanmax(array(list(args)), 0)
def test():
""" Unit test
"""
from numpy import all
threeDarray = array([[[1.,2.],[20.,3.]],[[3.,40.],[4.,5.]],[[5.,6.],[6.,7.]]])
grid1 = threeDarray[0]
grid2 = threeDarray[1]
grid3 = threeDarray[2]
correctResult = array([[ 5, 40],[20, 7]])
result = execute([threeDarray,0])
if not(all(result == correctResult)):
raise Exception
result = execute(grid1,grid2,grid3)
if not(all(result == correctResult)):
raise Exception
print "Max Test Complete"

View file

@ -1,26 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import where, isnan, NaN
def execute(precip, accum_precip):
sum = accum_precip.sum(axis=1)
sum[sum < -0.005] = NaN
return where(isnan(precip), sum, precip)

View file

@ -1,44 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# SOFTWARE HISTORY
#
# Date Ticket# Engineer Description
# ------------ ---------- ----------- --------------------------
# 18 Feb 2010 #4502 jelkins Initial Creation.
from numpy import array, nanmin
def execute(*args):
""" Return the min value at each grid point
@param args
a list of grids or a single 3D grid. If a single 3D variable is passed in,
will compute without considering the vertical coordinate information.
"""
if len(args) == 1 and isinstance(args[0], list):
if len(args[0]) == 2 and len(args[0][0].shape) == 3:
# we have been passed 3D data
return nanmin(args[0][0], 0)
else:
return execute(*args[0])
else:
return nanmin(array(list(args)), 0)

View file

@ -1,39 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import power
##
# This routine estimates 1000 to 500 mb thickness from 500 height
# and mean sea level pressure.
def execute(mslp, GH):
"Calculate thickness from mean sea level pressure and height."
# Constants from Fortran code.
# Wish these had more descriptive names; don't know what they are.
a = 0.4599042
b = 3.262312
c = 0.1902672
# calculate thickness
dZ = GH * a
denom = power(mslp, c)
denom -= b
dZ /= denom
return dZ

View file

@ -1,39 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import multiply
import Vector
def execute(*args):
""" Perform multiplication of any number of scalars or of a vector and a scalar.
"""
if type(args[0]) == tuple:
return vectorMultiply(args)
else:
return scalarMultiply(args)
def scalarMultiply(args):
return reduce(multiply, args)
def vectorMultiply(args):
return Vector.componentsTo(scalarMultiply((args[0][0], args[1])), scalarMultiply((args[0][1], args[1])))

View file

@ -1,74 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
#import sys
#sys.path.append("..")
from numpy import ndarray,NaN, float32
from Multiply import execute as Multiply
import Vector
##
# This routine computes the non-advective local change of an arbitrary
# C conservative parameter `a'.
#
# @param u: U component of wind.
# @param v: V component of wind.
# @param a: Arbitrary conservative parameter.
# @param dx: Grid spacing in the x-direction (m)
# @param dy: Grid spacing in the y-direction (m)
# @return: A vector representing the non-advective local change
def execute(u, v, a, dx, dy):
u = -u
# TODO: support scalars for dx and dy
# find 2x partial derivatives of u, v, and a with respect to x
dudy = (u[2:,1:-1] - u[0:-2,1:-1]) / dy[1:-1,1:-1]
dvdy = (v[2:,1:-1] - v[0:-2,1:-1]) / dy[1:-1,1:-1]
dady = (a[2:,1:-1] - a[0:-2,1:-1]) / dy[1:-1,1:-1]
# find 2x partial derivatives of u, v, and a with respect to y
dudx = (u[1:-1,2:] - u[1:-1,0:-2]) / dx[1:-1,1:-1]
dvdx = (v[1:-1,2:] - v[1:-1,0:-2]) / dx[1:-1,1:-1]
dadx = (a[1:-1,2:] - a[1:-1,0:-2]) / dx[1:-1,1:-1]
# create empty arrays for output
dadxdt = ndarray(u.shape, float32)
dadydt = ndarray(u.shape, float32)
# First output array:
# mark the output array edges as invalid
dadxdt[0,:] = NaN
dadxdt[-1,:] = NaN
dadxdt[1:-1,0] = NaN
dadxdt[1:-1,-1] = NaN
# calculate values for the middle
dadxdt[1:-1,1:-1] = -0.5 * (dudx*dadx + dvdx*dady)
# Second output array:
# mark the output array edges as invalid
dadydt[0,:] = NaN
dadydt[-1,:] = NaN
dadydt[1:-1,0] = NaN
dadydt[1:-1,-1] = NaN
# calculate values for the middle
dadydt[1:-1,1:-1] = -0.5 * (dudy*dadx + dvdy*dady)
return Vector.componentsTo(-dadxdt, dadydt)

View file

@ -1,30 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file Negate.py
import Vector
##
# Negates input
def execute(input):
if type(input) == tuple:
return Vector.componentTo(-input[0], -input[1])
else:
return -input

View file

@ -1,73 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
import gridslice
import meteolib
import numpy
from numpy import equal
from numpy import where
from numpy import zeros
from numpy import empty
from numpy import concatenate
from numpy import greater
from numpy import less
from unit import pascalToMilliBar
import gridslice
def execute1(pressure, sfcPress):
sfcPress = sfcPress.reshape(- 1, 1);
result = concatenate((sfcPress, pressure), 1)
return pascalToMilliBar(result)
def execute2(levels, staElev):
staElev = staElev.reshape(- 1, 1);
hft2m = 30.48
levelsInMeters = where(equal(levels, - 9999), - 9999, levels * 30.48)
GH = concatenate((staElev, levelsInMeters), 1)
return meteolib.ztopsa(GH)
def execute3(numProfLvlsStation, MB):
ret = zeros(numProfLvlsStation.shape, 'float32')
ret.fill(MB)
return ret
def execute4(prCloudStation,lowCldStation,midCldStation,hiCldStation):
prCloudMB = prCloudStation/100
CCPval = ccpExecute(lowCldStation,midCldStation,hiCldStation)
isCeiling = where(greater(CCPval, 0.5), CCPval, -9999)
prCloudMB[isCeiling == -9999] = -9999
prCloudMB[prCloudStation == -9999] = -9999
return prCloudMB;
def execute5(height, elevation):
# array height is in meters
# scalar elevation is in meters
pressure = where(equal(height, - 9999), - 9999, height + elevation)
pressure = meteolib.ztopsa(pressure)
return pressure
def execute6(numLevelsStation, MB):
ret = zeros(numLevelsStation.shape, 'float32')
ret.fill(MB)
return ret
def ccpExecute(lowCldStation,midCldStation,hiCldStation):
P = maximum(lowCldStation,midCldStation,hiCldStation)
return where(equal(P, -9999), -9999, P/100);

View file

@ -1,74 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
import numpy
from numpy import equal
from numpy import where
def execute1(precipTypeStation,POP_hour_bestCatStation):
tmp = where(equal(POP_hour_bestCatStation, 1), precipTypeStation, -9999)
tmp[tmp == 0] = 0
tmp[tmp == 1] = 71
tmp[tmp == 2] = 89
tmp[tmp == 3] = 79
return tmp.astype('float32')
def execute2(CRAIN, CFRZR, CICEP, CSNOW):
tmp = 1.0*CRAIN+2.0*CFRZR+4.0*CICEP+8.0*CSNOW
# convert numeric value into symbol
tmp[tmp <= 0] = 0
tmp[tmp > 9] = 0
tmp[tmp == 1] = 79
tmp[tmp == 2] = 89
tmp[tmp == 3] = 71
tmp[tmp == 4] = 77
tmp[tmp == 5] = 80
tmp[tmp == 6] = 67
tmp[tmp == 7] = 69
tmp[tmp == 8] = 47
tmp[tmp == 9] = 183
return tmp
def execute3(CRAIN6hr, CFRZR6hr, CICEP6hr, CSNOW6hr):
return execute2(CRAIN6hr, CFRZR6hr, CICEP6hr, CSNOW6hr)
def execute4(CRAIN12hr, CFRZR12hr, CICEP12hr, CSNOW12hr):
return execute2(CCRAIN12hr, CFRZR12hr, CICEP12hr, CSNOW12hr)
def execute5(CPOZP3hr, CPOFP3hr):
# probablilities are represented between 0 and 1
tmp= 2.0*CPOZP3hr+8.0*CPOFP3hr
tmp[tmp <= 0] = 0
tmp[tmp > 15] = 0
tmp[tmp == 1] = 79
tmp[tmp == 2] = 71
tmp[tmp == 3] = 71
tmp[tmp == 4] = 80
tmp[tmp == 5] = 183
tmp[tmp == 6] = 185
tmp[tmp == 7] = 185
tmp[tmp == 8] = 89
tmp[tmp == 9] = 47
tmp[tmp == 10] = 186
tmp[tmp == 11] = 186
tmp[tmp == 12] = 184
tmp[tmp == 13] = 188
tmp[tmp == 14] = 187
tmp[tmp == 15] = 187
return tmp

View file

@ -1,30 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from Average import execute as Average
from numpy import concatenate, where, zeros, shape, nan
import gridslice
def execute1(pvv):
pv = zeros((shape(pvv)[0], 1))
pv.fill(nan)
result = concatenate((pv, pvv),1)
return result

View file

@ -1,88 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
## @file PartialDerivative.py
from numpy import isscalar
from numpy import NaN
##
# Calculate d/dx and d/dy of Qty.
# This method assumes masked arrays as input and returns masked
# arrays as output. That saves work when the return value would just
# be masked again for use in other calculations. When in doubt, call
# execute() instead.
#
# The returned arrays are the same size as Qty, with all the outermost
# edges masked. d/dx could extend to the top and bottom, and
# d/dy to the left and right edge, but this has not been done, in order
# to maintain the same output as the old g2gkinematics.f function #7.
#
# @param Qty: Quantity to differentiate.
# @type Qty: 2D numpy array
# @param dx: Spacing between data points in X direction.
# @type dx: scalar or 2D numpy array
# @param dy: Spacing between data points in Y direction.
# @type dy: scalar or 2D numpy array
# @return: d/dx and d/dy of Qty.
# @rtype: tuple(dq/dx,dq/dy) of masked 2D numpy arrays.
def calculate(Qty, dx, dy):
"Calculate d/dx and d/dy of Qty."
if not isscalar(dx):
dx = dx[1:-1,1:-1]
if not isscalar(dy):
dy = dy[1:-1,1:-1]
cropped_dqdx = Qty[1:-1,2:] - Qty[1:-1,0:-2]
cropped_dqdx /= dx * 2
cropped_dqdy = Qty[2:,1:-1] - Qty[0:-2,1:-1]
cropped_dqdy /= dy * 2
dqdx = Qty + NaN
dqdy = Qty + NaN
dqdx[1:-1,1:-1] = cropped_dqdx
dqdy[1:-1,1:-1] = cropped_dqdy
return (dqdx, dqdy)
##
# Calculate d/dx and d/dy of Qty.
# This method takes unmasked arrays and returns unmasked arrays.
#
# @param Qty: Quantity to differentiate.
# @type Qty: 2D numpy array
# @param dx: Spacing between data points in X direction.
# @type dx: scalar or 2D numpy array
# @param dy: Spacing between data points in Y direction.
# @type dy: scalar or 2D numpy array
# @return: d/dx and d/dy of Qty.
# @rtype: tuple(dq/dx,dq/dy) of 2D numpy arrays.
def execute(Qty, dx, dy):
""
# assume dx and dy are never zero or near-inifinite
dqdx, dqdy = calculate(Qty, dx, dy)
return (dqdx, dqdy)

View file

@ -1,39 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import power
##
# Calculate theta from temperature and pressure.
# This function accepts numpy arrays
#
# @param P: Pressure in millibars
# @param T: Temperature in degrees Kelvin
# @return: Potential temperature in degrees Kelvin
#
def execute(P, T):
"Calculate theta (PoT) from temperature and pressure."
pComponent = 1000.0/P
pComponent = power(pComponent, 0.286)
theta = T * pComponent
# convert theta to an unmasked array with 1e37 for invalid values
return theta

View file

@ -1,60 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import power
from numpy import multiply
def execute(*args):
"""Combine the arguments into a polynomial result
@return result of arg1*arg2^arg3 + arg4*arg5^arg6 + arg7*arg8
"""
result = 0
terms = len(args) / 3
if len(args) % 3 != 0:
terms = terms + 1
for i in range(terms):
coefficient = args[i * 3]
variable = 1 if i * 3 + 1 >= len(args) else args[i * 3 + 1]
exponent = 1 if i * 3 + 2 >= len(args) else args[i * 3 + 2]
result += multiply(coefficient, power(variable, exponent))
return result
def test():
from numpy import array
if not(all(execute(array([1., 2.])) == array([1., 2.]))):
raise Exception
if not(all(execute(array([1., 2.]), array([3., 4.])) == array([3., 8.]))):
raise Exception
if not(all(execute(array([1., 2.]), array([3., 4.]), array([5., 6.])) == array([243., 8192.]))):
raise Exception
print "Poly Test Complete"

View file

@ -1,69 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
import IsenStability
import Vorticity
##
# Calculate the isentropic potential vorticity through a layer.
#
# @change: Adapted from calcpv.f 2008-18-06
# User Notes:
#
# 1. Stability is defined as -dP/d(theta). We calculate this through
# the layer from the isentropic surface 'n' to the surface above it,
# 'n+1'.
#
# 2. Since we are dealing with a layer, we calculate a mean absolute
# vorticity using the winds at the upper and lower layers.
#
# 3. The PV is then [mean abs. vort.]/[stability]
#
# @param p_up: Pressure on upper isentrope (mb)
# @param p_lo: Pressure on lower isentrope (mb)
# @param o_up: Upper isentrope (K) (usually scalar)
# @param o_lo: This (lower) isentrope (K) (usually scalar)
# @param Wind_up: tuple(U,V) of winds on upper isentrope (m/s)
# @param Wind_lo: tuple(U,V) of winds on lower isentrope (m/s)
# @param dx: Spacing between data points in X direction (m)
# @param dy: Spacing between data points in Y direction (m)
# @param coriolis: Coriolis parameters (/s)
# @return: isentropic potential vorticity array
#
def execute(p_up, p_lo, o_up, o_lo, vector_up, vector_lo, dx, dy, coriolis):
"Calculate the isentropic potential vorticity through a layer."
u_up, v_up = vector_up
u_lo, v_lo = vector_lo
# Calculate the absolute vorticity at each isentropic surface.
avort1 = Vorticity.execute(u_up, v_up, coriolis, dx, dy)
avort2 = Vorticity.execute(u_lo, v_lo, coriolis, dx, dy)
# Calculate the isentropic stability through the layer.
pvort = IsenStability.execute(p_up, p_lo, o_up, o_lo)
# Calculate isentropic potential vorticity
result = avort1 + avort2
result *= 0.5
result /= pvort
return result

View file

@ -1,87 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
#
# Software History
#
# 2013/1/17 DR 15655 Melissa Porricelli Modified final 'result'
# calculation to remove multiplication
# by 0.5. Displayed values were
# off by a factor of this amount
# in comparison to A1.
# A1 calc in pvpres.f.
###
from numpy import zeros
import Gradient
import Vorticity
##
# Calculate the isobaric potential vorticity through a layer.
#
# User Notes:
#
# 1. Stability is defined as -dP/d(theta). We calculate this through
# the layer from the isobaric surface 'n' to the surface above it,
# 'n+1'.
# 2. Since we are dealing with a layer, we calculate a mean absolute
# vorticity using the winds at the upper and lower layers.
# 3. The PV is then [mean abs. vort]/[stability] + theta->pres term
#
# originally from pvpres.f, by J.Ramer
# @change: Converted from Fortran on 2008-16-06
#
# @param t_up: Theta on upper isobaric sfc (K)
# @param t:lo: Theta on this isobaric sfc (K)
# @param p_up: Upper pressure (mb)
# @param p_lo: This (lower) pressure (mb)
# @param Wind_up: tuple(U,V) winds on upper surface (m/s)
# @param Wind_lo: tuple(U,V) winds on lower surface (m/s)
# @param dx: Spacing in X direction (m)
# @param dy: Spacing in Y direction (m)
# @param coriolis: Coriolis parameters (/s)
# @return: Isobaric potential vorticity
# @rtype: numpy array
def execute(t_up, t_lo, p_up, p_lo, vector_up, vector_lo, dx, dy, coriolis):
""
u_up, v_up = vector_up
u_lo, v_lo = vector_lo
# Calculate the absolute vorticity at each isobaric surface.
avort1 = Vorticity.execute(u_up, v_up, coriolis, dx, dy)
avort2 = Vorticity.execute(u_lo, v_lo, coriolis, dx, dy)
# Calculate the temperature gradient on each surface.
grad_lo = Gradient.execute(t_lo, dx, dy)
grad_up = Gradient.execute(t_up, dx, dy)
dtdx1, dtdy1 = grad_lo
dtdx2, dtdy2 = grad_up
# Calculate difference arrays.
dp = p_up - p_lo
dt = t_up - t_lo
du = u_up - u_lo
dv = v_up - v_lo
dtdx = dtdx1 + dtdx2
dtdy = dtdy1 + dtdy2
av = avort1 + avort2
result = (-0.5 * (av*dt + (du*dtdy - dv*dtdx)) / dp)
return result

View file

@ -1,28 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import power
def execute(arg1, arg2):
""" Raise the first arg to the second arg.
"""
return power(arg1,arg2)

View file

@ -1,26 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import where, NaN
# For modelsounding cloud layers, return prCloud when it is > bottom and < top
def execute(prCloudHgt, bottom, top):
prCloudHgt = where(prCloudHgt>bottom, prCloudHgt, NaN);
return where(prCloudHgt<top, prCloudHgt, NaN)

View file

@ -1,26 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import where, NaN
##
# Designed to replace sample_up and sample_down in design files for point data
#
def execute(data, qc):
return where(qc == 0, data, NaN)

View file

@ -1,53 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import exp
def execute1(T, DpT):
"Calculate Relative Humidity from Temperature and Dewpoint Temp."
RH = T-DpT
RH *= 0.0091379024
RH += 6106.396/T
RH -= 6106.396/DpT
RH = exp(RH)
RH *= 100
return RH
def execute2(P, T, SHx):
"Calculate Relative Humidity from Pressure, Temp, and Spec Humidity."
# Constants from the Fortran code.
a = 22.05565
b = 0.0091379024
c = 6106.396
epsilonx1k = 622.0
shxDenom = SHx * 0.378
shxDenom += epsilonx1k
tDenom = -b*T
tDenom += a
tDenom -= c/T
RH = P * SHx
RH /= shxDenom
RH /= exp(tDenom)
return RH

View file

@ -1,135 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from math import log, exp
from numpy import NaN, isnan
class Level:
def __init__(self, pr, ht, tp, td, wd, ws):
self.pr = pr
self.ht = ht
self.tp = tp
self.td = td
self.wd = wd
self.ws = ws
def getPr(self):
return self.pr
def getHt(self):
return self.ht
def getParam(self, index):
if(index == 1):
value = self.pr
elif(index == 2):
value = self.ht
elif(index == 3):
value = self.tp
elif(index == 4):
value = self.td
elif(index == 5):
value = self.wd
elif(index == 6):
value = self.ws
return value
def execute(prMan, htMan, tpMan, tdMan, wdMan, wsMan, prSigT, tpSigT, tdSigT, htSigW, wdSigW, wsSigW, numMand, numSigT, numSigW, validTime, wmoStaNum):
numMand[numMand <= -9999] = 0
numSigT[numSigT <= -9999] = 0
numSigW[numSigW <= -9999] = 0
levelsMap = {}
for c in range(0, len(numMand)):
key = (validTime[c],wmoStaNum[c])
levels = levelsMap.get(key,[])
levelsMap[key] = levels
for i in range(0, numMand[c]):
levels.append(Level(prMan[c][i], htMan[c][i], tpMan[c][i], tdMan[c][i], wdMan[c][i], wsMan[c][i]))
for i in range(0, numSigT[c]):
levels.append(Level(prSigT[c][i], NaN, tpSigT[c][i], tdSigT[c][i], NaN, NaN))
for levels in levelsMap.values():
levels.sort(key=Level.getPr)
for i in range(1,len(levels)):
level = levels[i]
# Interpolate missing Heights
if isnan(level.ht):
levelblo = levels[i-1]
tv1 = level.tp
tv2 = levelblo.tp
mtv = ((tv1 + tv2) / 2);
delta = log(levelblo.pr) - log(level.pr);
delta *= 29.2898 * mtv;
level.ht = levelblo.ht + delta
for c in range(0, len(numMand)):
key = (validTime[c],wmoStaNum[c])
levels = levelsMap.get(key,[])
levelsMap[key] = levels
for i in range(0, numSigW[c]):
levels.append(Level(NaN, htSigW[c][i], NaN, NaN, wdSigW[c][i], wsSigW[c][i]))
for levels in levelsMap.values():
levels.sort(key=Level.getPr)
levelsCopy = levels[:]
# Remove Pressure levels
for i in range(1,len(levels)):
if levelsCopy[i].pr == levelsCopy[i-1].pr:
levels.remove(levelsCopy[i])
levels.sort(key=Level.getHt)
levelsCopy = levels[:]
# Remove duplicate height levels
for i in range(1,len(levels)):
if levelsCopy[i].ht == levelsCopy[i-1].ht:
levels.remove(levelsCopy[i])
for i in range(1,len(levels)):
level = levels[i]
levelblo = levels[i-1]
# Interpolate missing Temps
if isnan(level.tp):
levelabv = None
for j in range(i,len(levels)):
if not(isnan(levels[j].tp)):
levelabv = levels[j]
break
if levelabv != None:
level.tp = levelblo.tp + (level.ht - levelblo.ht)*(levelabv.tp - levelblo.tp)/(levelabv.ht - levelblo.ht)
# Interpolate missing Dewpoints
if isnan(level.td):
levelabv = None
for j in range(i,len(levels)):
if not(isnan(levels[j].td)):
levelabv = levels[j]
break
if levelabv != None:
level.td = levelblo.td + (level.ht - levelblo.ht)*(levelabv.td - levelblo.td)/(levelabv.ht - levelblo.ht)
# Interpolate missing Pressure
if isnan(level.pr):
zA = levelblo.ht;
zB = level.ht;
dz = zB - zA;
tvA = levelblo.tp + 273.13;
tvB = level.tp + 273.13;
k = (tvA + tvB) / 2 * 28.2898;
pA = levelblo.pr;
k = (-dz / k) + log(pA);
level.pr = exp(k);
levelsList = []
for c in range(0, len(numMand)):
key = (validTime[c],wmoStaNum[c])
levels = levelsMap.get(key, None)
if (levels != None):
levelsList.append(levels)
return levelsList

View file

@ -1,37 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import array, NaN, float32
def execute(levelsList, index):
maxLen = 0
results = []
for i in range(len(levelsList)):
levels = levelsList[i]
values = []
for level in levels:
values.append(level.getParam(index))
if len(levels) > maxLen:
maxLen = len(levels)
results.append(values)
for result in results:
while len(result) < maxLen:
result.append(NaN)
return array(results, float32)

View file

@ -1,69 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from Add import execute as Add
from math import radians,cos,sin
from numpy import array,dot
import Vector
def execute(*args):
"""Rotate a vector
Rotate vector by a number of degrees, or transform the vector with a matrix.
The arguments after the vector need to be constants.
"""
if type(args[0]) != tuple:
raise ValueError("first argument must be a vector")
elif len(args) == 2:
return rotateDegrees(args)
else:
return vectorTransformMatrix(args)
def rotateDegrees(args):
v = args[0]
r = radians(args[1])
s = sin(r)
c = cos(r)
return vectorTransformMatrix(v,c,-s,s,c)
def vectorTransformMatrix(args):
uComponent,vComponent = args[0]
uComponentShape = uComponent.shape
vComponentShape = vComponent.shape
# in-place flatten the arrays
uComponent = uComponent.reshape((uComponent.size,))
vComponent = vComponent.reshape((vComponent.size,))
vector = array([uComponent, vComponent])
transform = array([[args[1], args[2]], [args[3], args[4]]],dtype=uComponent.dtype)
u,v = dot(transform, vector)
# resize the arrays back to appropriate size
u.resize(uComponentShape)
v.resize(vComponentShape)
return Vector.componentsTo(u, v)

View file

@ -1,33 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# Calculate Specific Humidity (g/Kg) from Pressure, Temperature and
# Relative Humidity.
# ----------------------------------------------------------------
from numpy import concatenate
from numpy import where
from numpy import equal
def execute(specHum, q2):
q2 = q2.reshape(-1, 1);
result = concatenate((q2, specHum), 1)
return where(equal(result, -9999), -9999, result*1000)

View file

@ -1,34 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import add
from numpy import array
from numpy import zeros
from Magnitude import execute as Magnitude
def execute(uStk, vStk):
res = zeros(uStk[0].shape, 'float32')
for i in range(1, len(uStk)):
u1 = uStk[i-1]
v1 = vStk[i-1]
u2 = uStk[i]
v2 = vStk[i]
res += Magnitude(u2-u1, v2-v1)
return res

View file

@ -1,67 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
import gridslice
from numpy import ndarray
def execute(*args):
#defineNumpySlice(vc, param, targetLevel, sense)
# vc - cube/list for output parameter
# param - cube of data that is the same as targetLevel
# targetLevel - constant value level
# sense - corresponds to vc for vertical coordinate and whether to search for highest/lowest occurence
#
#createNumpySlice(vc, s3d, targetLevel, sense, hyb-optional)
# vc - cube/list of data that is the same as targetLevel
# s3d - cube for output parameter
# targetLevel - 2d grid of data values to find in vc to interpolate s3d
# sense - corresponds to vc for vertical coordinate and whether to search for highest/lowest occurence
if len(args) == 3:
# cube, target level, sense
if type(args[1]) == ndarray:
# target level is 2d grid
rval = gridslice.createNumpySlice(args[0][1], args[0][0], args[1], int(args[2]))
return rval
else:
# target level is single value
rval = gridslice.defineNumpySlice(args[0][1], args[0][0], args[1], int(args[2]))
return rval
else:
if type(args[2]) == ndarray:
# cube, cube, grid, sense
rval = gridslice.createNumpySlice(args[0][0], args[1][0], args[2], int(args[3]))
return rval
else:
# cube, cube, level, sense
# need to make 2 slice calls, one to create the common level slice
# and a second to make the output slice
levelSense = -1
if abs(int(args[3])) > 1:
levelSense = -2
# create slice at common cube level
rval = gridslice.defineNumpySlice(args[0][1], args[0][0], args[2], levelSense)
rval = gridslice.createNumpySlice(args[1][1], args[1][0], rval, int(args[3]))
return rval

View file

@ -1,43 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
import gridslice
from numpy import zeros, float32
def execute(*args):
#createNumpySlice(vc, s3d, targetLevel, sense, hyb-optional)
# vc - cube/list of data that is the same as targetLevel
# s3d - cube for output parameter
# targetLevel - 2d grid of data values to find in vc to interpolate s3d
# sense - corresponds to vc for vertical coordinate and whether to search for highest/lowest occurence
if len(args) == 4:
# cube, grid, sense, hybrid
grid = args[1]
if isinstance(grid, float):
grid = zeros((args[0][1].shape[1], args[0][1].shape[2]),float32)
grid.fill(args[1])
rval = gridslice.createNumpySlice(args[0][1], args[0][0], grid, int(args[2]), int(args[3]))
return rval
else:
# cube, cube, grid, sense, hybrid
rval = gridslice.createNumpySlice(args[1][0], args[0][0], args[2], int(args[3]), int(args[4]))
return rval

View file

@ -1,43 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
# ----------------------------------------------------------------
#
#
# ----------------------------------------------------------------
import numpy
import U
import V
from numpy import hypot
def execute(windSpeed, windDir, accum_windSpeed24, accum_windDir24):
U0 = U.calculate(windSpeed, windDir)
V0 = V.calculate(windSpeed, windDir)
U24 = U.calculate(accum_windSpeed24, accum_windDir24)
V24 = V.calculate(accum_windSpeed24, accum_windDir24)
DU = U0 - U24
DV = V0 - V24
wSp = hypot(DU,DV)
return wSp

View file

@ -1,59 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# Calculate Specific Humidity (g/Kg) from Pressure, Temperature and
# Relative Humidity.
# ----------------------------------------------------------------
from numpy import exp
from numpy.ma.core import masked_greater
from numpy.ma.core import masked_values
from numpy.ma.core import filled
def execute(*args):
if len(args) == 2:
"Calculate specific humidity (SHx) from pressure and vapor pressure"
# From the old C version:
# formula is 1000*epsilon*e/(p-(1-epsilon)*e), where
# p is pressure, e is vapor pressure, and epsilon is the
# constant 0.622. For evaluation, 622.0 (1000*epsilon) is
# divided into the constants in the denominator.
# p is input field 0, e is input field 1. The 100.0 in
# the first constant is converting mb to pascals.
P=args[0]
VPP=args[1]
oneMinusEpsilon = 0.378
epsilonTimes1000 = 622
num = 100.0/epsilonTimes1000 * VPP
dnm = P - ((oneMinusEpsilon/epsilonTimes1000) * VPP)
SHx = num / dnm
return SHx
else:
"Calculate specific humidity (SHx) from pressure, temperature, and \
relative humidity."
# Get masked arrays so results for out-of-range values aren't generated
Pma = args[0]
Tma = args[1]
RHma = args[2]
# do the math
eee = RHma*exp(28.48859-0.0091379024*Tma-6106.396/Tma)
SHx = eee/(Pma-0.00060771703*eee)
return SHx

View file

@ -1,26 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import std, dstack
def execute(*args):
if len(args) == 1 and isinstance(args[0], list):
return execute(*args[0])
return dstack(args).std(axis=2)

View file

@ -1,27 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import where, NaN, float32
def execute(arg1, arg2):
if(arg2.dtype == float32):
return where(arg1 >= 0, arg2, NaN)
else:
return where(arg1 >= 0, arg2, -9999)

View file

@ -1,55 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
from numpy import sqrt
from numpy import ndarray
from numpy import float32
##
# Routine to calculate sweat index from the total totals,
# 850 dewpoint, and wind components at 850 and 500.
#
# @param tt: Total Totals
# @param td8: 850mb dewpoint (K)
# @param u8: 850mb u-component (K)
# @param v8: 850mb v-component (K)
# @param u5: 500mb u-component (K)
# @param v5: 500mb v-component (K)
# @return: Sweat Index
def execute(tt, td8, u8, v8, u5, v5):
s8 = sqrt(u8*u8+v8*v8)
s5 = sqrt(u5*u5+v5*v5)
s = s8*s5
z_mask = s == 0
nz_mask = s != 0
sdir = ndarray(tt.shape, float32)
sdir[z_mask] = 0
sdir[nz_mask] = (u5[nz_mask]*v8[nz_mask]-v5[nz_mask]*u8[nz_mask])/s[nz_mask]
tt49= ndarray(tt.shape, float32)
tt49_mask = tt>=49
tt49[tt<49] = 0
tt49[tt49_mask] = tt[tt49_mask]-49.0
result = 12*(td8-273.15)
result += 20*tt49
result += 2*1.944*s8
result += s5*1.944
result += 125*(sdir+0.2)
return result

View file

@ -1,51 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# Returns tempFromTenths if present, otherwise returns temperature
#
# ----------------------------------------------------------------
import numpy
from numpy import concatenate
from numpy import isnan
from unit import celciusToKelvin
# build an array that contains the higher resolution tempFromTenths where possible, otherwise use the lower resolution temperature
def combineTemps(temperature,tempFromTenths):
return numpy.where(isnan(tempFromTenths), temperature, tempFromTenths)
##
# Returns tempFromTenths if present, otherwise returns temperature
#
# @param temperature : temperature
# @param tempFromTenths : tempFromTenths
# @return:
# @rtype:
#
def execute1(temperature,tempFromTenths):
return celciusToKelvin(combineTemps(temperature,tempFromTenths))
def execute2(temperature, temp2):
temp2 = temp2.reshape(-1, 1);
result = concatenate((temp2, temperature), 1)
return result
def execute3(temperatureStation,tempFromTenthsStation):
return celciusToKelvin(combineTemps(temperatureStation, tempFromTenthsStation))

View file

@ -1,33 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
# ----------------------------------------------------------------
#
#
# ----------------------------------------------------------------
import numpy
from numpy import isnan
def execute(temperature, tempFromTenths, accum_temperature24, accum_tempFromTenths24):
T0 = numpy.where(isnan(tempFromTenths), temperature, tempFromTenths)
T0 = (T0 * 1.8) + 32
T24 = numpy.where(isnan(accum_tempFromTenths24), accum_temperature24, accum_tempFromTenths24)
T24 =(T24 * 1.8) + 32
return T0 - T24

View file

@ -1,28 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
from numpy import where
def execute(TP3hrtPlus0, TP3hrtMinus10800):
TP3hrtPlus0 = where((TP3hrtPlus0 > 4e13),0,TP3hrtPlus0);
TP3hrtMinus10800 = where((TP3hrtMinus10800 > 4e13),0,TP3hrtMinus10800);
result = TP3hrtPlus0+TP3hrtMinus10800
return result

View file

@ -1,37 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
# ----------------------------------------------------------------
# calculate virtual temperature from the pressure,
# temperature, and relative humidity.
# ----------------------------------------------------------------
from numpy import clip
from numpy import exp
def execute(P,T,RH):
"Calculate virtual temperature from temperature(K), relative \
humidity(0 to 100) and Pressure."
rhqc=clip(RH,1.0,100.0)
k = T
eee = exp(21.0827887-0.0091379024*k-6106.396/k)
tv = T * P / (P - RH * eee)
return tv

View file

@ -1,253 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
##
## @file TW.py
import gridslice
from numpy import abs
from numpy import bool
from numpy import clip
from numpy import empty
from numpy import exp
from numpy import log
from numpy import rank
from numpy import sqrt
from numpy import zeros
from numpy import where
## Constants declared in the original Fortran.
# I wish they had descriptive names.
## @var c1: constant used in multiple functions.
c1 = 0.0091379024
## @var c2: constant used in multiple functions.
c2 = 6106.3960
##
# Calculate wet bulb temperature in degrees C
# from pressure, temperature, and relative humidity.
#
# @param P: Pressure in millibars
# @param T: Temperature in degrees Kelvin
# @param RH: Relative humidity from 0 to 100
# @return: wet bulb temperature in degrees C as a numpy array
#
def execute1(P, T, RH):
# Find dry-bulb temp from temp and relative humidity
RH = clip(RH,1.0,100.0)
b = T * c1
b += c2/T
b -= log(RH/100.0)
val = b*b
val -= 223.1986
td = b - sqrt(val)
td /= 0.0182758048
TW = myTW(T, td, P)
return TW
def execute2(GH3D, T3D, RH3D, TILT):
P = zeros(T3D[0][0].shape, 'float32')
for i in range(len(RH3D[0])):
T = T3D[0][i]
RH = RH3D[0][i]
P.fill(RH3D[1][0][i])
TW = execute1(P, T, RH)
T3D[0][i] = TW
vc3d = GH3D[0]
slice3d = T3D[0]
slice = gridslice.createNumpySlice(vc3d, slice3d, TILT, 2)
slice = where(slice == 9.99999993e+36, - 999999.0, slice)
return slice
##
# This function takes arrays of values and filters them against expected values
# for terrestrial weather. Where the dewpoint temperature is greater than the
# observed temperature or the saturation vapor pressure is below 10 (millibars?),
# it returns a TW[i,j] estimate of (T[i,j]+Td[i,j])/2. Where T[i,j] is less than 100
# degrees Kelvin, it returns a TW[i,j] estimate of T[i,j]. Other cells are passed on
# to the calcTW() function for iterative estimation of TW.
#
# @param T: Observed temperature in degrees Kelvin.
# @param Td: Calculated dewpoint temperature in degrees Kelvin.
# @param P: Pressure in millibars
# @return: TW, the wet-bulb temperature
# @rtype: assumed to be a numpy array.
def myTW(T, Td, P):
td_lt_t = Td < T
t_ge_cold = T>=100.0
typical = td_lt_t & t_ge_cold
satVP = zeros(T.shape, dtype=T.dtype)
satVP[typical] = calcSatVP(T[typical])
norm_svp = satVP <= 10.0
sane = typical & norm_svp
veryCold = td_lt_t & ~t_ge_cold
td_satvp_bad = ~( td_lt_t & norm_svp ) # includes NaN cells
if rank(P)==0:
sel = lambda val, mask : val
else:
sel = lambda val, mask : val[mask]
TW = empty(T.shape, dtype=T.dtype)
TW[td_satvp_bad] = (T[td_satvp_bad] + Td[td_satvp_bad])/2
TW[veryCold] = T[veryCold]
TW[sane] = calcTW(T[sane], Td[sane], sel(P,sane), satVP[sane])
return(TW)
##
# Calculate saturation vapor pressure from T.
#
# @param T: Temperature in degrees K
# @return: ew, saturation vapor pressure (in millibars?)
#
def calcSatVP(T):
"Calculate saturation vapor pressure from T."
c0 = 26.66082
ew = -c1*T
ew += c0
ew -= c2/T
return ew
##
# Calculate TW iteratively from the input parameters. This is called
# from myTW(); the parameter arrays have already been tested for
# "reasonable" values.
#
# @param T: temperature in degrees K
# @param Td: dry-bulb temperature in degrees K
# @param P: pressure (in millibars?)
# @param satVP: saturation vapor pressure (in millibars?)
# @return: TW, the estimated wet-bulb temperature.
# @rtype: assumed to be a numpy array
#
def calcTW(T, Td, P, satVP):
""
f = 0.0006355
miniscule = 1e-5
satVPmin = -50.0
satVPmax = 10.0
# One more special case: "ridiculously small" dewpoint vapor pressure.
Tdx = Td
dewptVP = calcSatVP(Tdx)
dvpLtM50 = dewptVP < satVPmin
while( any(dvpLtM50) ):
Tdx[dvpLtM50] = Tdx[dvpLtM50] + 10
dewptVP[dvpLtM50] = calcSatVP(Tdx[dvpLtM50])
# Update dvpLtM50.
# The fancy indexing is to skip cells that are already False,
# since they won't become True as a result of the loop.
dvpLtM50[dvpLtM50] = dewptVP[dvpLtM50] < satVPmin
fp = P * f
if rank(fp) == 0:
sel = lambda val, mask : val
else:
sel = lambda val, mask : val[mask]
# we need e**satVP and e**dewptVP for calculating s and the initial TW estimate
# (e is the base of the natural logarithm).
expSatVP = exp(satVP)
expDewptVP = exp(dewptVP)
s = expSatVP-expDewptVP
s /= T-Tdx
TW = T*fp
TW += Tdx*s
TW /= fp+s
# new, non-exp'd satVP for the loop below
satVP = calcSatVP(TW)
# create some boolean arrays to use in the for loop.
# creating them once here allows us to use the out parameter of
# the numpy boolean functions to avoid repeatedly allocating space.
# satVPinRange = empty(shape(TW), dtype=bool)
# svrTemp = empty(shape(TW), dtype=bool)
# At each step of the iteration, esat(Tw)-esat(Td) is compared to
# (T-Tw)*p/(eps*L). When that difference is less than one part in
# 10000 of esat(Tw), or ten iterations have been done, the iteration stops.
# This is basically trying to find the value of Tw where de is 0. The
# value s is the derivative of de with respect to Tw, a fairly standard
# numerical technique for finding the zero value of a function.
for iteration in xrange(10):
satVPinRange = (satVPmin<=satVP) & (satVP<=satVPmax)
if not any(satVPinRange):
break
expSatVP = exp(satVP[satVPinRange])
de = calcDe(T[satVPinRange], TW[satVPinRange], expSatVP,
expDewptVP[satVPinRange], sel(fp,satVPinRange))
mag = abs(de/expSatVP)
significant = miniscule <= mag
if any(significant):
satVPinRange[satVPinRange] = significant
# now we can mix TW and de to obtain a new estimate for TW:
TW[satVPinRange] = reestimateTW( TW[satVPinRange],
expSatVP[significant],
de[significant],
sel(fp, satVPinRange) )
# Get a new satVP to see if the iteration can stop.
# Only cells with a different TW will have a different satVP.
satVP[satVPinRange] = calcSatVP(TW[satVPinRange])
else:
break
return TW
##
# Calculate the value of the de variable used in calcTW.
#
# @param T: The temperature in degrees Kelvin.
# @param TW: The current estimate of the wet-bulb temperature in degrees Kelvin.
# @param expSatVP: The current saturation vapor pressure estimate
# @param expDewptVP: The current dewpoint vapor pressure estimate
# @param fp: The pressure times the f constant.
# @return: the calculated value for de
#
def calcDe(T, TW, expSatVP, expDewptVP, fp):
"Calculate de from temp, estimated "
de = T-TW
de *= fp
de += expDewptVP
de -= expSatVP
return de
def reestimateTW(TW, expSatVP, de, fp):
"Find a new estimate for TW based on the current estimate, \
the saturation vapor pressure, de, and fp."
s = -c2/(TW*TW)
s += c1
s *= expSatVP
s -= fp
newTW = TW - de/s
return newTW

View file

@ -1,163 +0,0 @@
##
# This software was developed and / or modified by Raytheon Company,
# pursuant to Contract DG133W-05-CQ-1067 with the US Government.
#
# U.S. EXPORT CONTROLLED TECHNICAL DATA
# This software product contains export-restricted data whose
# export/transfer/disclosure is restricted by U.S. law. Dissemination
# to non-U.S. persons whether in the United States or abroad requires
# an export license or other authorization.
#
# Contractor Name: Raytheon Company
# Contractor Address: 6825 Pine Street, Suite 340
# Mail Stop B8
# Omaha, NE 68106
# 402.291.0100
#
# See the AWIPS II Master Rights File ("Master Rights File.pdf") for
# further licensing information.
###
import numpy as np
import AdiabaticTemp
# create an array of base pressures
pBases = np.array([200,350,500,600,700,850,1000], dtype=np.float32)
# create an array of splitpoints for pressure bins
pSplits = np.array([250,400,550,650,750,900], dtype=np.float32)
tmax = 333
tmin = 193
# create a lookup table of adiabatic temperatures for chosen base pressures
parr, tarr = np.ogrid[0:len(pBases), tmin:tmax]
parr = pBases[parr]
tarr = np.array(tarr, dtype=np.float32)
lookupTable = AdiabaticTemp.execute(tarr, parr)
##
# Calculate the saturation temperature of an equivalent temperature at
# the specified pressure. Another way of expressing it is to find the
# temperature TofTe such that AdiabaticTemp.execute(TofTe,P)=T.
#
# @param T: Adiabatic temperature of TofTe at P (K)
# @type T: Numpy array of float32
# @param P: The pressure at which the calculation is made (mb)
# @type P: Numpy array of float32
# @return: TofTe
# @rtype: Numpy array
def execute(T,P):
"Find Tans such that AdiabaticTemp.execute(Tans,P) == T."
if np.shape(T) != np.shape(P):
bcastOnes = np.ones(T.shape, dtype=np.float32) * np.ones(np.shape(P), dtype=np.float32)
T = T * bcastOnes
P = P * bcastOnes
pBaseIdx = np.searchsorted(pSplits, P)
# find cells that are too low or too high
lowTMask = T<lookupTable[pBaseIdx,1]
goodhighTMask = ((T<=lookupTable[pBaseIdx,-1]) & ~(T<=tmax)) # NaNs in T are True in highTMask
Tgood = np.array(T, dtype=np.float32)
Tgood[goodhighTMask] = tmax
highTMask = ~(Tgood<=lookupTable[pBaseIdx,-1])
goodMask = ~( lowTMask | highTMask)
Tgood = Tgood[goodMask]
pBaseIdx = pBaseIdx[goodMask]
# find t1 and t2, lower/upper integer estimate arrays
t1 = np.ones(Tgood.shape, dtype=np.integer) * tmin
t2 = np.array(Tgood, dtype=np.integer)
TMask = t2-t1>=3
t1Mask = np.ones(np.shape(TMask), dtype=np.bool)
t2Mask = np.ones(np.shape(TMask), dtype=np.bool)
# Tgood cells just above the lower threshold can result in False TMask cells.
# Don't allow t1Mask or t2Mask to be True where TMask is False.
t1Mask[:] = TMask
t2Mask[:] = TMask
while(np.any(TMask)):
# t = (int)((t1+t2)/2), for cells where TMask is True
tEst = np.array((t1[TMask]+t2[TMask])/2, dtype=np.integer)
telu = lookupTable[pBaseIdx[TMask], tEst-tmin]
# where adiabaticTemp(t,pBase) is > T, adjust t2 down.
cmpMask = telu > Tgood[TMask]
t2Mask[TMask] = cmpMask
t2[t2Mask] = tEst[cmpMask]
# where adiabaticTemp(t,pBase) is < T, adjust t1 up.
cmpMask = telu < Tgood[TMask]
t1Mask[TMask] = cmpMask
t1[t1Mask] = tEst[cmpMask]
# where we hit it exactly, set t1 and t2 to one-off.
cmpMask = (telu == Tgood[TMask])
t1Mask[TMask] = cmpMask
t2Mask[TMask] = cmpMask
t1[t1Mask] = tEst[cmpMask]-1
t2[t2Mask] = tEst[cmpMask]+1
# Clean up t1Mask and t2Mask for next iteration.
# Without this, t?[t?Mask] could expect too many values.
t1Mask[TMask] = False
t2Mask[TMask] = False
# figure out which cells still need the estimate narrowed.
TMask[TMask] = t2[TMask]-t1[TMask]>=3
# weight t1 and t2 by the ratio of base pressure to P
# (t1 and t2 become floating-pt here)
pBase = pBases[pBaseIdx]
weight = np.sqrt(pBase/P[goodMask])
t1 = (1-weight) * lookupTable[pBaseIdx, t1-tmin] + weight * t1.astype(T.dtype)
t2 = (1-weight) * lookupTable[pBaseIdx, t2-tmin] + weight * t2.astype(T.dtype)
# now find a new weight based on the difference between T and
# the adiabatic temperature of t1 and t2.
if np.isscalar(P):
PP = P
else:
PP = P[goodMask]
diff1 = T[goodMask] - AdiabaticTemp.execute(t1,PP)
diff2 = AdiabaticTemp.execute(t2,PP) - T[goodMask]
weight = diff2/(diff1+diff2)
Tans = weight*t1 + (1-weight)*t2
diff = AdiabaticTemp.execute(Tans,PP) - T[goodMask]
# Iterate to find result to nearest .01 degree K
TMask[:] = True # all cells still to calculate
for loopcount in xrange(10):
if np.any(TMask):
# see which estimates are high and low
t1Mask[TMask] = diff[TMask] < -0.01
t2Mask[TMask] = diff[TMask] > 0.01
# adjust down the high estimates
t2[t2Mask] = Tans[t2Mask]
diff2[t2Mask] = diff[t2Mask]
# adjust up the low estimates
t1[t1Mask] = Tans[t1Mask]
diff1[t1Mask] = -diff[t1Mask]
# calculate a new weight, estimated T, and diff for active cells
if np.isscalar(P):
PPP = P
else:
PPP = PP[TMask]
weight[TMask] = diff2[TMask]/(diff1[TMask]+diff2[TMask])
Tans[TMask] = weight[TMask]*t1[TMask] + (1-weight[TMask])*t2[TMask]
diff[TMask] = AdiabaticTemp.execute(Tans[TMask], PPP) - T[TMask]
# eliminate cells we've finished
TMask[TMask] = t1Mask[TMask] | t2Mask[TMask]
else:
break
# Create an empty fully masked array for the result
result = np.empty(np.shape(T), dtype=np.float32)
result[:] = np.NaN
result[goodMask] = Tans
# replace the very low input values with the original T
result[lowTMask] = T[lowTMask]
return result

Some files were not shown because too many files have changed in this diff Show more