001 /**
002 *
003 */
004 package org.wdssii.singleradar;
005
006 import org.apache.commons.logging.Log;
007 import org.apache.commons.logging.LogFactory;
008 import org.wdssii.core.Algorithm;
009 import org.wdssii.core.BuilderFactory;
010 import org.wdssii.core.DataEncoder;
011 import org.wdssii.core.DataType;
012 import org.wdssii.core.IndexRecord;
013 import org.wdssii.core.PolarGrid;
014 import org.wdssii.core.RadialSet;
015 import org.wdssii.core.VCP;
016
017 /**
018 * @author lakshman
019 *
020 */
021 public class VilAlgorithm extends Algorithm {
022 private static Log log = LogFactory.getLog(VilAlgorithm.class);
023
024 private String inputType = "Reflectivity";
025
026 private boolean vilEnabled = true;
027
028 private boolean compositeEnabled = true;
029
030 private boolean rapidUpdateEnabled = true;
031
032 private float azimuthalSpacingDegrees = 1.0f;
033
034 private float gateWidthKms = 1.0f;
035
036 private int numGates = 500;
037
038 private VCP vcp = null;
039
040 private PolarGrid[] elevations = null;
041
042 private int latest = -1;
043
044 public String getInputType() {
045 return inputType;
046 }
047
048 public void setInputType(String inputType) {
049 this.inputType = inputType;
050 }
051
052 public boolean isCompositeEnabled() {
053 return compositeEnabled;
054 }
055
056 public void setCompositeEnabled(boolean compositeEnabled) {
057 this.compositeEnabled = compositeEnabled;
058 }
059
060 public boolean isRapidUpdateEnabled() {
061 return rapidUpdateEnabled;
062 }
063
064 public void setRapidUpdateEnabled(boolean rapidUpdateEnabled) {
065 this.rapidUpdateEnabled = rapidUpdateEnabled;
066 }
067
068 public boolean isVilEnabled() {
069 return vilEnabled;
070 }
071
072 public void setVilEnabled(boolean vilEnabled) {
073 this.vilEnabled = vilEnabled;
074 }
075
076 public float getAzimuthalSpacingDegrees() {
077 return azimuthalSpacingDegrees;
078 }
079
080 public void setAzimuthalSpacingDegrees(float azimuthalSpacing) {
081 this.azimuthalSpacingDegrees = (azimuthalSpacing);
082 }
083
084 public float getGateWidthKms() {
085 return gateWidthKms;
086 }
087
088 public void setGateWidthKms(float gateWidth) {
089 this.gateWidthKms = (gateWidth);
090 }
091
092 public int getNumGates() {
093 return numGates;
094 }
095
096 public void setNumGates(int numGates) {
097 this.numGates = (numGates);
098 }
099
100 /**
101 * Workhorse method of algorithm; will be called by Index on every new
102 * record.
103 */
104 public synchronized void handleRecord(IndexRecord rec) {
105 try {
106 if (log.isDebugEnabled()) {
107 log.debug("Received: " + rec);
108 }
109 if (!rec.getDataType().equals(inputType)) {
110 return;
111 }
112 RadialSet rs = (RadialSet) BuilderFactory.createDataType(rec);
113 PolarGrid pg = new PolarGrid(rs, getAzimuthalSpacingDegrees(),
114 getGateWidthKms()*1000, getNumGates(), true);
115 update(pg);
116 if (isReadyToCompute()) {
117 String[] subtypes = getSubTypes();
118 PolarGrid composite = computeComposite();
119 if (compositeEnabled) {
120 DataEncoder.writeDataAndNotify(composite, getOutputDir(), subtypes);
121 }
122 if (vilEnabled) {
123 PolarGrid vil = computeVil(composite);
124 DataEncoder.writeDataAndNotify(vil, getOutputDir(), subtypes);
125 }
126 }
127 } catch (Exception e) {
128 log.error("Logic error", e);
129 // continue
130 }
131 }
132
133 private PolarGrid computeVil(PolarGrid composite) {
134 PolarGrid example = composite;
135 PolarGrid vil = new PolarGrid(example,"VIL", 0, 0);
136 vil.addAttributes(example);
137 vil.setAttribute("Unit", "kg/m^2");
138 int numRadials = vil.getNumRadials();
139 int numGates = vil.getNumGates();
140 float[][] comp = composite.getValues();
141 float[][] result = vil.getValues();
142 float[] midbeam_sines = computeMidbeamSines();
143 for (int gate = 0; gate < numGates; gate++) {
144 for (int radial = 0; radial < numRadials; radial++) {
145 if (comp[radial][gate] <= 0) {
146 continue;
147 }
148 float temporary1 = 0;
149 for (int elev = 0; elev < elevations.length; elev++) {
150 // compute Z
151 float dbz = elevations[elev].getValues()[radial][gate];
152 if (dbz == DataType.MissingData) {
153 continue;
154 }
155 float range = vil.getRangeToFirstGateKms() + (gate + 0.5f) * vil.getGateWidthKms();
156 // depth of beam
157 float hAgl;
158 if (elev == 0) {
159 hAgl = heightAboveGround(range, midbeam_sines[elev]);
160 } else {
161 float hAgl_this = heightAboveGround(range,
162 midbeam_sines[elev]);
163 float hAgl_below = heightAboveGround(range,
164 midbeam_sines[elev - 1]);
165 hAgl = hAgl_this - hAgl_below;
166 }
167 temporary1 += getPowPow(dbz) * hAgl;
168 }// elev
169 if (temporary1 <= 1)
170 result[radial][gate] = DataType.MissingData;
171 else
172 result[radial][gate] = temporary1;
173 }// gate
174 }// radial
175
176 return vil;
177 }
178
179 private static final float minDBZ = -30;
180
181 private static final float maxDBZ = 55;
182
183 private static final float deltaDBZ = 0.5f;
184
185 private static final float[] lookup = makePowPowLookup();
186
187 private static float[] makePowPowLookup() {
188 int n = Math.round((maxDBZ - minDBZ) / deltaDBZ) + 1;
189 float[] lookup = new float[n];
190 final double puissance = 4.0 / 7.0;
191 for (int i = 0; i < n; ++i) {
192 float dbz = minDBZ + i * deltaDBZ;
193 double Z = Math.pow(10, (dbz / 10));
194 lookup[i] = (float) (0.00000344 * Math.pow(Z, puissance));
195 }
196 return lookup;
197 }
198
199 private float getPowPow(float dbz) {
200 int bin = Math.round((dbz - minDBZ) / deltaDBZ);
201 if (bin < 0)
202 bin = 0;
203 else if (bin >= lookup.length)
204 bin = lookup.length - 1;
205 return lookup[bin];
206 }
207
208 private static final float FACTOR = (float) (2.0 * 1.21 * 6371.0);
209
210 private float heightAboveGround(float range, float sine_angle) {
211 float zhgt;
212 float hgt;// , val;
213 hgt = (range * sine_angle) + ((range * range) / FACTOR);
214 zhgt = hgt * 1000;
215 return (zhgt);
216 }
217
218 private float[] computeMidbeamSines() {
219 float[] midbeam_sines = new float[elevations.length];
220 for (int elev = 0; elev < elevations.length; elev++) {
221 float this_angle = elevations[elev].getElevation() + 0.5f
222 * elevations[elev].getAngularResolution(); // current
223 midbeam_sines[elev] = (float) Math.sin(Math.toRadians(this_angle));
224 }
225 return midbeam_sines;
226 }
227
228 private String[] getSubTypes() {
229 float elev = elevations[latest].getElevation();
230 String subType = getSubTypeForElevation(elev);
231 String[] subtypes;
232 if (isVolumeComplete()) {
233 subtypes = new String[] { subType, "vol" };
234 } else {
235 subtypes = new String[] { subType };
236 }
237 return subtypes;
238 }
239
240 private void update(PolarGrid pg) {
241 // get vcp
242 int newvcp = Integer.parseInt((String) pg.getAttribute("vcp"));
243 if (vcp == null || newvcp != vcp.getVcp()) {
244 if (log.isInfoEnabled()) {
245 log.info("VCP changed to " + newvcp
246 + " resetting elevation history");
247 }
248 elevations = null;
249 vcp = VCP.getVCP(newvcp);
250 elevations = new PolarGrid[vcp.getLastTilt()];
251 }
252
253 // find appropriate elevation
254 int tiltNo = vcp.getTiltNumber(pg.getElevation());
255 latest = tiltNo - 1;
256 elevations[latest] = pg;
257 if (log.isInfoEnabled()) {
258 log.info("updated with " + latest + "th elevation at "
259 + pg.getElevation());
260 }
261 }
262
263 private boolean isVolumeComplete() {
264 return (elevations != null && latest >= 0 && latest == (elevations.length - 1));
265 }
266
267 private boolean isReadyToCompute() {
268 if (elevations == null || elevations.length == 0) {
269 log.debug("No elevations yet");
270 return false;
271 }
272 for (int i = 0; i < elevations.length; ++i) {
273 if (elevations[i] == null) {
274 if (log.isDebugEnabled()) {
275 log.debug(i + "th elevation not received");
276 }
277 return false;
278 }
279 }
280 if (!rapidUpdateEnabled && !isVolumeComplete()) {
281 log.debug("The latest elevation received is not last in volume");
282 return false;
283 }
284 log.debug("Ready to compute");
285 return true;
286 }
287
288 private PolarGrid computeComposite() {
289 PolarGrid example = elevations[latest];
290 PolarGrid composite = new PolarGrid(example, example.getTypeName(), 0, 0);
291 composite.addAttributes(example);
292 int numRadials = composite.getNumRadials();
293 int numGates = composite.getNumGates();
294 float[][] result = composite.getValues();
295 for (int e = 0; e < elevations.length; ++e) {
296 float[][] values = elevations[e].getValues();
297 for (int i = 0; i < numRadials; ++i) {
298 for (int j = 0; j < numGates; ++j) {
299 if (values[i][j] > result[i][j]) {
300 result[i][j] = values[i][j];
301 }
302 }
303 }
304 }
305 log.debug("composite computed");
306 return composite;
307 }
308
309 public static void main(String[] args) throws Exception {
310 VilAlgorithm alg = new VilAlgorithm();
311 alg.setupAndExecute(args);
312 }
313
314 }