archive-gh-me/MeteoInfo
1
0
Fork 0
mirror of https://github.com/meteoinfo/MeteoInfo.git synced 2025-12-08 20:36:05 +00:00
Code Issues Projects Releases Wiki Activity

add velocity_potential and stream_function functions

Browse Source
This commit is contained in:
wyq 2024-06-09 16:57:35 +08:00
parent 8acd7e0837
commit b11475d83f
12 changed files with 223 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
meteoinfo-common/src/main/java/org/meteoinfo/common/util/GlobalUtil.java
View File

@ -67,7 +67,7 @@ import java.util.zip.ZipInputStream;
public static String getVersion(){
String version = GlobalUtil.class.getPackage().getImplementationVersion();
if (version == null || version.equals("")) {
version = "3.8.10";
version = "3.8.11";
}
return version;
}

19
meteoinfo-dataframe/src/main/java/org/meteoinfo/dataframe/Column.java
View File

@ -13,6 +13,8 @@ import org.meteoinfo.ndarray.util.DataTypeUtil;
import java.text.DecimalFormat;
import java.time.format.DateTimeFormatter;
import java.util.Formatter;
import java.util.Locale;
/**
*
@ -149,10 +151,17 @@ public class Column {
*/
public static Column factory(String name, Array array){
DataType dtype = array.getDataType();
if (dtype == DataType.OBJECT && (array.getObject(0) instanceof LocalDateTime)){
DateTimeColumn col = new DateTimeColumn(name);
col.updateFormat(array);
return col;
switch (dtype) {
case DATE:
DateTimeColumn col = new DateTimeColumn(name);
col.updateFormat(array);
return col;
case OBJECT:
if (array.getObject(0) instanceof LocalDateTime) {
col = new DateTimeColumn(name);
col.updateFormat(array);
return col;
}
}
return new Column(name, dtype);
}
@ -261,7 +270,7 @@ public class Column {
}
return ((LocalDateTime) o).format(this.dateTimeFormatter);
} else
return String.format(format, o);
return String.format(Locale.US, format, o);
}
}

2
meteoinfo-dataframe/src/main/java/org/meteoinfo/dataframe/DataFrame.java
View File

