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 cross_section function in meteolib package

Browse Source
This commit is contained in:
wyq 2024-08-08 15:27:03 +08:00
parent c7605e7b15
commit e4004a2a4a
20 changed files with 283 additions and 375 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/GraphicFactory.java
View File

@ -6216,6 +6216,10 @@ public class GraphicFactory {
if (Double.isNaN(v)) {
continue;
}
pb = (PolygonBreak) ls.findLegendBreak(v);
if (pb == null) {
continue;
}
x1 = x_s.getDouble(i * colNum + j);
x2 = x_s.getDouble(i * colNum + j + 1);
x3 = x_s.getDouble((i + 1) * colNum + j);
@ -6228,7 +6232,6 @@ public class GraphicFactory {
points.add(new PointD(x2, y_s.getDouble(i * colNum + j + 1)));
points.add((PointD) points.get(0).clone());
ps.setPoints(points);
pb = (PolygonBreak) ls.findLegendBreak(v);
Graphic graphic = new Graphic(ps, pb);
gc.add(graphic);
}

3
meteoinfo-chart/src/main/java/org/meteoinfo/chart/graphic/MeshGraphic.java
View File

@ -3,6 +3,7 @@ package org.meteoinfo.chart.graphic;
import org.joml.Vector3f;
import org.meteoinfo.chart.jogl.Transform;
import org.meteoinfo.common.Extent3D;
import org.meteoinfo.geometry.legend.ColorBreak;
import org.meteoinfo.geometry.legend.LegendManage;
import org.meteoinfo.geometry.colors.TransferFunction;
import org.meteoinfo.geometry.graphic.GraphicCollection3D;
@ -329,7 +330,7 @@ public class MeshGraphic extends GraphicCollection3D {
if (Float.isNaN(vertexValue[i])) {
color = legendScheme.getLegendBreak(0).getColor().getRGBComponents(null);
} else {
color = legendScheme.findLegendBreak(vertexValue[i]).getColor().getRGBComponents(null);
color = legendScheme.findLegendBreakAlways(vertexValue[i]).getColor().getRGBComponents(null);
}
System.arraycopy(color, 0, vertexColor, i * 4, 4);
}

11
meteoinfo-data/src/main/java/org/meteoinfo/data/dimarray/Dimension.java
View File

@ -254,9 +254,9 @@ public class Dimension {
}
/**
* Get dimension identifer
* Get dimension identifier
*
* @return Dimension identifer
* @return Dimension identifier
*/
public int getDimId() {
return dimId;
@ -482,6 +482,13 @@ public class Dimension {
return BigDecimalUtil.sub(dimValue.getDouble(1), dimValue.getDouble(0));
}
/**
* Flip the dimension data array
*/
public void flip() {
this.dimValue = this.dimValue.flip(0).copy();
}
/**
* Extract dimension

55
meteoinfo-geometry/src/main/java/org/meteoinfo/geometry/legend/LegendScheme.java
View File

@ -53,8 +53,8 @@ package org.meteoinfo.geometry.legend;
private ExtendFraction extendFraction = ExtendFraction.NONE;
private List<ColorBreak> legendBreaks;
private boolean hasNoData;
private double minValue;
private double maxValue;
private double minValue = -Double.MAX_VALUE;
private double maxValue = Double.MAX_VALUE;
private double undef;
private Map<Object, ColorBreak> uniqueValueMap;
private ColorMap colorMap;
@ -549,7 +549,7 @@ package org.meteoinfo.geometry.legend;
* @param v Value
* @return Legend break
*/
public ColorBreak findLegendBreak(Number v){
public ColorBreak findLegendBreakAlways(Number v){
switch (this.legendType) {
case SINGLE_SYMBOL:
return this.legendBreaks.get(0);
@ -583,12 +583,55 @@ package org.meteoinfo.geometry.legend;
}
}
/**
* Find legend break by value
* @param v Value
* @return Legend break
*/
public ColorBreak findLegendBreak(Number v){
switch (this.legendType) {
case SINGLE_SYMBOL:
return this.legendBreaks.get(0);
case UNIQUE_VALUE:
if (this.uniqueValueMap == null || this.uniqueValueMap.size() != this.legendBreaks.size())
this.updateUniqueValueMap();
if (this.uniqueValueMap.containsKey(v)) {
return this.uniqueValueMap.get(v);
} else {
return null;
}
default:
if (v.doubleValue() < this.minValue || v.doubleValue() > this.maxValue) {
return null;
}
double sv, ev;
for (ColorBreak cb : this.legendBreaks){
sv = Double.parseDouble(cb.getStartValue().toString());
ev = Double.parseDouble(cb.getEndValue().toString());
if (sv == ev){
if (v.doubleValue() == sv)
return cb;
} else {
if (v.doubleValue() >= sv && v.doubleValue() < ev){
return cb;
}
}
}
return this.legendBreaks.get(this.getBreakNum() - 1);
}
}
/**
* Get legend break index by value
* @param v Value
* @return Legend break index
*/
public int legendBreakIndex(double v) {
if (v < this.minValue || v > this.maxValue) {
return -1;
}
double sv, ev;
for (int i = 0; i < this.legendBreaks.size(); i ++){
ColorBreak cb = this.legendBreaks.get(i);
@ -603,11 +646,7 @@ package org.meteoinfo.geometry.legend;
}
}
}
if (v >= this.getMaxValue())
return this.getBreakNum() - 1;
else
return 0;
return this.getBreakNum() - 1;
}
/**

36
meteoinfo-lab/milconfig.xml
View File

@ -1,32 +1,32 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<MeteoInfo File="milconfig.xml" Type="configurefile">
<Path OpenPath="D:\Working\MIScript\Jython\mis\map\projection">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\contour"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\traj"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\map\maskout"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\hdf"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\grib"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\netcdf"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\grads"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
<Path OpenPath="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\surf">
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\calc"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io"/>
<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\projection"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\radar"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\interpolation"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\slice"/>
<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\plot_types\3d\jogl"/>
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\surf"/>
</Path>
<File>
<OpenedFiles>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\plot\plot_cdata_3.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\io\radar\radar_cma_base_colorbar.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\map\projection\aeqd_proj_2.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\meteo\interpolation\cross_section_1.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\slice\slice_2d.py"/>
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\interpolate\griddata_kriging_lcc.py"/>
</OpenedFiles>
<RecentFiles>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\plot\plot_cdata_3.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\io\radar\radar_cma_base_colorbar.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\map\projection\aeqd_proj_2.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\meteo\interpolation\cross_section_1.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\3d\jogl\slice\slice_2d.py"/>
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\interpolate\griddata_kriging_lcc.py"/>
</RecentFiles>
</File>
<Font>

BIN
meteoinfo-lab/pylib/mipylib/dataset/dimvariable$py.class
View File

Binary file not shown.

5
meteoinfo-lab/pylib/mipylib/dataset/dimvariable.py
View File

@ -112,7 +112,7 @@ class DimVariable(object):
indices = tuple(indices)
if len(indices) != self.ndim:
print 'indices must be ' + str(self.ndim) + ' dimensions!'
print('indices must be ' + str(self.ndim) + ' dimensions!')
return None
if not self.proj is None and not self.proj.isLonLat():
@ -253,8 +253,9 @@ class DimVariable(object):
if sidx > eidx and step > 0:
step = -step
else:
print k
print(k)
return None
if isrange:
if sidx >= dimlen or eidx >= dimlen:
print('Index out of range!')

35
meteoinfo-lab/pylib/mipylib/geolib/_graphic.py
View File

@ -12,13 +12,20 @@ class GeoGraphicCollection(object):
:param geographic: (*JGeoGraphicCollection*) Java GeoGraphicCollection object.
"""
self._geographic = geographic
self.shapetype = geographic.getShapeType()
self.proj = geographic.getProjInfo()
def get_graphics(self):
def get_graphics(self, xshift=0):
"""
Get graphics.
:param xshift: (*float*) X shift
:return: Graphics.
"""
return self._geographic
if xshift == 0:
return self._geographic
else:
return self._geographic.xshift(xshift)
@property
def shapes(self):
@ -82,6 +89,30 @@ class GeoGraphicCollection(object):
"""
self._geographic.moveLabel(label, x, y)
@property
def visible(self):
return self._geographic.isVisible()
@visible.setter
def visible(self, val):
self._geographic.setVisible(val)
@property
def legend(self):
"""
Get legend scheme.
"""
return self._geographic.getLegendScheme()
@legend.setter
def legend(self, legend):
"""
Set legend scheme.
:param legend: (*LegendScheme*) Legend scheme.
"""
self._geographic.setLegendScheme(legend)
def update_legend(self, ltype, fieldname):
"""
Update legend scheme.

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

Binary file not shown.

4
meteoinfo-lab/pylib/mipylib/meteolib/interpolate/__init__.py
View File

@ -1,4 +1,6 @@
from .one_dimension import *
from slices import *
__all__ = []
__all__ += one_dimension.__all__
__all__ += one_dimension.__all__
__all__ += slices.__all__

52
meteoinfo-lab/pylib/mipylib/meteolib/interpolate/slices.py Normal file
View File

@ -0,0 +1,52 @@
"""
Tools for interpolating to a vertical slice/cross section through data.
"""
from org.meteoinfo.math.interpolate import InterpUtil, InterpolationMethod
from mipylib.numeric.core import NDArray
__all__ = ['cross_section']
def cross_section(data, start, end, steps=100, interp_type='linear'):
r"""Obtain an interpolated cross-sectional slice through gridded data.
This function takes a vertical cross-sectional slice along a geodesic through the given data
on a regular grid, which is given as an `DimArray` so that we can utilize its coordinate and
projection metadata.
Parameters
----------
data: `DimArray`
Three- (or higher) dimensional field(s) to interpolate. The DimArray must include both an x and
y coordinate dimension.
start: (2, ) array-like
A latitude-longitude pair designating the start point of the cross section (units are
degrees north and degrees east).
end: (2, ) array-like
A latitude-longitude pair designating the end point of the cross section (units are
degrees north and degrees east).
steps: int, optional
The number of points along the geodesic between the start and the end point
(including the end points) to use in the cross section. Defaults to 100.
interp_type: str, optional
The interpolation method, either 'linear' or 'nearest'. Defaults to 'linear'.
Returns
-------
List of `NDArray`
The interpolated cross section, new x coordinates and new y coordinates
"""
x = data.dimvalue(-1)
y = data.dimvalue(-2)
interp_type = InterpolationMethod.valueOf(interp_type.upper())
xyslice = start + end
if data.ndim == 2:
rs = InterpUtil.sliceXY(x._array, y._array, data._array, xyslice, steps, interp_type)
else:
z = data.dimvalue(-3)
rs = InterpUtil.sliceXY(x._array, y._array, z._array, data._array, xyslice, steps, interp_type)
return NDArray(rs[0]), NDArray(rs[1]), NDArray(rs[2])

BIN
meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes$py.class
View File

Binary file not shown.

12
meteoinfo-lab/pylib/mipylib/plotlib/_mapaxes.py
View File

@ -507,7 +507,7 @@ class MapAxes(Axes):
encoding = kwargs.pop('encoding', None)
layer = migeo.georead(fn, encoding)
islayer = True
elif isinstance(args[0], MILayer):
elif isinstance(args[0], (MILayer, GeoGraphicCollection)):
layer = args[0]
islayer = True
@ -516,7 +516,7 @@ class MapAxes(Axes):
layer.visible = visible
xshift = kwargs.pop('xshift', 0)
zorder = kwargs.pop('zorder', None)
if layer.layer_type == LayerTypes.IMAGE_LAYER:
if isinstance(layer, MILayer) and layer.layer_type == LayerTypes.IMAGE_LAYER:
interpolation = kwargs.pop('interpolation', None)
graphics = layer.get_graphics(xshift, interpolation)
else:
@ -1186,7 +1186,7 @@ class MapAxes(Axes):
:param zorder: (*int*) Z-order of created layer for display.
:param select: (*boolean*) Set the return layer as selected layer or not.
:returns: (*VectoryLayer*) Polygon VectoryLayer created from array data.
:returns: (*Graphics*) Polygon graphics created from array data.
"""
n = len(args)
if n <= 2:
@ -1202,8 +1202,10 @@ class MapAxes(Axes):
if a.ndim == 2 and x.ndim == 1:
x, y = np.meshgrid(x, y)
ls = plotutil.getlegendscheme(args, a.min(), a.max(), **kwargs)
vmin = kwargs.pop('vmin', a.min())
vmax = kwargs.pop('vmax', a.max())
ls = plotutil.getlegendscheme(args, vmin, vmax, **kwargs)
ls = ls.convertTo(ShapeTypes.POLYGON)
if not kwargs.has_key('edgecolor'):
kwargs['edgecolor'] = None

228
meteoinfo-math/src/main/java/org/meteoinfo/math/interpolate/InterpUtil.java
View File

@ -2186,147 +2186,12 @@ public class InterpUtil {
return r;
}
/**
* Get value index in a dimension array
*
* @param dim Dimension array
* @param v The value
* @return value index
*/
public static int getDimIndex(Array dim, Number v) {
dim = dim.copyIfView();
switch (dim.getDataType()) {
case BYTE:
return Arrays.binarySearch((byte[]) dim.getStorage(), v.byteValue());
case INT:
return Arrays.binarySearch((int[]) dim.getStorage(), v.intValue());
case SHORT:
return Arrays.binarySearch((short[]) dim.getStorage(), v.shortValue());
case LONG:
return Arrays.binarySearch((long[]) dim.getStorage(), v.longValue());
case FLOAT:
return Arrays.binarySearch((float[]) dim.getStorage(), v.floatValue());
case DOUBLE:
return Arrays.binarySearch((double[]) dim.getStorage(), v.doubleValue());
}
int n = (int) dim.getSize();
if (v.doubleValue() < dim.getDouble(0) || v.doubleValue() > dim.getDouble(n - 1)) {
return -1;
}
int idx = n - 1;
for (int i = 1; i < n; i++) {
if (v.doubleValue() < dim.getDouble(i)) {
idx = i - 1;
break;
}
}
return idx;
}
/**
* Get grid array x/y value index
* @param xdim X coordinate array
* @param ydim Y coordinate array
* @param x X value
* @param y Y value
* @return X/Y index
*/
public static int[] gridIndex(Array xdim, Array ydim, double x, double y) {
if (xdim.getRank() == 1) {
int xn = (int) xdim.getSize();
int yn = (int) ydim.getSize();
int xIdx = getDimIndex(xdim, x);
if (xIdx == -1 || xIdx == -(xn + 1)) {
return null;
} else if (xIdx < 0) {
xIdx = -xIdx - 2;
}
int yIdx = getDimIndex(ydim, y);
if (yIdx == -1 || yIdx == -(yn + 1)) {
return null;
} else if (yIdx < 0) {
yIdx = -yIdx - 2;
}
if (xIdx == xn - 1) {
xIdx = xn - 2;
}
if (yIdx == yn - 1) {
yIdx = yn - 2;
}
return new int[]{yIdx, xIdx};
} else {
int xIdx = -1, yIdx = -1;
int[] shape = xdim.getShape();
int yn = shape[0];
int xn = shape[1];
Index index = new Index2D(shape);
double x1, x2, y1, y2;
for (int i = 0; i < yn - 1; i++) {
for (int j = 0; j < xn - 1; j++) {
index = index.set(i, j);
y1 = ydim.getDouble(index);
index = index.set(i+1, j);
y2 = ydim.getDouble(index);
if (y >= y1 && y < y2) {
index = index.set(i, j);
x1 = xdim.getDouble(index);
index = index.set(i, j+1);
x2 = xdim.getDouble(index);
if (x >= x1 && x < x2) {
yIdx = i;
xIdx = j;
}
}
}
}
if (yIdx >= 0 && xIdx >= 0)
return new int[]{yIdx, xIdx};
else
return null;
}
}
/**
* Get grid array x/y value index
* @param xdim X coordinate array
* @param ydim Y coordinate array
* @param x X value
* @param y Y value
* @return X/Y index
*/
public static int[] gridIndex(double[][] xdim, double[][] ydim, double x, double y) {
int xIdx = -1, yIdx = -1;
int yn = xdim.length;
int xn = xdim[0].length;
for (int i = 0; i < yn - 1; i++) {
for (int j = 0; j < xn - 1; j++) {
if (y >= ydim[i][j] && y < ydim[i+1][j]) {
if (x >= xdim[i][j] && x < xdim[i][j+1]) {
yIdx = i;
xIdx = j;
}
}
}
}
if (yIdx >= 0 && xIdx >= 0)
return new int[]{yIdx, xIdx};
else
return null;
}
private static double bilinear(Array data, Index dindex, Array xdim, Array ydim, double x, double y) {
xdim = xdim.copyIfView();
ydim = ydim.copyIfView();
double iValue = Double.NaN;
int[] xyIdx = gridIndex(xdim, ydim, x, y);
int[] xyIdx = ArrayUtil.gridIndex(xdim, ydim, x, y);
if (xyIdx == null) {
return iValue;
}
@ -2390,7 +2255,7 @@ public class InterpUtil {
xdim = xdim.copyIfView();
ydim = ydim.copyIfView();
int[] xyIdx = gridIndex(xdim, ydim, x, y);
int[] xyIdx = ArrayUtil.gridIndex(xdim, ydim, x, y);
if (xyIdx == null) {
return Double.NaN;
}
@ -2510,10 +2375,6 @@ public class InterpUtil {
InterpolationMethod method) {
xa = xa.copyIfView();
ya = ya.copyIfView();
za = za.copyIfView();
data = data.copyIfView();
RectInterpolator3D interpolator3D = RectInterpolator3D.factory(xa, ya, za, data, method);
double x1 = xySlice.get(0).doubleValue();
double y1 = xySlice.get(1).doubleValue();
@ -2532,14 +2393,87 @@ public class InterpUtil {
int xn = dx == 0 ? 1 : (int) Math.ceil((x2 - x1) / dx);
int yn = dy == 0 ? 1 : (int) Math.ceil(Math.abs(y2 - y1) / dy);
int rn = Math.max(xn, yn);
Array xs = ArrayUtil.lineSpace(x1, x2, rn, true);
Array ys = ArrayUtil.lineSpace(y1, y2, rn, true);
int zn = (int)za.getSize();
Array xs2d = ArrayUtil.repeat(xs.reshape(new int[]{1, rn}), Arrays.asList(zn), 0);
Array ys2d = ArrayUtil.repeat(ys.reshape(new int[]{1, rn}), Arrays.asList(zn), 0);
Array zs2d = ArrayUtil.repeat(za.reshape(new int[]{zn, 1}), Arrays.asList(rn), 1);
return sliceXY(xa, ya, za, data, xySlice, rn, method);
}
/**
* Slice 3D array data by x/y cross-section
* @param xa X coordinate array - 1D
* @param ya Y coordinate array - 1D
* @param za Z coordinate array - 1D
* @param data Data array - 3D
* @param xySlice X/Y slice points - [x1, y1, x2, y2]
* @param steps Number of points
* @param method Interpolation method - nearest or linear
* @return x/y cross-section slice data array
*/
public static Array[] sliceXY(Array xa, Array ya, Array za, Array data, List<Number> xySlice,
int steps, InterpolationMethod method) {
xa = xa.copyIfView();
ya = ya.copyIfView();
za = za.copyIfView();
data = data.copyIfView();
RectInterpolator3D interpolator3D = RectInterpolator3D.factory(xa, ya, za, data, method);
double x1 = xySlice.get(0).doubleValue();
double y1 = xySlice.get(1).doubleValue();
double x2 = xySlice.get(2).doubleValue();
double y2 = xySlice.get(3).doubleValue();
if (x1 > x2) {
double temp = x2;
x2 = x1;
x1 = temp;
temp = y2;
y2 = y1;
y1 = temp;
}
Array xs = ArrayUtil.lineSpace(x1, x2, steps, true);
Array ys = ArrayUtil.lineSpace(y1, y2, steps, true);
int zn = (int) za.getSize();
Array xs2d = ArrayUtil.repeat(xs.reshape(new int[]{1, steps}), Arrays.asList(zn), 0);
Array ys2d = ArrayUtil.repeat(ys.reshape(new int[]{1, steps}), Arrays.asList(zn), 0);
Array zs2d = ArrayUtil.repeat(za.reshape(new int[]{zn, 1}), Arrays.asList(steps), 1);
Array r = interpolator3D.interpolate(xs2d, ys2d, zs2d);
return new Array[]{r, xs, ys, za, xs2d, ys2d, zs2d};
}
/**
* Slice 2D array data by x/y cross-section
* @param xa X coordinate array - 1D
* @param ya Y coordinate array - 1D
* @param data Data array - 2D
* @param xySlice X/Y slice points - [x1, y1, x2, y2]
* @param steps Number of points
* @param method Interpolation method - nearest or linear
* @return x/y cross-section slice data array
*/
public static Array[] sliceXY(Array xa, Array ya, Array data, List<Number> xySlice,
int steps, InterpolationMethod method) {
xa = xa.copyIfView();
ya = ya.copyIfView();
data = data.copyIfView();
RectInterpolator interpolator = RectInterpolator.factory(xa, ya, data, method);
double x1 = xySlice.get(0).doubleValue();
double y1 = xySlice.get(1).doubleValue();
double x2 = xySlice.get(2).doubleValue();
double y2 = xySlice.get(3).doubleValue();
if (x1 > x2) {
double temp = x2;
x2 = x1;
x1 = temp;
temp = y2;
y2 = y1;
y1 = temp;
}
Array xs = ArrayUtil.lineSpace(x1, x2, steps, true);
Array ys = ArrayUtil.lineSpace(y1, y2, steps, true);
Array r = interpolator.interpolate(xs, ys);
return new Array[]{r, xs, ys};
}
}

111
meteoinfo-math/src/main/java/org/meteoinfo/math/interpolate/RectInterpolator.java
View File

@ -26,108 +26,19 @@ public abstract class RectInterpolator {
}
/**
* Get value index in a dimension array
*
* @param dim Dimension array
* @param v The value
* @return value index
* Factory method
* @param xa X coordinate array - 1D
* @param ya Y coordinate array - 1D
* @param va Value array - 2D or more than 2D
* @param method Interpolation method
* @return RectInterpolator
*/
public static int getDimIndex(Array dim, Number v) {
dim = dim.copyIfView();
switch (dim.getDataType()) {
case BYTE:
return Arrays.binarySearch((byte[]) dim.getStorage(), v.byteValue());
case INT:
return Arrays.binarySearch((int[]) dim.getStorage(), v.intValue());
case SHORT:
return Arrays.binarySearch((short[]) dim.getStorage(), v.shortValue());
case LONG:
return Arrays.binarySearch((long[]) dim.getStorage(), v.longValue());
case FLOAT:
return Arrays.binarySearch((float[]) dim.getStorage(), v.floatValue());
case DOUBLE:
return Arrays.binarySearch((double[]) dim.getStorage(), v.doubleValue());
}
int n = (int) dim.getSize();
if (v.doubleValue() < dim.getDouble(0) || v.doubleValue() > dim.getDouble(n - 1)) {
return -1;
}
int idx = n - 1;
for (int i = 1; i < n; i++) {
if (v.doubleValue() < dim.getDouble(i)) {
idx = i - 1;
break;
}
}
return idx;
}
/**
* Get grid array x/y value index
* @param xdim X coordinate array
* @param ydim Y coordinate array
* @param x X value
* @param y Y value
* @return X/Y index
*/
public static int[] gridIndex(Array xdim, Array ydim, double x, double y) {
if (xdim.getRank() == 1) {
int xn = (int) xdim.getSize();
int yn = (int) ydim.getSize();
int xIdx = getDimIndex(xdim, x);
if (xIdx == -1 || xIdx == -(xn + 1)) {
return null;
} else if (xIdx < 0) {
xIdx = -xIdx - 2;
}
int yIdx = getDimIndex(ydim, y);
if (yIdx == -1 || yIdx == -(yn + 1)) {
return null;
} else if (yIdx < 0) {
yIdx = -yIdx - 2;
}
if (xIdx == xn - 1) {
xIdx = xn - 2;
}
if (yIdx == yn - 1) {
yIdx = yn - 2;
}
return new int[]{yIdx, xIdx};
public static RectInterpolator factory(Array xa, Array ya, Array va,
InterpolationMethod method) {
if (method == InterpolationMethod.NEAREST) {
return new RectNearestInterpolator(xa, ya, va);
} else {
int xIdx = -1, yIdx = -1;
int[] shape = xdim.getShape();
int yn = shape[0];
int xn = shape[1];
Index index = new Index2D(shape);
double x1, x2, y1, y2;
for (int i = 0; i < yn - 1; i++) {
for (int j = 0; j < xn - 1; j++) {
index = index.set(i, j);
y1 = ydim.getDouble(index);
index = index.set(i+1, j);
y2 = ydim.getDouble(index);
if (y >= y1 && y < y2) {
index = index.set(i, j);
x1 = xdim.getDouble(index);
index = index.set(i, j+1);
x2 = xdim.getDouble(index);
if (x >= x1 && x < x2) {
yIdx = i;
xIdx = j;
}
}
}
}
if (yIdx >= 0 && xIdx >= 0)
return new int[]{yIdx, xIdx};
else
return null;
return new RectLinearInterpolator(xa, ya, va);
}
}

46
meteoinfo-math/src/main/java/org/meteoinfo/math/interpolate/RectInterpolator3D.java
View File

@ -4,6 +4,7 @@ import org.meteoinfo.ndarray.Array;
import org.meteoinfo.ndarray.DataType;
import org.meteoinfo.ndarray.Index;
import org.meteoinfo.ndarray.Index3D;
import org.meteoinfo.ndarray.math.ArrayUtil;
import java.util.Arrays;
@ -45,45 +46,6 @@ public abstract class RectInterpolator3D {
}
}
/**
* Get value index in a dimension array
*
* @param dim Dimension array
* @param v The value
* @return value index
*/
public static int getDimIndex(Array dim, Number v) {
dim = dim.copyIfView();
switch (dim.getDataType()) {
case BYTE:
return Arrays.binarySearch((byte[]) dim.getStorage(), v.byteValue());
case INT:
return Arrays.binarySearch((int[]) dim.getStorage(), v.intValue());
case SHORT:
return Arrays.binarySearch((short[]) dim.getStorage(), v.shortValue());
case LONG:
return Arrays.binarySearch((long[]) dim.getStorage(), v.longValue());
case FLOAT:
return Arrays.binarySearch((float[]) dim.getStorage(), v.floatValue());
case DOUBLE:
return Arrays.binarySearch((double[]) dim.getStorage(), v.doubleValue());
}
int n = (int) dim.getSize();
if (v.doubleValue() < dim.getDouble(0) || v.doubleValue() > dim.getDouble(n - 1)) {
return -1;
}
int idx = n - 1;
for (int i = 1; i < n; i++) {
if (v.doubleValue() < dim.getDouble(i)) {
idx = i - 1;
break;
}
}
return idx;
}
/**
* Get grid array x/y value index
* @param xdim X coordinate array
@ -100,21 +62,21 @@ public abstract class RectInterpolator3D {
int xn = (int) xdim.getSize();
int yn = (int) ydim.getSize();
int zn = (int) zdim.getSize();
int xIdx = getDimIndex(xdim, x);
int xIdx = ArrayUtil.getDimIndex(xdim, x);
if (xIdx == -1 || xIdx == -(xn + 1)) {
return null;
} else if (xIdx < 0) {
xIdx = -xIdx - 2;
}
int yIdx = getDimIndex(ydim, y);
int yIdx = ArrayUtil.getDimIndex(ydim, y);
if (yIdx == -1 || yIdx == -(yn + 1)) {
return null;
} else if (yIdx < 0) {
yIdx = -yIdx - 2;
}
int zIdx = getDimIndex(zdim, z);
int zIdx = ArrayUtil.getDimIndex(zdim, z);
if (zIdx == -1 || zIdx == -(zn + 1)) {
return null;
} else if (zIdx < 0) {

2
meteoinfo-math/src/main/java/org/meteoinfo/math/interpolate/RectLinearInterpolator.java
View File

@ -24,7 +24,7 @@ public class RectLinearInterpolator extends RectInterpolator{
@Override
double cellValue(Index dindex, double x, double y) {
double iValue = Double.NaN;
int[] xyIdx = gridIndex(xa, ya, x, y);
int[] xyIdx = ArrayUtil.gridIndex(xa, ya, x, y);
if (xyIdx == null) {
return iValue;
}

5
meteoinfo-math/src/main/java/org/meteoinfo/math/interpolate/RectNearestInterpolator.java
View File

@ -2,6 +2,7 @@ package org.meteoinfo.math.interpolate;
import org.meteoinfo.ndarray.Array;
import org.meteoinfo.ndarray.Index;
import org.meteoinfo.ndarray.math.ArrayUtil;
public class RectNearestInterpolator extends RectInterpolator{
/**
@ -17,7 +18,7 @@ public class RectNearestInterpolator extends RectInterpolator{
@Override
double cellValue(Index dindex, double x, double y) {
int[] xyIdx = gridIndex(xa, ya, x, y);
int[] xyIdx = ArrayUtil.gridIndex(xa, ya, x, y);
if (xyIdx == null) {
return Double.NaN;
}
@ -49,7 +50,7 @@ public class RectNearestInterpolator extends RectInterpolator{
@Override
double interpolate(double x, double y) {
int[] xyIdx = gridIndex(xa, ya, x, y);
int[] xyIdx = ArrayUtil.gridIndex(xa, ya, x, y);
if (xyIdx == null) {
return Double.NaN;
}

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

@ -5053,20 +5053,17 @@ public class ArrayUtil {
if (idx < n - 1) {
index.setDim(idx, indices[idx]);
iterIndex(ii, index, indices, idx + 1);
if (indices[idx] < index.getShape(idx) - 1)
if (indices[idx] < index.getShape(idx) - 1) {
index.setDim(idx, indices[idx] + 1);
else
index.setDim(idx, indices[idx]);
iterIndex(ii, index, indices, idx + 1);
iterIndex(ii, index, indices, idx + 1);
}
} else {
index.setDim(idx, indices[idx]);
ii.add((Index) index.clone());
//System.out.println(index);
if (indices[idx] < index.getShape(idx) - 1)
index.setDim(idx, indices[idx] + 1);
else
index.setDim(idx, indices[idx]);
ii.add((Index) index.clone());
ii.add((Index) index.clone());
//System.out.println(index);
}
}

37
meteoinfo-projection/src/main/java/org/meteoinfo/projection/info/AzimuthEquidistant.java
View File

@ -56,41 +56,6 @@ public class AzimuthEquidistant extends ProjectionInfo {
// </editor-fold>
// <editor-fold desc="Methods">
@Override
public void updateBoundary() {
double epsilon = 1e-10;
double cenLon = this.getCenterLon();
double minLon = cenLon - 180 + epsilon;
double maxLon = cenLon + 180 - epsilon;
double minLat = -90;
double maxLat = 90;
List<PointD> points = new ArrayList<>();
double lon = minLon;
double lat = minLat;
while (lon < maxLon) {
points.add(new PointD(lon, lat));
lon += 1;
}
lon = maxLon;
while (lat < maxLat) {
points.add(new PointD(lon, lat));
lat += 1;
}
lat = maxLat;
while (lon > minLon) {
points.add(new PointD(lon, lat));
lon -= 1;
}
lon = minLon;
while (lat > minLat) {
points.add(new PointD(lon, lat));
lat -= 1;
}
lat = minLat;
points.add(new PointD(lon, lat));
PolygonShape ps = new PolygonShape();
ps.setPoints(points);
this.boundary = ProjectionUtil.projectPolygonShape(ps, KnownCoordinateSystems.geographic.world.WGS1984, this);
}
// </editor-fold>
}