add hcurl and hdivg functions in meteo module

This commit is contained in:
wyq 2020-11-24 09:37:05 +08:00
parent 4dbf66e9ec
commit adfa6ea5a7
4 changed files with 145 additions and 7 deletions

View File

@ -1,7 +1,6 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<MeteoInfo File="milconfig.xml" Type="configurefile">
<Path OpenPath="D:\Working\MIScript\Jython\mis\io\burf">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\weather"/>
<Path OpenPath="D:\Working\MIScript\Jython\mis\meteo">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\polar"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\topology"/>
@ -11,22 +10,23 @@
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\dataframe"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\toolbox"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\toolbox\miml"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\toolbox\miml\cluster"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\burf"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
</Path>
<File>
<OpenedFiles>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\polar\polar_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\plot3_multi_color_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\io\burf\bufr_6.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\vort_digv.py"/>
</OpenedFiles>
<RecentFiles>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\polar\polar_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\plot3_multi_color_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\io\burf\bufr_6.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\meteo\vort_digv.py"/>
</RecentFiles>
</File>
<Font>

View File

@ -5,7 +5,7 @@
# Note: Jython, some functions code revised from MetPy
#-----------------------------------------------------
from org.meteoinfo.math import ArrayUtil
from org.meteoinfo.math import ArrayMath
from org.meteoinfo.math.meteo import MeteoMath
import mipylib.numeric as np
@ -15,7 +15,7 @@ import constants as constants
__all__ = [
'cumsimp','dewpoint','dewpoint2rh','dewpoint_rh','dry_lapse','ds2uv','equivalent_potential_temperature',
'exner_function','flowfun','h2p',
'exner_function','flowfun','h2p','hcurl','hdivg',
'moist_lapse','p2h','potential_temperature','qair2rh','rh2dewpoint',
'sigma_to_pressure','tc2tf',
'temperature_from_potential_temperature','tf2tc','uv2ds','pressure_to_height_std',
@ -63,6 +63,74 @@ def ds2uv(d, s):
else:
r = MeteoMath.ds2uv(d, s)
return r[0], r[1]
def hcurl(u, v, x=None, y=None):
'''
Calculates the vertical component of the curl (ie, vorticity). The data should be lon/lat projection.
:param u: (*array*) U component 2D array.
:param v: (*array*) V component 2D array.
:param x: (*array*) X coordinate.
:param y: (*array*) Y coordinate.
:returns: Array of the vertical component of the curl.
'''
ny, nx = u.shape
if x is None:
if isinstance(u, DimArray):
x = u.dimvalue(1)
else:
x = np.arange(nx)
elif isinstance(x, (list, tuple)):
x = np.array(x)
if y is None:
if isinstance(v, DimArray):
y = u.dimvalue(0)
else:
y = np.arange(ny)
elif isinstance(y, (list, tuple)):
y = np.array(y)
r = ArrayMath.hcurl(u.asarray(), v.asarray(), x.asarray(), y.asarray())
if isinstance(u, DimArray):
return DimArray(NDArray(r), u.dims, u.fill_value, u.proj)
else:
return NDArray(r)
def hdivg(u, v, x=None, y=None):
'''
Calculates the horizontal divergence using finite differencing. The data should be lon/lat projection.
:param u: (*array*) U component array.
:param v: (*array*) V component array.
:param x: (*array*) X coordinate.
:param y: (*array*) Y coordinate.
:returns: Array of the horizontal divergence.
'''
ny, nx = u.shape
if x is None:
if isinstance(u, DimArray):
x = u.dimvalue(1)
else:
x = np.arange(nx)
elif isinstance(x, (list, tuple)):
x = np.array(x)
if y is None:
if isinstance(v, DimArray):
y = u.dimvalue(0)
else:
y = np.arange(ny)
elif isinstance(y, (list, tuple)):
y = np.array(y)
r = ArrayMath.hdivg(u.asarray(), v.asarray(), x.asarray(), y.asarray())
if isinstance(u, DimArray):
return DimArray(NDArray(r), u.dims, u.fill_value, u.proj)
else:
return NDArray(r)
def p2h(press):
"""

View File

@ -7356,6 +7356,41 @@ public class ArrayMath {
return gData;
}
/**
* Calculates the vertical component of the curl (ie, vorticity)
*
* @param uData U component
* @param vData V component
* @param xx X dimension value
* @param yy Y dimension value
* @return Curl
*/
public static Array hcurl(Array uData, Array vData, Array xx, Array yy) {
xx = xx.copyIfView();
yy = yy.copyIfView();
int rank = uData.getRank();
int[] shape = uData.getShape();
Array lonData = Array.factory(DataType.DOUBLE, shape);
Array latData = Array.factory(DataType.DOUBLE, shape);
Index index = lonData.getIndex();
int[] current;
for (int i = 0; i < lonData.getSize(); i++) {
current = index.getCurrentCounter();
lonData.setDouble(index, xx.getDouble(current[rank - 1]));
latData.setDouble(index, yy.getDouble(current[rank - 2]));
index.incr();
}
Array dv = cdiff(vData, rank - 1);
Array dx = mul(cdiff(lonData, rank - 1), Math.PI / 180);
Array du = cdiff(mul(uData, cos(mul(latData, Math.PI / 180))), rank - 2);
Array dy = mul(cdiff(latData, rank - 2), Math.PI / 180);
Array gData = div(sub(div(dv, dx), div(du, dy)), mul(cos(mul(latData, Math.PI / 180)), 6.37e6));
return gData;
}
/**
* Calculates the horizontal divergence using finite differencing
*
@ -7388,6 +7423,41 @@ public class ArrayMath {
return gData;
}
/**
* Calculates the horizontal divergence using finite differencing
*
* @param uData U component
* @param vData V component
* @param xx X dimension value
* @param yy Y dimension value
* @return Divergence
*/
public static Array hdivg(Array uData, Array vData, Array xx, Array yy) {
xx = xx.copyIfView();
yy = yy.copyIfView();
int rank = uData.getRank();
int[] shape = uData.getShape();
Array lonData = Array.factory(DataType.DOUBLE, shape);
Array latData = Array.factory(DataType.DOUBLE, shape);
Index index = lonData.getIndex();
int[] current;
for (int i = 0; i < lonData.getSize(); i++) {
current = index.getCurrentCounter();
lonData.setDouble(index, xx.getDouble(current[rank - 1]));
latData.setDouble(index, yy.getDouble(current[rank - 2]));
index.incr();
}
Array du = cdiff(uData, rank - 1);
Array dx = mul(cdiff(lonData, rank - 1), Math.PI / 180);
Array dv = cdiff(mul(vData, cos(mul(latData, Math.PI / 180))), rank - 2);
Array dy = mul(cdiff(latData, rank - 2), Math.PI / 180);
Array gData = div(add(div(du, dx), div(dv, dy)), mul(cos(mul(latData, Math.PI / 180)), 6.37e6));
return gData;
}
/**
* Take magnitude value from U/V grid data
*