add hypot, s_, index_exp functions

2025-12-08 20:36:05 +00:00 · 2023-07-31 14:33:52 +08:00 · 2023-07-31 14:33:52 +08:00 · c16a807b40
commit c16a807b40
parent dffb2da247
13 changed files with 176 additions and 28 deletions
--- a/.idea/misc.xml
+++ b/.idea/misc.xml
@ -1,4 +1,3 @@
-<?xml version="1.0" encoding="UTF-8"?>
 <project version="4">
  <component name="ExternalStorageConfigurationManager" enabled="true" />
  <component name="MavenProjectsManager">
--- a/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/GraphicFactory.java
+++ b/meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/GraphicFactory.java
@ -1312,7 +1312,7 @@ public class GraphicFactory {
        double c;
        ColorBreak cb;
        boolean fixZ = false;
-        double z = 0;
+        double x, y, z = 0;
        if (zdata.getSize() == 1 && xdata.getSize() > 1) {
            fixZ = true;
            z = zdata.getDouble(0);
@ -1322,13 +1322,18 @@ public class GraphicFactory {
        IndexIterator zIter = zdata.getIndexIterator();
        IndexIterator cIter = cdata.getIndexIterator();
        while (xIter.hasNext()) {
-            ps = new PointZShape();
-            if (fixZ) {
-                ps.setPoint(new PointZ(xIter.getDoubleNext(), yIter.getDoubleNext(), z));
-            } else {
-                ps.setPoint(new PointZ(xIter.getDoubleNext(), yIter.getDoubleNext(), zIter.getDoubleNext()));
-            }
+            x = xIter.getDoubleNext();
+            y = yIter.getDoubleNext();
            c = cIter.getDoubleNext();
+            if (!fixZ) {
+                z = zIter.getDoubleNext();
+            }
+            if (Double.isNaN(x) || Double.isNaN(y) || Double.isNaN(z) || Double.isNaN(c)) {
+                continue;
+            }
+
+            ps = new PointZShape();
+            ps.setPoint(new PointZ(x, y, z));
            cb = ls.findLegendBreak(c);
            graphics.add(new Graphic(ps, cb));
        }
--- a/meteoinfo-common/src/main/java/org/meteoinfo/common/util/GlobalUtil.java
+++ b/meteoinfo-common/src/main/java/org/meteoinfo/common/util/GlobalUtil.java
@ -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.6.3";
+             version = "3.6.4";
         }
         return version;
     }
--- a/meteoinfo-data/src/main/java/org/meteoinfo/data/meteodata/DataInfo.java
+++ b/meteoinfo-data/src/main/java/org/meteoinfo/data/meteodata/DataInfo.java
@ -723,5 +723,21 @@ import org.meteoinfo.projection.ProjectionInfo;
         return null;
     }

+     /**
+      * Get the data is Radial (Radar) or not
+      * @return Is Radial or not
+      */
+     public boolean isRadial() {
+         Attribute ra = findGlobalAttribute("featureType");
+         if (ra != null) {
+             String va = ra.getStringValue();
+             if (va.equalsIgnoreCase("RADIAL")) {
+                 return true;
+             }
+         }
+
+         return false;
+     }
+
     // </editor-fold>
 }
--- a/meteoinfo-data/src/main/java/org/meteoinfo/data/meteodata/radar/Transform.java
+++ b/meteoinfo-data/src/main/java/org/meteoinfo/data/meteodata/radar/Transform.java
@ -1,5 +1,7 @@
 package org.meteoinfo.data.meteodata.radar;