@ -3166,7 +3166,7 @@ public class DataFrame implements Iterable {
if (formats.get(i) == null) {
vstr = this.getValue(j, i).toString();
} else {
vstr = String.format(formats.get(i), this.getValue(j, i));
vstr = String.format(Locale.US, formats.get(i), this.getValue(j, i));
}
if (line.isEmpty()) {
line = vstr;

36
meteoinfo-lab/milconfig.xml
View File

@ -1,32 +1,34 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<MeteoInfo File="milconfig.xml" Type="configurefile">
<Path OpenPath="D:\Working\MIScript\Jython\mis\io\radar\cinrad">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\dataframe\series"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\dataframe"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\contour"/>
<Path OpenPath="D:\Working\MIScript\Jython\mis\meteo\calc">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\image"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\calc"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\imshow"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\gridshow"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\plot"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\radar"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\radar\cinrad"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\dataframe"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\wind"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\arrow"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\matlab"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\calc"/>
</Path>
<File>
<OpenedFiles>
<OpenedFile File="D:\Working\MIScript\Jython\mis\io\radar\cinrad\radar_pa_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\io\radar\cinrad\radar_cc_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\io\radar\cinrad\radar_ca_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_3.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_2.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\velocity_potential_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\calc\velocity_potential_2.py"/>
</OpenedFiles>
<RecentFiles>
<RecentFile File="D:\Working\MIScript\Jython\mis\io\radar\cinrad\radar_pa_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\io\radar\cinrad\radar_cc_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\io\radar\cinrad\radar_ca_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_3.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\stream_function_2.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\velocity_potential_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\meteo\calc\velocity_potential_2.py"/>
</RecentFiles>
</File>
<Font>
@ -34,5 +36,5 @@
</Font>
<LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
<Figure DoubleBuffering="true"/>
<Startup MainFormLocation="-7,0" MainFormSize="1340,820"/>
<Startup MainFormLocation="-7,0" MainFormSize="1398,806"/>
</MeteoInfo>

BIN
meteoinfo-lab/pylib/mipylib/meteolib/calc/kinematics$py.class
View File

Binary file not shown.

126
meteoinfo-lab/pylib/mipylib/meteolib/calc/kinematics.py
View File

@ -14,11 +14,133 @@ from .. import constants
__all__ = [
'cdiff', 'divergence', 'vorticity', 'advection', 'absolute_vorticity', 'potential_vorticity_barotropic',
'ageostrophic_wind', 'frontogenesis', 'geostrophic_wind', 'montgomery_streamfunction',
'potential_vorticity_baroclinic', 'shearing_deformation', 'storm_relative_helicity',
'stretching_deformation', 'total_deformation', 'inertial_advective_wind', 'q_vector'
'potential_vorticity_baroclinic', 'shearing_deformation', 'stream_function', 'storm_relative_helicity',
'stretching_deformation', 'total_deformation', 'inertial_advective_wind', 'q_vector', 'velocity_potential'
]
def _dx_atmos(longitude, latitude):
dx = np.gradient(longitude, axis=1)
dx = dx * (np.pi / 180) * constants.Re * np.cos(latitude * np.pi / 180)
return dx
def _dy_atmos(latitude):
dy = np.gradient(latitude, axis=0)
dy = dy * (np.pi / 180) * constants.Re
return dy
def _divh_atmos(longitude, latitude, u, v):
dx = _dx_atmos(longitude, latitude)
dy = _dy_atmos(latitude)
du = np.gradient(u, axis=1)
dv = np.gradient(v, axis=0)
divh = du / dx + dv / dy - v * np.tan(latitude * np.pi / 180) / constants.Re
divh[np.abs(latitude)==90] = np.nan
return divh
def _curlz_atmos(longitude, latitude, u, v):
dx = _dx_atmos(longitude, latitude)
dy = _dy_atmos(latitude)
du = np.gradient(u, axis=0)
dv = np.gradient(v, axis=1)
curlz = dv / dx - du / dy + u * np.tan(latitude * np.pi / 180) / constants.Re
curlz[np.abs(latitude)==90] = np.nan
return curlz
def _grad_atmos(longitude, latitude, H):
dx = _dx_atmos(longitude, latitude)
dy = _dy_atmos(latitude)
grad_y, grad_x = np.gradient(H)
grad_x = grad_x / dx
grad_y = grad_y / dy
return grad_x, grad_y
def velocity_potential(longitude, latitude, u, v, loop_max=1e10, epsilon=1e-10):
"""
Calculate velocity potential using Richardson iterative method.
Parameters
----------
lon : (M, N) `array`
longitude array
lat : (M, N) `array`
latitude array
u : (M, N) `array`
x component of the wind
v : (M, N) `array`
y component of the wind
Returns
-------
A 3-item tuple of arrays
velocity potential and divergence wind component
"""
if isinstance(loop_max, float):
loop_max = int(loop_max)
divh = _divh_atmos(longitude, latitude, u, v) # vertical divergence
dx2 = _dx_atmos(longitude, latitude)**2 # square of latitude gradient
dy2 = _dy_atmos(latitude)**2 # square of longitude gradient
phi = np.NDArray(MeteoMath.richardsonIteration(divh._array, dx2._array, dy2._array,
loop_max, epsilon))
phi = -phi
# divergence wind
DphiDx, DphiDy = _grad_atmos(longitude, latitude, phi)
Uphi = DphiDx
Vphi = DphiDy
# make boundary NaN
phi[0,:] = np.nan; Uphi[0,:] = np.nan; Vphi[0,:] = np.nan;
phi[-1,:] = np.nan; Uphi[-1,:] = np.nan; Vphi[-1,:] = np.nan;
phi[:,0] = np.nan; Uphi[:,0] = np.nan; Vphi[:,0] = np.nan;
phi[:,-1] = np.nan; Uphi[:,-1] = np.nan; Vphi[:,-1] = np.nan;
return phi, Uphi, Vphi
def stream_function(longitude, latitude, u, v, loop_max=int(1e10), epsilon=1e-10):
"""
Calculate stream function using Richardson iterative method.
Parameters
----------
lon : (M, N) `array`
longitude array
lat : (M, N) `array`
latitude array
u : (M, N) `array`
x component of the wind
v : (M, N) `array`
y component of the wind
Returns
-------
A 3-item tuple of arrays
stream function and vorticity wind component
"""
if isinstance(loop_max, float):
loop_max = int(loop_max)
curlz = _curlz_atmos(longitude, latitude, u, v) # vorticity
dx2 = _dx_atmos(longitude, latitude)**2 # square of latitude gradient
dy2 = _dy_atmos(latitude)**2 # square of longitude gradient
psi = np.NDArray(MeteoMath.richardsonIteration(curlz._array, dx2._array, dy2._array,
loop_max, epsilon))
psi = -psi
# vorticity wind
DpsiDx, DpsiDy = _grad_atmos(longitude, latitude, psi)
Upsi = -DpsiDy
Vpsi = DpsiDx
# make boundary NaN
psi[0,:] = np.nan; Upsi[0,:] = np.nan; Vpsi[0,:] = np.nan;
psi[-1,:] = np.nan; Upsi[-1,:] = np.nan; Vpsi[-1,:] = np.nan;
psi[:,0] = np.nan; Upsi[:,0] = np.nan; Vpsi[:,0] = np.nan;
psi[:,-1] = np.nan; Upsi[:,-1] = np.nan; Vpsi[:,-1] = np.nan;
return psi, Upsi, Vpsi
def cdiff(a, dimidx):
"""
Performs a centered difference operation on an array in a specific direction

3
meteoinfo-lab/pylib/mipylib/meteolib/calc/tools.py
View File

@ -438,6 +438,7 @@ def first_derivative(f, axis=None, x=None, delta=None):
This uses 3 points to calculate the derivative, using forward or backward at the edges of
the grid as appropriate, and centered elsewhere. The irregular spacing is handled
explicitly, using the formulation as specified by [Bowen2005]_.
Parameters
----------
f : array-like
@ -454,10 +455,12 @@ def first_derivative(f, axis=None, x=None, delta=None):
delta : array-like, optional
Spacing between the grid points in `f`. Should be one item less than the size
of `f` along `axis`.
Returns
-------
array-like
The first derivative calculated along the selected axis`
See Also
--------
second_derivative

1
meteoinfo-lab/pylib/mipylib/meteolib/constants/default.py
View File

@ -8,6 +8,7 @@ Cp_d = 1005 #Specific heat at constant pressure for dry air (J kg^-1)
epsilon = Mw / Md
kappa = 0.286
degCtoK = 273.15 # Temperature offset between K and C (deg C)
Re = earth_avg_radius = 6371008.7714 # Earth average radius with meters
g = 9.8 # Gravitational acceleration (m / s^2)
sat_pressure_0c = 6.112 #Saturation pressure at 0 degree (hPa)
T_BASE = 300.

BIN
meteoinfo-lab/pylib/mipylib/meteolib/meteo$py.class
View File

Binary file not shown.

41
meteoinfo-math/src/main/java/org/meteoinfo/math/meteo/MeteoMath.java
View File

@ -11,6 +11,7 @@ import org.meteoinfo.ndarray.Array;
import org.meteoinfo.ndarray.DataType;
import org.meteoinfo.ndarray.Index;
import org.meteoinfo.ndarray.IndexIterator;
import org.meteoinfo.ndarray.math.ArrayUtil;
import java.util.Arrays;
@ -1130,4 +1131,44 @@ public class MeteoMath {
return r;
}
/**
* Richardson's iteration method for calculating Poisson equation
*
* @param array Input array - 2D
* @param dx2 X gradient square
* @param dy2 Y gradient square
* @param loopMax Maximum loop number
* @param epsilon Minimum error value
* @return Poisson equation output
*/
public static Array richardsonIteration(Array array, Array dx2, Array dy2, long loopMax, double epsilon) {
int[] shape = array.getShape();
int m = shape[0];
int n = shape[1];
Array phi = ArrayUtil.zeros(shape, DataType.DOUBLE);
Array res = ArrayUtil.full(shape, -9999., DataType.DOUBLE);
for (int k=0; k < loopMax; k++) {
for (int i = 1; i < m - 1; i++) {
for (int j = 1; j < n - 1; j++) {
// calculate residual
res.setDouble(i * n + j, (phi.getDouble((i + 1) * n + j) +
phi.getDouble((i - 1) * n + j) - 2 * phi.getDouble(i * n + j)) /
dx2.getDouble(i * n + j) + (phi.getDouble(i * n + j + 1) +
phi.getDouble(i * n + j - 1) - 2 * phi.getDouble(i * n + j)) /
dy2.getDouble(i * n + j) + array.getDouble(i * n + j));
// iterator
phi.setDouble(i * n + j, phi.getDouble(i * n + j) + res.getDouble(i * n + j) /
(2 / dx2.getDouble(i * n + j) + 2 / dy2.getDouble(i * n + j)));
}
}
if (ArrayMath.max(res).doubleValue() < epsilon) {
break;
}
}
return phi;
}
}

24
meteoinfo-ndarray/src/main/java/org/meteoinfo/ndarray/math/ArrayUtil.java
View File

@ -1062,15 +1062,11 @@ public class ArrayUtil {
* @param dtype Data type
* @return Array Result array
*/
public static Array full(List<Integer> shape, Object fillValue, DataType dtype) {
int[] ashape = new int[shape.size()];
for (int i = 0; i < shape.size(); i++) {
ashape[i] = shape.get(i);
}
public static Array full(int[] shape, Object fillValue, DataType dtype) {
if (dtype == null) {
dtype = ArrayMath.getDataType(fillValue);
}
Array a = Array.factory(dtype, ashape);
Array a = Array.factory(dtype, shape);
for (int i = 0; i < a.getSize(); i++) {
a.setObject(i, fillValue);
@ -1079,6 +1075,22 @@ public class ArrayUtil {
return a;
}
/**
* Return a new array of given shape and type, filled with fill value.
*
* @param shape Shape
* @param fillValue Fill value
* @param dtype Data type
* @return Array Result array
*/
public static Array full(List<Integer> shape, Object fillValue, DataType dtype) {
int[] ashape = new int[shape.size()];
for (int i = 0; i < shape.size(); i++) {
ashape[i] = shape.get(i);
}
return full(ashape, fillValue, dtype);
}
/**
* Return a new array of given shape and type, filled with fill value.
*

2
pom.xml
View File

@ -35,7 +35,7 @@
<properties>
<project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
<java.version>1.8</java.version>
<revision>3.8.10</revision>
<revision>3.8.11</revision>
<maven.compiler.source>8</maven.compiler.source>
<maven.compiler.target>8</maven.compiler.target>
<maven.compiler.release>8</maven.compiler.release>