Skip to content

Commit abfa2d8

Browse files
TomBlockTomBlock
committed
implemented 2d/3d interpolation
1 parent 6494ece commit abfa2d8

File tree

7 files changed

+185
-103
lines changed

7 files changed

+185
-103
lines changed

‎post-processing-tool/src/main/java/com/bc/fiduceo/post/plugin/era5/BilinearInterpolator.java‎

Lines changed: 16 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,10 +4,14 @@ class BilinearInterpolator {
44

55
private final double a;
66
private final double b;
7+
private final int xMin;
8+
private final int yMin;
79

8-
BilinearInterpolator(double a, double b) {
10+
BilinearInterpolator(double a, double b, int xMin, int yMin) {
911
this.a = a;
1012
this.b = b;
13+
this.xMin = xMin;
14+
this.yMin = yMin;
1115
}
1216

1317
double interpolate(float c00, float c10, float c01, float c11) {
@@ -16,7 +20,17 @@ class BilinearInterpolator {
1620
return lerp(interp0, interp1, b);
1721
}
1822

19-
static double lerp(double c0, double c1, double t) {
23+
int getXMin() {
24+
return xMin;
25+
}
26+
27+
int getYMin() {
28+
return yMin;
29+
}
30+
31+
private static double lerp(double c0, double c1, double t) {
2032
return c0 + (c1 - c0) * t;
2133
}
34+
35+
// probably clever to add the reading coordinates in the rea-5 raster
2236
}

‎post-processing-tool/src/main/java/com/bc/fiduceo/post/plugin/era5/Era5PostProcessing.java‎

Lines changed: 27 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -30,30 +30,6 @@ class Era5PostProcessing extends PostProcessing {
3030
matchupFields = null;
3131
}
3232

33-
// package access for testing only tb 2020-11-17
34-
static GeoRect getGeoRegion(Array lonArray, Array latArray) {
35-
// @todo 1 tb/tb check for anti-meridian and handle 2020-11-17
36-
final float[] lonMinMax = getMinMax(lonArray);
37-
final float[] latMinMax = getMinMax(latArray);
38-
39-
return new GeoRect(lonMinMax[0], lonMinMax[1], latMinMax[0], latMinMax[1]);
40-
}
41-
42-
// package access for testing only tb 2020-11-17
43-
static Rectangle getEra5RasterPosition(GeoRect geoRect) {
44-
final float lonMin = geoRect.getLonMin();
45-
final int xMin = getEra5LonMin(lonMin);
46-
47-
final float lonMax = geoRect.getLonMax();
48-
final int xMax = getEra5LonMax(lonMax);
49-
50-
// remember: y axis runs top->down so we need to invert the coordinates tb 2020-11-17
51-
final int yMax = getEra5LatMax(geoRect.getLatMin());
52-
final int yMin = getEra5LatMin(geoRect.getLatMax());
53-
54-
return new Rectangle(xMin, yMin, xMax - xMin + 1, yMax - yMin + 1);
55-
}
56-
5733
private static int getEra5LatMin(float latMax) {
5834
final double shiftedLat = latMax + EPS;
5935
final double scaledLatMax = Math.ceil(shiftedLat * 4) / 4;
@@ -88,6 +64,10 @@ static InterpolationContext getInterpolationContext(Array lonArray, Array latArr
8864

8965
final Index lonIdx = lonArray.getIndex();
9066
final Index latIdx = latArray.getIndex();
67+
int xMin = Integer.MAX_VALUE;
68+
int xMax = Integer.MIN_VALUE;
69+
int yMin = Integer.MAX_VALUE;
70+
int yMax = Integer.MIN_VALUE;
9171
for (int y = 0; y < shape[0]; y++) {
9272
for (int x = 0; x < shape[1]; x++) {
9373
lonIdx.set(y, x);
@@ -101,16 +81,35 @@ static InterpolationContext getInterpolationContext(Array lonArray, Array latArr
10181
// + calculate latitude delta -> b
10282
// + create BilinearInterpolator(a, b)
10383
// + store to context at (x, y)
104-
final double era5LonMin = getEra5LonMin(lon) * 0.25 - 180.0;
105-
final double era5LatMin = 90.0 - getEra5LatMin(lat) * 0.25;
84+
final int era5_X_min = getEra5LonMin(lon);
85+
final int era5_Y_min = getEra5LatMin(lat);
86+
if (era5_X_min < xMin) {
87+
xMin = era5_X_min;
88+
}
89+
if (era5_X_min > xMax) {
90+
xMax = era5_X_min;
91+
}
92+
if (era5_Y_min < yMin) {
93+
yMin = era5_Y_min;
94+
}
95+
if (era5_Y_min > yMax) {
96+
yMax = era5_Y_min;
97+
}
98+
99+
final double era5LonMin = era5_X_min * 0.25 - 180.0;
100+
final double era5LatMin = 90.0 - era5_Y_min * 0.25;
106101

107102
// we have a quarter degree raster and need to normalize the distance tb 2020-11-20
108103
final double lonDelta = (lon - era5LonMin) * 4.0;
109104
final double latDelta = (era5LatMin - lat) * 4.0;
110105

111-
final BilinearInterpolator interpolator = new BilinearInterpolator(lonDelta, latDelta);
106+
final BilinearInterpolator interpolator = new BilinearInterpolator(lonDelta, latDelta, era5_X_min, era5_Y_min);
112107
context.set(x, y, interpolator);
113108
}
109+
110+
// add 2 to width and height to have always 4 points for the interpolation tb 2020-11-30
111+
final Rectangle era5Rect = new Rectangle(xMin, yMin, xMax - xMin + 2, yMax - yMin + 2);
112+
context.setEra5Region(era5Rect);
114113
}
115114
return context;
116115
}
@@ -137,7 +136,7 @@ private static float[] getMinMax(Array floatArray) {
137136

138137

139138
@Override
140-
protected void prepare(NetcdfFile reader, NetcdfFileWriter writer) throws IOException, InvalidRangeException {
139+
protected void prepare(NetcdfFile reader, NetcdfFileWriter writer) {
141140
final Dimension matchupCountDimension = reader.findDimension(FiduceoConstants.MATCHUP_COUNT);
142141
if (matchupCountDimension == null) {
143142
throw new RuntimeException("Expected dimension not present in file: " + FiduceoConstants.MATCHUP_COUNT);

‎post-processing-tool/src/main/java/com/bc/fiduceo/post/plugin/era5/InterpolationContext.java‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,13 @@
11
package com.bc.fiduceo.post.plugin.era5;
22

3+
import java.awt.*;
4+
35
class InterpolationContext {
46

57
private final BilinearInterpolator[][] interpolators;
68
private final int width;
79
private final int height;
10+
private Rectangle era5Region;
811

912
InterpolationContext(int x, int y) {
1013
width = x;
@@ -28,6 +31,14 @@ public void set(int x, int y, BilinearInterpolator interpolator) {
2831
interpolators[y][x] = interpolator;
2932
}
3033

34+
Rectangle getEra5Region() {
35+
return era5Region;
36+
}
37+
38+
void setEra5Region(Rectangle era5Region) {
39+
this.era5Region = era5Region;
40+
}
41+
3142
private void checkBoundaries(int x, int y) {
3243
if (x < 0 || x >= width || y < 0 || y >= height) {
3344
throw new IllegalArgumentException("Access interpolator out of raster: " + x + ", " + y);

‎post-processing-tool/src/main/java/com/bc/fiduceo/post/plugin/era5/SatelliteFields.java‎

Lines changed: 79 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,6 @@
11
package com.bc.fiduceo.post.plugin.era5;
22

33
import com.bc.fiduceo.FiduceoConstants;
4-
import com.bc.fiduceo.core.GeoRect;
54
import com.bc.fiduceo.reader.ReaderUtils;
65
import com.bc.fiduceo.util.NetCDFUtils;
76
import com.bc.fiduceo.util.TimeUtils;
@@ -108,34 +107,98 @@ void compute(Configuration config, NetcdfFile reader, NetcdfFileWriter writer) t
108107
// + prepare interpolation context
109108
final int numMatches = NetCDFUtils.getDimensionLength(FiduceoConstants.MATCHUP_COUNT, reader);
110109
final int[] nwpShape = getNwpShape(satFieldsConfig, lonArray.getShape());
110+
nwpShape[0] = 1; // we read matchup layer by layer
111111
final int[] nwpOffset = getNwpOffset(lonArray.getShape(), nwpShape);
112112

113-
final Index index = era5TimeArray.getIndex();
113+
final Index timeIndex = era5TimeArray.getIndex();
114114
for (int m = 0; m < numMatches; m++) {
115115
nwpOffset[0] = m;
116-
nwpShape[0] = 1; // @todo 1 tb/tb adapt to 3d variables 2020-11-24
116+
117117
final Array lonLayer = lonArray.section(nwpOffset, nwpShape).reduce();
118118
final Array latLayer = latArray.section(nwpOffset, nwpShape).reduce();
119119

120-
final GeoRect geoRegion = Era5PostProcessing.getGeoRegion(lonLayer, latLayer);
121-
final Rectangle era5RasterPosition = Era5PostProcessing.getEra5RasterPosition(geoRegion);
120+
final int[] shape = lonLayer.getShape();
121+
final int width = shape[1];
122+
final int height = shape[0];
123+
124+
//final Rectangle era5RasterPosition = Era5PostProcessing.getEra5RasterPosition(geoRegion);
122125
final InterpolationContext interpolationContext = Era5PostProcessing.getInterpolationContext(lonLayer, latLayer);
126+
final Rectangle layerRegion = interpolationContext.getEra5Region();
123127

124-
index.set(m);
125-
final int era5Time = era5TimeArray.getInt(index);
128+
timeIndex.set(m);
129+
final int era5Time = era5TimeArray.getInt(timeIndex);
126130

127131
// iterate over variables
128132
// + assemble variable file name
129-
// - read variable data extract
130-
// - interpolate (2d, 3d per layer)
133+
// + read variable data extract
134+
// + interpolate (2d, 3d per layer)
131135
// - store to target raster
132136
final Set<String> variableKeys = variables.keySet();
133137
for (final String variableKey : variableKeys) {
134138
final Variable variable = variableCache.get(variableKey, era5Time);
135-
final Array subset = readSubset(numLayers, era5RasterPosition, variableKey, variable);
136-
137-
// @todo 1 tb/tb continue here 2020-11-25
138-
final BilinearInterpolator bilinearInterpolator = interpolationContext.get(0, 0);
139+
final Array subset = readSubset(numLayers, layerRegion, variableKey, variable);
140+
final Index subsetIndex = subset.getIndex();
141+
142+
// final TemplateVariable templateVariable = variables.get(variableKey);
143+
// final Variable targetVariable = writer.findVariable(templateVariable.getName());
144+
// final Array targetArray = Array.factory(subset.getDataType(), subset.getShape());
145+
// final Index targetIndex = targetArray.getIndex();
146+
147+
final int rank = subset.getRank();
148+
if (rank == 2) {
149+
for (int y = 0; y < height; y++) {
150+
for (int x = 0; x < width; x++) {
151+
final BilinearInterpolator interpolator = interpolationContext.get(x, y);
152+
final int offsetX = interpolator.getXMin() - layerRegion.x;
153+
final int offsetY = interpolator.getYMin() - layerRegion.y;
154+
155+
subsetIndex.set(offsetY, offsetX);
156+
final float c00 = subset.getFloat(subsetIndex);
157+
158+
subsetIndex.set(offsetY, offsetX + 1);
159+
final float c10 = subset.getFloat(subsetIndex);
160+
161+
subsetIndex.set(offsetY + 1, offsetX);
162+
final float c01 = subset.getFloat(subsetIndex);
163+
164+
subsetIndex.set(offsetY + 1, offsetX + 1);
165+
final float c11 = subset.getFloat(subsetIndex);
166+
167+
final double interpolate = interpolator.interpolate(c00, c01, c10, c11);
168+
// targetIndex.set(y, x);
169+
// targetArray.setFloat(targetIndex, (float) interpolate);
170+
}
171+
}
172+
} else if (rank == 3) {
173+
for (int z = 0; z < numLayers; z++) {
174+
for (int y = 0; y < height; y++) {
175+
for (int x = 0; x < width; x++) {
176+
final BilinearInterpolator interpolator = interpolationContext.get(x, y);
177+
final int offsetX = interpolator.getXMin() - layerRegion.x;
178+
final int offsetY = interpolator.getYMin() - layerRegion.y;
179+
180+
subsetIndex.set(z, offsetY, offsetX);
181+
final float c00 = subset.getFloat(subsetIndex);
182+
183+
subsetIndex.set(z, offsetY, offsetX + 1);
184+
final float c10 = subset.getFloat(subsetIndex);
185+
186+
subsetIndex.set(z, offsetY + 1, offsetX);
187+
final float c01 = subset.getFloat(subsetIndex);
188+
189+
subsetIndex.set(z, offsetY + 1, offsetX + 1);
190+
final float c11 = subset.getFloat(subsetIndex);
191+
192+
final double interpolate = interpolator.interpolate(c00, c01, c10, c11);
193+
// targetIndex.set(z, y, x);
194+
// targetArray.setFloat(targetIndex, (float) interpolate);
195+
}
196+
}
197+
}
198+
} else {
199+
throw new IllegalStateException("Unexpected variable rank: " + rank + " " + variableKey);
200+
}
201+
// writer.write(targetVariable, targetArray);
139202
}
140203
}
141204
} finally {
@@ -150,17 +213,16 @@ private Array readSubset(int numLayers, Rectangle era5RasterPosition, String var
150213
final int[] origin = new int[]{0, era5RasterPosition.y, era5RasterPosition.x};
151214
final int[] shape = new int[]{1, era5RasterPosition.height, era5RasterPosition.width};
152215
subset = variable.read(origin, shape);
153-
subset = subset.reduce();
154-
subset = NetCDFUtils.scaleIfNecessary(variable, subset);
155216
} else if (rank == 4) {
156217
final int[] origin = new int[]{0, 0, era5RasterPosition.y, era5RasterPosition.x};
157218
final int[] shape = new int[]{1, numLayers, era5RasterPosition.height, era5RasterPosition.width};
158219
subset = variable.read(origin, shape);
159-
subset = subset.reduce();
160-
subset = NetCDFUtils.scaleIfNecessary(variable, subset);
161220
} else {
162221
throw new IOException("variable rank invalid: " + variableKey);
163222
}
223+
224+
subset = subset.reduce();
225+
subset = NetCDFUtils.scaleIfNecessary(variable, subset);
164226
return subset;
165227
}
166228

‎post-processing-tool/src/test/java/com/bc/fiduceo/post/plugin/era5/BilinearInterpolatorTest.java‎

Lines changed: 20 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -8,71 +8,79 @@ public class BilinearInterpolatorTest {
88

99
@Test
1010
public void testInterpolate_flat() {
11-
final BilinearInterpolator interpolator = new BilinearInterpolator(0.34, 0.67);
11+
final BilinearInterpolator interpolator = new BilinearInterpolator(0.34, 0.67, 56, 123);
1212

1313
final double interpolate = interpolator.interpolate(1.f, 1.f, 1.f, 1.f);
1414
assertEquals(1.0, interpolate, 1e-8);
1515
}
1616

1717
@Test
1818
public void testInterpolate_rise_left_right() {
19-
final BilinearInterpolator interpolator = new BilinearInterpolator(0.5, 0.5);
19+
final BilinearInterpolator interpolator = new BilinearInterpolator(0.5, 0.5, 57, 124);
2020

2121
final double interpolate = interpolator.interpolate(-1.f, 1.f, -1.f, 1.f);
2222
assertEquals(0.0, interpolate, 1e-8);
2323
}
2424

2525
@Test
2626
public void testInterpolate_rise_top_down() {
27-
final BilinearInterpolator interpolator = new BilinearInterpolator(0.5, 0.5);
27+
final BilinearInterpolator interpolator = new BilinearInterpolator(0.5, 0.5, 58, 125);
2828

2929
final double interpolate = interpolator.interpolate(-1.f, -1.f, 1.f, 1.f);
3030
assertEquals(0.0, interpolate, 1e-8);
3131
}
3232

3333
@Test
3434
public void testInterpolate_standard() {
35-
final BilinearInterpolator interpolator = new BilinearInterpolator(0.18, 0.74);
35+
final BilinearInterpolator interpolator = new BilinearInterpolator(0.18, 0.74, 55, 126);
3636

3737
final double interpolate = interpolator.interpolate(18f, 19f, 18.4f, 19.13f);
3838
assertEquals(18.44003565673828, interpolate, 1e-8);
3939
}
4040

4141
@Test
4242
public void testInterpolate_boundaries() {
43-
final BilinearInterpolator left = new BilinearInterpolator(0.0, 0.67);
43+
final BilinearInterpolator left = new BilinearInterpolator(0.0, 0.67, 60, 127);
4444
double interpolate = left.interpolate(10f, 11f, 12f, 13f);
4545
assertEquals(11.34, interpolate, 1e-8);
4646

47-
final BilinearInterpolator right = new BilinearInterpolator(1.0, 0.67);
47+
final BilinearInterpolator right = new BilinearInterpolator(1.0, 0.67, 60, 127);
4848
interpolate = right.interpolate(10f, 11f, 12f, 13f);
4949
assertEquals(12.34, interpolate, 1e-8);
5050

51-
final BilinearInterpolator top = new BilinearInterpolator(0.66, 0.0);
51+
final BilinearInterpolator top = new BilinearInterpolator(0.66, 0.0, 60, 127);
5252
interpolate = top.interpolate(10f, 11f, 12f, 13f);
5353
assertEquals(10.66, interpolate, 1e-8);
5454

55-
final BilinearInterpolator bottom = new BilinearInterpolator(0.66, 1.0);
55+
final BilinearInterpolator bottom = new BilinearInterpolator(0.66, 1.0, 60, 127);
5656
interpolate = bottom.interpolate(10f, 11f, 12f, 13f);
5757
assertEquals(12.66, interpolate, 1e-8);
5858
}
5959

6060
@Test
6161
public void testInterpolate_corners() {
62-
final BilinearInterpolator upperLeft = new BilinearInterpolator(0.0, 0.0);
62+
final BilinearInterpolator upperLeft = new BilinearInterpolator(0.0, 0.0, 61, 128);
6363
double interpolate = upperLeft.interpolate(10f, 11f, 12f, 13f);
6464
assertEquals(10.0, interpolate, 1e-8);
6565

66-
final BilinearInterpolator upperRight = new BilinearInterpolator(1.0, 0.0);
66+
final BilinearInterpolator upperRight = new BilinearInterpolator(1.0, 0.0, 61, 128);
6767
interpolate = upperRight.interpolate(10f, 11f, 12f, 13f);
6868
assertEquals(11.0, interpolate, 1e-8);
6969

70-
final BilinearInterpolator lowerLeft= new BilinearInterpolator(0.0, 1.0);
70+
final BilinearInterpolator lowerLeft= new BilinearInterpolator(0.0, 1.0, 61, 128);
7171
interpolate = lowerLeft.interpolate(10f, 11f, 12f, 13f);
7272
assertEquals(12.0, interpolate, 1e-8);
7373

74-
final BilinearInterpolator lowerRight = new BilinearInterpolator(1.0, 1.0);
74+
final BilinearInterpolator lowerRight = new BilinearInterpolator(1.0, 1.0, 61, 128);
7575
interpolate = lowerRight.interpolate(10f, 11f, 12f, 13f);
7676
assertEquals(13.0, interpolate, 1e-8);
7777
}
78+
79+
@Test
80+
public void testConstructAndGetCoordinates() {
81+
final BilinearInterpolator interpolator = new BilinearInterpolator(0.0, 0.0, 62, 129);
82+
83+
assertEquals(62, interpolator.getXMin());
84+
assertEquals(129, interpolator.getYMin());
85+
}
7886
}

0 commit comments

Comments
 (0)