+import org.meteoinfo.ndarray.Array;
+import org.meteoinfo.ndarray.DataType;
 import org.meteoinfo.ndarray.math.ArrayMath;

 public class Transform {
@ -41,7 +43,7 @@ public class Transform {
        return new double[]{x, y, z};
    }

-    private static double azimuth(float x, float y) {
+    private static double xyToAzimuth(float x, float y) {
        double az = Math.PI / 2 - Math.atan2(x + y, 1);
        if (az < 0) {
            az = 2 * Math.PI + az;
@ -62,8 +64,39 @@ public class Transform {
                Math.cos(Math.sqrt(x * x + y * y)  / R));
        double elevation = Math.acos((R + z) * Math.sin(Math.sqrt(x * x + y * y) / R) / ranges) *
                180. / Math.PI;
-        double azimuth = azimuth(x, y);
+        double azimuth = xyToAzimuth(x, y);

        return new double[]{azimuth, ranges, elevation};
    }
+
+    /**
+     * Convert cartesian coordinate to antenna coordinate
+     * @param xa x coordinate array in meters
+     * @param ya y coordinate array in meters
+     * @param za z coordinate array in meters
+     * @param h Altitude of the instrument, above sea level, units:m
+     * @return Antenna coordinate from the radar
+     */
+    public static Array[] cartesianToAntenna(Array xa, Array ya, Array za, float h) {
+        xa = xa.copyIfView();
+        ya = ya.copyIfView();
+        za = za.copyIfView();
+
+        Array ranges = Array.factory(DataType.DOUBLE, xa.getShape());
+        Array azimuth = Array.factory(DataType.DOUBLE, xa.getShape());
+        Array elevation = Array.factory(DataType.DOUBLE, xa.getShape());
+        float x, y, z;
+
+        for (int i = 0; i < xa.getSize(); i++) {
+            x = xa.getFloat(i);
+            y = ya.getFloat(i);
+            z = za.getFloat(i);
+            double[] rr = cartesianToAntenna(x, y, z, h);
+            azimuth.setDouble(i, rr[0]);
+            ranges.setDouble(i, rr[1]);
+            elevation.setDouble(i, rr[2]);
+        }
+
+        return new Array[]{azimuth, ranges, elevation};
+    }
 }
--- a/meteoinfo-lab/milconfig.xml
+++ b/meteoinfo-lab/milconfig.xml
@ -1,36 +1,36 @@
 <?xml version="1.0" encoding="UTF-8" standalone="no"?>
 <MeteoInfo File="milconfig.xml" Type="configurefile">
-  <Path OpenPath="D:\Working\MIScript\Jython\mis\common_math\linalg">
-    <RecentFolder Folder="D:\Working\MIScript\Jython"/>
-    <RecentFolder Folder="D:\Working\MIScript"/>
-    <RecentFolder Folder="D:\Working\MIScript\cuace_dust"/>
-    <RecentFolder Folder="D:\Working\MIScript\cuace_dust\py"/>
-    <RecentFolder Folder="D:\Working\MIScript\cuace_dust\py\plot"/>
+  <Path OpenPath="D:\Working\MIScript\Jython\mis\satellite\gpm">
    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\array\complex"/>
    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\traj"/>
    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\map"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\radar"/>
    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\webmap"/>
-    <RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math"/>
    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\interpolate"/>
    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\linalg"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\burf"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\io"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\radar"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\satellite\cloudsat"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\satellite"/>
+    <RecentFolder Folder="D:\Working\MIScript\Jython\mis\satellite\gpm"/>
  </Path>
  <File>
    <OpenedFiles>
      <OpenedFile File="D:\MyProgram\java\MeteoInfoDev\toolbox\meteoview3d\_reload.py"/>
      <OpenedFile File="D:\MyProgram\java\MeteoInfoDev\toolbox\meteoview3d\mainGUI.py"/>
-      <OpenedFile File="D:\Working\MIScript\Jython\mis\io\radar\radar_bz2_3d_1.py"/>
-      <OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\linalg\dot_4.py"/>
-      <OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\linalg\dot_1.py"/>
+      <OpenedFile File="D:\Working\MIScript\Jython\mis\io\radar\radar_bin_proj.py"/>
+      <OpenedFile File="D:\Working\MIScript\Jython\mis\satellite\cloudsat\CloudSAT_Swath.py"/>
+      <OpenedFile File="D:\Working\MIScript\Jython\mis\satellite\gpm\GPM_DPRGMI_3d.py"/>
    </OpenedFiles>
    <RecentFiles>
      <RecentFile File="D:\MyProgram\java\MeteoInfoDev\toolbox\meteoview3d\_reload.py"/>
      <RecentFile File="D:\MyProgram\java\MeteoInfoDev\toolbox\meteoview3d\mainGUI.py"/>
-      <RecentFile File="D:\Working\MIScript\Jython\mis\io\radar\radar_bz2_3d_1.py"/>
-      <RecentFile File="D:\Working\MIScript\Jython\mis\common_math\linalg\dot_4.py"/>
-      <RecentFile File="D:\Working\MIScript\Jython\mis\common_math\linalg\dot_1.py"/>
+      <RecentFile File="D:\Working\MIScript\Jython\mis\io\radar\radar_bin_proj.py"/>
+      <RecentFile File="D:\Working\MIScript\Jython\mis\satellite\cloudsat\CloudSAT_Swath.py"/>
+      <RecentFile File="D:\Working\MIScript\Jython\mis\satellite\gpm\GPM_DPRGMI_3d.py"/>
    </RecentFiles>
  </File>
  <Font>
--- a/meteoinfo-lab/pylib/mipylib/init$py.class
+++ b/meteoinfo-lab/pylib/mipylib/init$py.class
--- a/meteoinfo-lab/pylib/mipylib/dataset/midata$py.class
+++ b/meteoinfo-lab/pylib/mipylib/dataset/midata$py.class
--- a/meteoinfo-lab/pylib/mipylib/dataset/midata.py
+++ b/meteoinfo-lab/pylib/mipylib/dataset/midata.py
@ -118,6 +118,8 @@ def addfile(fname, access='r', dtype='netcdf', keepopen=False, **kwargs):
        meteodata.openData(fname, keepopen)
        if meteodata.getDataInfo().getDataType() == MeteoDataType.RADAR:
            datafile = RadarDataFile(meteodata, access=access)
+        elif meteodata.getDataInfo().isRadial():
+            datafile = RadarDataFile(meteodata, access=access)
        else:
            datafile = DimDataFile(meteodata, access=access)
        return datafile
--- a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class
+++ b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric$py.class
--- a/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py
+++ b/meteoinfo-lab/pylib/mipylib/numeric/core/numeric.py
@ -42,7 +42,7 @@ __all__ = [
    'arctan2','atan2','ave_month','average','histogram','broadcast_to','cdiff','ceil',
    'concatenate','conj','conjugate','corrcoef','cos','cosh','cylinder','degrees','delete','delnan','diag','diff',
    'datatable','dot','empty','empty_like','exp','eye','flatnonzero','floor',
-    'fmax','fmin','full','hcurl','hdivg','hstack','identity','indices','interp2d','interpn','isarray',
+    'fmax','fmin','full','hcurl','hdivg','hstack','hypot','identity','indices','interp2d','interpn','isarray',
    'isclose','isfinite','isinf','isnan','isscalar','linspace','log','log10','logical_not','logspace',
    'magnitude','max','maximum','mean','median','meshgrid','min','minimum','monthname',
    'moveaxis','newaxis','ones','ones_like','outer','peaks','pol2cart','power','radians','reciprocal','reshape',
@ -1055,6 +1055,30 @@ def arctan2(x1, x2):
            return r
    else:
        return math.atan2(x1, x2)
+
+def hypot(x1, x2):
+    """
+    Given the “legs” of a right triangle, return its hypotenuse.
+
+    Equivalent to sqrt(x1**2 + x2**2), element-wise. If x1 or x2 is scalar_like (i.e., unambiguously
+    cast-able to a scalar type), it is broadcast for use with each element of the other argument.
+
+    Parameters
+    ----------
+    x1, x2 : array_like
+        Leg of the triangle(s). If x1.shape != x2.shape, they must be broadcastable to a common
+        shape (which becomes the shape of the output).
+
+    Returns
+    -------
+    z : array
+      The hypotenuse of the triangle(s). This is a scalar if both x1 and x2 are scalars.
+
+    """
+    x1 = asarray(x1)
+    x2 = asarray(x2)
+    z = sqrt(x1 * x1 + x2 * x2)
+    return z
        
 def exp(x):
    """
--- a/meteoinfo-lab/pylib/mipylib/numeric/lib/index_tricks$py.class
+++ b/meteoinfo-lab/pylib/mipylib/numeric/lib/index_tricks$py.class
--- a/meteoinfo-lab/pylib/mipylib/numeric/lib/index_tricks.py
+++ b/meteoinfo-lab/pylib/mipylib/numeric/lib/index_tricks.py
@ -3,7 +3,7 @@ import math
 from mipylib.numeric import core as _nx
 from ..core.numerictypes import ScalarType, find_common_type

-__all__ = ['r_','c_','mgrid']
+__all__ = ['r_','c_','mgrid','s_','index_exp']

 class nd_grid(object):
    """
@ -341,4 +341,73 @@ class CClass(AxisConcatenator):
        AxisConcatenator.__init__(self, -1, ndmin=2, trans1d=0)


-c_ = CClass()
+c_ = CClass()
+
+
+# You can do all this with slice() plus a few special objects,
+# but there's a lot to remember. This version is simpler because
+# it uses the standard array indexing syntax.
+#
+# Written by Konrad Hinsen <hinsen@cnrs-orleans.fr>
+# last revision: 1999-7-23
+#
+# Cosmetic changes by T. Oliphant 2001
+#
+#
+
+class IndexExpression:
+    """
+    A nicer way to build up index tuples for arrays.
+
+    .. note::
+       Use one of the two predefined instances `index_exp` or `s_`
+       rather than directly using `IndexExpression`.
+
+    For any index combination, including slicing and axis insertion,
+    ``a[indices]`` is the same as ``a[np.index_exp[indices]]`` for any
+    array `a`. However, ``np.index_exp[indices]`` can be used anywhere
+    in Python code and returns a tuple of slice objects that can be
+    used in the construction of complex index expressions.
+
+    Parameters
+    ----------
+    maketuple : bool
+        If True, always returns a tuple.
+
+    See Also
+    --------
+    index_exp : Predefined instance that always returns a tuple:
+       `index_exp = IndexExpression(maketuple=True)`.
+    s_ : Predefined instance without tuple conversion:
+       `s_ = IndexExpression(maketuple=False)`.
+
+    Notes
+    -----
+    You can do all this with `slice()` plus a few special objects,
+    but there's a lot to remember and this version is simpler because
+    it uses the standard array indexing syntax.
+
+    Examples
+    --------
+    >>> np.s_[2::2]
+    slice(2, None, 2)
+    >>> np.index_exp[2::2]
+    (slice(2, None, 2),)
+
+    >>> np.array([0, 1, 2, 3, 4])[np.s_[2::2]]
+    array([2, 4])
+
+    """
+
+    def __init__(self, maketuple):
+        self.maketuple = maketuple
+
+    def __getitem__(self, item):
+        if self.maketuple and not isinstance(item, tuple):
+            return (item,)
+        else:
+            return item
+
+
+index_exp = IndexExpression(maketuple=True)
+s_ = IndexExpression(maketuple=False)