001 /**
002 *
003 */
004 package org.wdssii.nse;
005
006 import java.util.ArrayList;
007 import java.util.Collections;
008 import java.util.Comparator;
009 import java.util.Date;
010 import java.util.List;
011
012 import org.apache.commons.logging.Log;
013 import org.apache.commons.logging.LogFactory;
014 import org.wdssii.core.Algorithm;
015 import org.wdssii.core.BuilderFactory;
016 import org.wdssii.core.DataEncoder;
017 import org.wdssii.core.DataType;
018 import org.wdssii.core.IndexRecord;
019 import org.wdssii.core.LatLonGrid;
020 import org.wdssii.core.Location;
021 import org.wdssii.core.UdUnits;
022
023 import ucar.units.Unit;
024
025 /**
026 *
027 * Near-storm environment algorithm.
028 *
029 * @author lakshman
030 *
031 */
032 public class NseAlgorithm extends Algorithm {
033 private static Log log = LogFactory.getLog(NseAlgorithm.class);
034
035 private String heightProductName = "HGT-ISBL";
036 private String temperatureProductName = "TMP-ISBL";
037 private String relativeHumidityProductName = "RH-ISBL";
038 private String gridRelativeUWindProductName = "UGRD-ISBL";
039 private String gridRelativeVWindProductName = "VGRD-ISBL";
040 private String verticalVelocityProductName = "VVEL-ISBL";
041 private String terrainElevationProductName = "HGT-SFC";
042 private String mslPressureProductName = "MSLMA-MSL";
043 private String sfcGridRelativeUWindProductName = "UGRD-HTGL";
044 private String sfcGridRelativeVWindProductName = "VGRD-HTGL";
045 private String sfcPressureProductName = "PRES-SFC";
046 private String sfcTemperatureProductName = "TMP-HTGL";
047 private String sfcDewPointDepression = "DEPR-HTGL";
048 private String subTypeSuffix = "_analysis";
049 private int numPressureLevels = 31;
050
051 public static void main(String[] args) throws Exception {
052 NseAlgorithm alg = new NseAlgorithm();
053 alg.setup(args);
054 alg.init();
055 alg.execute();
056 }
057
058 public static class InProduct {
059 private String productName;
060 private Unit unit;
061 private List<LatLonGrid> llgrids = new ArrayList<LatLonGrid>();
062
063 public InProduct(String productName, String unit) {
064 this.productName = productName;
065 if (unit != null) {
066 this.unit = UdUnits.parse(unit);
067 }
068 }
069
070 public void update(LatLonGrid grid) {
071 if (this.getTime().before(grid.getTime())) {
072 log.info("New " + productName
073 + " grid received ... deleting older grids");
074 llgrids.clear();
075 }
076 if (unit != null) {
077 grid.convertToUnit(unit);
078 }
079 llgrids.add(grid);
080 log.info("Received " + llgrids.size() + " products for "
081 + productName);
082 }
083
084 public int getNumGrids() {
085 return llgrids.size();
086 }
087
088 public Date getTime() {
089 if (llgrids.size() == 0) {
090 return new Date(0);
091 }
092 return llgrids.get(0).getTime();
093 }
094
095 public void sortByPressureLevel() {
096 Collections.sort(llgrids, new Comparator<LatLonGrid>() {
097 @Override
098 public int compare(LatLonGrid o1, LatLonGrid o2) {
099 Double p1 = new Double((String) o1
100 .getAttribute("PressureLevel"));
101 Double p2 = new Double((String) o2
102 .getAttribute("PressureLevel"));
103 // higher pressures first
104 return p2.compareTo(p1);
105 }
106 });
107 }
108
109 public float getPressureLevel(int index){
110 return Float.valueOf((String) llgrids.get(index).getAttribute("PressureLevel"));
111 }
112 }
113
114 private InProduct[] in3DProducts;
115 private InProduct[] inSfcProducts;
116
117 private void init() {
118 in3DProducts = new InProduct[] {
119 new InProduct(heightProductName, "gpm"),
120 new InProduct(temperatureProductName, "degreeC"),
121 new InProduct(relativeHumidityProductName, null), // will get changed to dewpoint
122 new InProduct(gridRelativeUWindProductName, "MetersPerSecond"), // will get changed to u
123 new InProduct(gridRelativeVWindProductName, "MetersPerSecond"), // will get changed to v
124 new InProduct(verticalVelocityProductName, "MetersPerSecond") };
125 inSfcProducts = new InProduct[] {
126 new InProduct(terrainElevationProductName, "m"),
127 new InProduct(mslPressureProductName, "mb"),
128 new InProduct(sfcGridRelativeUWindProductName, "MetersPerSecond"), // will get changed to u
129 new InProduct(sfcGridRelativeVWindProductName, "MetersPerSecond"), // will get changed to v
130 new InProduct(sfcPressureProductName, "mb"),
131 new InProduct(sfcTemperatureProductName, "degreeC"),
132 new InProduct(sfcDewPointDepression, "degreeC"), }; // will get changed to dewpoint
133 }
134
135 private void runNSE(){
136 // set up pressure levels
137 float[] pressureLevels = new float[numPressureLevels];
138 for (int i=0; i < numPressureLevels; ++i){
139 float ht = in3DProducts[0].getPressureLevel(i);
140 float temp = in3DProducts[1].getPressureLevel(i);
141 float dewp = in3DProducts[2].getPressureLevel(i);
142 float u = in3DProducts[3].getPressureLevel(i);
143 float v = in3DProducts[4].getPressureLevel(i);
144 float w = in3DProducts[5].getPressureLevel(i);
145 if (ht != temp || dewp != temp || u != temp || v != temp || w != temp){
146 log.error("The pressure levels in the grids do not match ...");
147 return;
148 }
149 pressureLevels[i] = ht;
150 }
151
152 // set up NSE profile
153 NSEProfile profile = new NSEProfile(pressureLevels);
154 profile.sfcElevation = inSfcProducts[0].llgrids.get(0).getValues();
155 profile.mslPressure = inSfcProducts[1].llgrids.get(0).getValues();
156 profile.sfcUWind = inSfcProducts[2].llgrids.get(0).getValues();
157 profile.sfcVWind = inSfcProducts[3].llgrids.get(0).getValues();
158 profile.sfcPressure = inSfcProducts[4].llgrids.get(0).getValues();
159 profile.sfcTemperature = inSfcProducts[5].llgrids.get(0).getValues();
160 profile.sfcDewPoint = inSfcProducts[6].llgrids.get(0).getValues();
161
162 for (int i=0; i < numPressureLevels; ++i){
163 profile.height[i] = in3DProducts[0].llgrids.get(i).getValues();
164 profile.temperatureC[i] = in3DProducts[1].llgrids.get(i).getValues();
165 profile.dewpoint[i] = in3DProducts[2].llgrids.get(i).getValues();
166 profile.u[i] = in3DProducts[3].llgrids.get(i).getValues();
167 profile.v[i] = in3DProducts[4].llgrids.get(i).getValues();
168 profile.w[i] = in3DProducts[5].llgrids.get(i).getValues();
169 }
170 profile.finishInit();
171
172 // write out isotherm products
173 LatLonGrid output = new LatLonGrid(inSfcProducts[0].llgrids.get(0));
174 float[][] temp2DArray = output.getValues();
175
176 profile.fillHeightOfIsotherm(0, temp2DArray);
177 writeOutputGrid(output, temp2DArray, "HeightOf0C", "nseanalysis", "m");
178
179 profile.fillHeightOfIsotherm(-10, temp2DArray);
180 writeOutputGrid(output, temp2DArray, "HeightOf-10C", "nseanalysis", "m");
181
182 profile.fillHeightOfIsotherm(-20, temp2DArray);
183 writeOutputGrid(output, temp2DArray, "HeightOf-20C", "nseanalysis", "m");
184
185 profile.fillHeightOfWetBulb(0, temp2DArray);
186 writeOutputGrid(output, temp2DArray, "HeightOfWetBulbZero", "nseanalysis", "m");
187
188 profile.fillHeightOfWetBulb(4, temp2DArray);
189 writeOutputGrid(output, temp2DArray, "HeightOfWetBulbFour", "nseanalysis", "m");
190
191 profile.fillMeanWindSpeed(500, 300, temp2DArray);
192 writeOutputGrid(output, temp2DArray, "WindSpeedMean500-300mb", "nseanalysis", "m/s");
193
194 profile.fillMeanValuePressureLayerVerticalVelocity(850, 500, temp2DArray);
195 writeOutputGrid(output, temp2DArray, "VerticalVelocityMean_850-500mb", "nseanalysis", "dimensionless");
196
197 profile.fillMeanValuePressureLayerTemperatureC(500, 400, temp2DArray);
198 writeOutputGrid(output, temp2DArray, "TemperatureMean_500-400mb", "nseanalysis", "degreeC");
199
200 profile.fillNumberOfUndergroundGridpoints(temp2DArray);
201 writeOutputGrid(output, temp2DArray, "NumUndergroundGridpoints", "nseanalysis", "dimensionless");
202
203 profile.fillSurfaceRelativeHumidityPercent(temp2DArray);
204 writeOutputGrid(output, temp2DArray, "SfcRH", "nseanalysis", "%");
205
206 writeOutputGrid(output, profile.sfcDewPoint, "SfcDewPoint", "nseanalysis", "degreeC");
207
208 writeOutputGrid(output, profile.sfcVirtualTemperature, "SfcTv", "nseanalysis", "degreeC");
209
210 writeOutputGrid(output, profile.sfcThetaE, "SfcThetaE", "nseanalysis", "degreeK");
211
212 writeOutputGrid(output, profile.sfcTheta, "SfcTheta", "nseanalysis", "degreeK");
213
214 writeOutputGrid(output, profile.sfcMixingRatio, "SfcMixingRatio", "nseanalysis", "GramsPerKilogram");
215
216 writeOutputGrid(output, profile.sfcLclTemperature, "SfcLCLTemp", "nseanalysis", "degreeC");
217
218 writeOutputGrid(output, profile.sfcLclPressure, "SfcLCLPress", "nseanalysis", "mb");
219
220 writeOutputGrid(output, profile.sfcWetBulbTemperatureK, "SfcWetBulb", "nseanalysis", "degreeK");
221 }
222
223 private void writeOutputGrid(LatLonGrid output, float[][] values, String typeName, String subType, String unit) {
224 float[][] orig = output.getValues();
225 output.setValues(values);
226 try {
227 output.setTypeName(typeName);
228 output.setAttribute("SubType", subType);
229 output.setAttribute("Unit", unit);
230 DataEncoder.writeDataAndNotify(output, getOutputDir());
231 } finally {
232 output.setValues(orig);
233 }
234 }
235
236 private boolean updateProduct(IndexRecord rec) {
237 String dataType = rec.getDataType();
238 for (InProduct prod : in3DProducts) {
239 if (prod.productName.equals(dataType)) {
240 LatLonGrid grid = (LatLonGrid) BuilderFactory
241 .createDataType(rec);
242 prod.update(grid);
243 return true;
244 }
245 }
246 for (InProduct prod : inSfcProducts) {
247 if (prod.productName.equals(dataType)) {
248 LatLonGrid grid = (LatLonGrid) BuilderFactory
249 .createDataType(rec);
250 prod.update(grid);
251 return true;
252 }
253 }
254 log.debug("Not interested in " + dataType);
255 return false;
256 }
257
258 private boolean suffixMatches(IndexRecord rec) {
259 String subType = rec.getSubType();
260 if (!subType.endsWith(subTypeSuffix)) {
261 if (log.isDebugEnabled()) {
262 String dataType = rec.getDataType();
263 log.debug("Ignoring " + dataType + ":" + subType
264 + " since the subType does not end with "
265 + subTypeSuffix);
266 }
267 return false;
268 }
269 return true;
270 }
271
272 private boolean haveAllProducts() {
273 // do we have all the grids we need?
274 for (InProduct prod : in3DProducts) {
275 if (prod.llgrids.size() != numPressureLevels) {
276 log.debug("Have received only " + prod.llgrids.size() + " of "
277 + numPressureLevels + " pressure levels for "
278 + prod.productName);
279 return false;
280 }
281 }
282 for (InProduct prod : inSfcProducts) {
283 if (prod.llgrids.size() != 1) {
284 log.debug("Have not yet received " + prod.productName);
285 return false;
286 }
287 }
288 log.info("Have received all the products required.");
289
290 return true;
291 }
292
293 private void clearAllProducts() {
294 for (InProduct prod : in3DProducts) {
295 prod.llgrids.clear();
296 }
297 for (InProduct prod : inSfcProducts) {
298 prod.llgrids.clear();
299 }
300 }
301
302 private void prepareWindProducts(List<LatLonGrid> ugrid,
303 List<LatLonGrid> vgrid) {
304 log.info("Converting grid-relative winds in "
305 + ugrid.get(0).getTypeName() + " north-relative");
306 log.info("Converting grid-relative winds in "
307 + vgrid.get(0).getTypeName() + " north-relative");
308
309 // Here are some constants for the Lambert Conformal proejection:
310 final double WindRotationConstant = 0.422618;
311 final double XMeridian = -95.0;
312
313 Location nwcorner = ugrid.get(0).getLocation();
314 double w_lat = nwcorner.getLongitude();
315 float lonres = ugrid.get(0).getDeltaLon();
316 int numLat = ugrid.get(0).getNumLat();
317 int numLon = ugrid.get(0).getNumLon();
318 for (int k = 0; k < ugrid.size(); ++k) {
319 float[][] u = ugrid.get(k).getValues();
320 float[][] v = vgrid.get(k).getValues();
321 for (int i = 0; i < numLat; ++i) {
322 for (int j = 0; j < numLon; ++j) {
323 if (u[i][j] == DataType.MissingData
324 || v[i][j] == DataType.MissingData) {
325 u[i][j] = v[i][j] = DataType.MissingData;
326 } else {
327 double ang = WindRotationConstant
328 * (w_lat + j * lonres - XMeridian) * 0.017453;
329 double sinx = Math.sin(ang);
330 double cosx = Math.cos(ang);
331 double ug = u[i][j];
332 double vg = v[i][j];
333 double u_new = cosx * ug + sinx * vg;
334 double v_new = -sinx * ug + cosx * vg;
335 // don't need a decimal precision more than 2 places
336 int u_trunc = (int) (u_new * 100);
337 int v_trunc = (int) (v_new * 100);
338 u[i][j] = u_trunc * 0.01f;
339 v[i][j] = v_trunc * 0.01f;
340 }
341
342 }
343 }
344 ugrid.get(k).setTypeName("U");
345 vgrid.get(k).setTypeName("V");
346 }
347 }
348
349 public static float calcDewpointFromRH(float temp, float rh) {
350 float dewp = DataType.MissingData;
351
352 if ((temp != DataType.MissingData) && (rh != DataType.MissingData)) {
353 float RH = rh * 0.01f;
354 // from Bolton (MWR 1980):
355 float es = 6.112f * (float) Math.exp((17.67f * temp) / (temp + 243.5f));
356 float e = es * RH;
357 dewp = (243.5f * (float) Math.log(e) - 440.8f) / (19.48f - (float) Math.log(e));
358
359 }
360
361 return dewp;
362
363 }
364
365 private void prepareRelativeHumidity(List<LatLonGrid> tempGrids,
366 List<LatLonGrid> rhGrids) {
367 log.info("Converting relative humidity in "
368 + rhGrids.get(0).getTypeName()
369 + " to dewpoint based on temperature from "
370 + tempGrids.get(0).getTypeName());
371 int numLat = rhGrids.get(0).getNumLat();
372 int numLon = rhGrids.get(0).getNumLon();
373 for (int k = 0; k < rhGrids.size(); ++k) {
374 float[][] rhvalues = rhGrids.get(k).getValues();
375 float[][] tempvalues = tempGrids.get(k).getValues();
376 for (int i = 0; i < numLat; ++i) {
377 for (int j = 0; j < numLon; ++j) {
378 float rh = rhvalues[i][j];
379 float temp = tempvalues[i][j];
380 rhvalues[i][j] = calcDewpointFromRH(temp, rh);
381 }
382 }
383 rhGrids.get(k).setTypeName("DewPoint");
384 }
385 }
386
387 private void prepareDewPointDepression(List<LatLonGrid> tempGrids,
388 List<LatLonGrid> ddGrids) {
389 log.info("Converting dewpoint depression in "
390 + ddGrids.get(0).getTypeName()
391 + " to dewpoint based on temperature from "
392 + tempGrids.get(0).getTypeName());
393 int numLat = ddGrids.get(0).getNumLat();
394 int numLon = ddGrids.get(0).getNumLon();
395 for (int k = 0; k < ddGrids.size(); ++k) {
396 float[][] ddvalues = ddGrids.get(k).getValues();
397 float[][] tempvalues = tempGrids.get(k).getValues();
398 for (int i = 0; i < numLat; ++i) {
399 for (int j = 0; j < numLon; ++j) {
400 float dd = ddvalues[i][j];
401 float temp = tempvalues[i][j];
402 ddvalues[i][j] = (dd != DataType.MissingData && temp != DataType.MissingData)? DataType.MissingData : (temp-dd);
403 }
404 }
405 ddGrids.get(k).setTypeName("DewPoint");
406 }
407 }
408
409 @Override
410 public void handleRecord(IndexRecord rec) {
411 if (suffixMatches(rec) && updateProduct(rec) && haveAllProducts()) {
412 try {
413 for (InProduct prod : in3DProducts) {
414 prod.sortByPressureLevel();
415 }
416 prepareWindProducts(in3DProducts[3].llgrids,
417 in3DProducts[4].llgrids);
418 prepareWindProducts(inSfcProducts[2].llgrids,
419 inSfcProducts[3].llgrids);
420 prepareRelativeHumidity(in3DProducts[1].llgrids,
421 in3DProducts[2].llgrids);
422 prepareDewPointDepression(inSfcProducts[5].llgrids,
423 inSfcProducts[6].llgrids);
424
425 runNSE();
426
427 } finally {
428 clearAllProducts();
429 }
430 }
431 }
432 }