001 /**
002 *
003 */
004 package org.wdssii.dualpol.icing;
005
006 import java.io.File;
007 import java.io.FileWriter;
008 import java.io.IOException;
009 import java.io.PrintWriter;
010 import java.util.ArrayList;
011 import java.util.Date;
012 import java.util.List;
013
014 import org.apache.commons.logging.Log;
015 import org.apache.commons.logging.LogFactory;
016 import org.wdssii.core.DataType;
017 import org.wdssii.core.LatLonGrid;
018 import org.wdssii.core.NetcdfDataEncoder;
019 import org.wdssii.core.WDSSII;
020 import org.wdssii.core.LatLonGrid.Pixel;
021 import org.wdssii.dualpol.IDLNetcdfGridReader;
022
023
024 /**
025 * @author lakshman
026 *
027 */
028 public class Statistics {
029 private static Log log = LogFactory.getLog(Statistics.class);
030 public static void main(String[] args) throws IOException {
031 WDSSII.getInstance(); // sets timezone etc.
032 if ( args.length < 3 ){
033 System.err.println("Usage: java " + Statistics.class.getCanonicalName() + " in_pireps_dir in_icing_dir output_dir");
034 return;
035 }
036
037 Statistics stats = new Statistics();
038 stats.writeOutputs(args[0], args[1], args[2]);
039 }
040
041 public void writeOutputs(String pireps_dir, String icing_dir, String output_dir) throws IOException {
042 PilotReport[] pireps = new PirepsReader().readReports(pireps_dir, "150km");
043 IcingFile[] icings = IcingFile.getAllFiles(icing_dir, "KOUN");
044 new File(output_dir).mkdirs();
045
046 final long maxTimeDifferenceMsecs = 30 * 60 * 60 * 1000; // 30 minutes
047 int num_written = 0;
048 for (PilotReport r : pireps){
049 IcingFile nearest_icing = getClosestInTime(icings, r.getTime(), maxTimeDifferenceMsecs );
050 if (log.isInfoEnabled()){
051 log.info(r + "->" + nearest_icing);
052 }
053 if ( nearest_icing != null ){
054 String outName = "stat_" + NetcdfDataEncoder.getFormattedDate(r.getTime()) + "_" + num_written + ".txt";
055 File f = new File(output_dir, outName);
056 IcingGrids icingGrids = new IcingGrids( new IDLNetcdfGridReader().fromNetcdfFile(nearest_icing.file) );
057 writeStatsFile( r, icingGrids , f );
058 ++num_written;
059 }
060 }
061 }
062
063
064 private void writeStatsFile(PilotReport report, IcingGrids grid, File outFile ) throws IOException {
065 StringBuilder comment = new StringBuilder("#");
066 StringBuilder data = new StringBuilder();
067 String sep = "\t";
068
069 comment.append("Time").append(sep);
070 data.append( NetcdfDataEncoder.getFormattedDate(report.getTime())).append(sep);
071
072 comment.append("Lat").append(sep);
073 data.append( report.getLocation().getLatitude() ).append(sep);
074
075 comment.append("Lon").append(sep);
076 data.append( report.getLocation().getLongitude() ).append(sep);
077
078 comment.append("Ht").append(sep);
079 data.append( report.getLocation().getHeightKms() ).append(sep);
080
081 comment.append("Icing").append(sep);
082 data.append( report.getIcingIntensity() ).append(sep);
083
084 for (int i=0; i < IcingGrids.gridTypes.length; ++i){
085 LatLonGrid[] grids = grid.getGrids(i);
086 // get closest grid-point in space
087 LatLonGrid curr= getClosestGrid( report.getLocation().getHeightKms(), grids, 1.0 );
088 if ( curr == null ){
089 log.warn("Skipping " + report + " as it is above/below target area");
090 return;
091 }
092 LatLonGrid.Pixel pixel = curr.getPixel(report.getLocation());
093 if ( pixel == null ){
094 log.warn("Skipping " + report + " as it is outside target area");
095 return;
096 }
097
098 // grids above & below this grid point
099 LatLonGrid above = getAboveGrid( report.getLocation().getHeightKms(), grids );
100 LatLonGrid below = getBelowGrid( report.getLocation().getHeightKms(), grids );
101 if ( above == null ){
102 above = curr;
103 }
104 if ( below == null ){
105 below = curr;
106 }
107
108 // stats
109 addTextureStatistics(comment, data, sep, pixel, curr);
110 addVerticalDifferenceStatistics(comment, data, sep, pixel, curr, above, below);
111 addVerticalColumnStatistics(comment, data, sep, pixel, grids);
112 }
113
114 // now write the stuff out
115 PrintWriter writer = new PrintWriter(new FileWriter(outFile));
116 writer.println(comment.toString());
117 writer.println(data.toString());
118 writer.close();
119 }
120
121 private void addTextureStatistics(StringBuilder comment,
122 StringBuilder data, String sep, Pixel pixel, LatLonGrid curr) {
123 float max = 0;
124 float min = 999;
125 float sumx = 0;
126 float sumx2 = 0;
127 int N = 0;
128
129 for (int i = pixel.x - 2; i <= pixel.x + 2; ++i) for (int j = pixel.y - 2; j <= pixel.y + 2; ++j){
130 try {
131 float value = curr.getValues()[i][j];
132 if ( value == DataType.MissingData ){
133 continue;
134 }
135 if ( value > max ){
136 value = max;
137 }
138 if ( value < min ){
139 value = min;
140 }
141 sumx += value;
142 sumx2 += (value*value);
143 N ++;
144 } catch (Exception e){
145 // continue
146 }
147 }
148
149 float mean = (N > 0)? (sumx / N) : DataType.MissingData;
150 float var = (N > 1)? ((sumx2 - (sumx*sumx)/N)/(N-1)) : 0;
151 if ( N == 0 ){
152 min = mean = max = DataType.MissingData;
153 }
154 String typeName = curr.getTypeName();
155 comment.append("hor_avg" + typeName).append(sep);
156 data.append( mean ).append(sep);
157 comment.append("hor_max" + typeName).append(sep);
158 data.append( max ).append(sep);
159 comment.append("hor_min" + typeName).append(sep);
160 data.append( min ).append(sep);
161 comment.append("hor_var" + typeName).append(sep);
162 data.append( var ).append(sep);
163 comment.append("hor_N" + typeName).append(sep);
164 data.append( N ).append(sep);
165 }
166
167 private void addVerticalColumnStatistics(StringBuilder comment,
168 StringBuilder data, String sep, LatLonGrid.Pixel loc,
169 LatLonGrid[] grids) {
170 float max = 0;
171 float min = 999;
172 float sumx = 0;
173 float sumx2 = 0;
174 int N = 0;
175
176 for (int g=0; g < grids.length; ++g){
177 float value = grids[g].getValue(loc);
178 if ( value == DataType.MissingData ){
179 continue;
180 }
181 if ( value > max ){
182 value = max;
183 }
184 if ( value < min ){
185 value = min;
186 }
187 sumx += value;
188 sumx2 += (value*value);
189 N ++;
190 }
191
192 float mean = (N > 0)? (sumx / N) : DataType.MissingData;
193 float var = (N > 1)? ((sumx2 - (sumx*sumx)/N)/(N-1)) : 0;
194 if ( N == 0 ){
195 min = mean = max = DataType.MissingData;
196 }
197 String typeName = grids[0].getTypeName();
198 comment.append("vert_avg" + typeName).append(sep);
199 data.append( mean ).append(sep);
200 comment.append("vert_max" + typeName).append(sep);
201 data.append( max ).append(sep);
202 comment.append("vert_min" + typeName).append(sep);
203 data.append( min ).append(sep);
204 comment.append("vert_var" + typeName).append(sep);
205 data.append( var ).append(sep);
206 comment.append("vert_N" + typeName).append(sep);
207 data.append( N ).append(sep);
208 }
209
210 private static float diff(float a, float b){
211 if ( a == DataType.MissingData ){
212 a = 0;
213 }
214 if ( b == DataType.MissingData ){
215 b = 0;
216 }
217 return (a - b);
218 }
219
220 private void addVerticalDifferenceStatistics(StringBuilder comment, StringBuilder data, String sep,
221 LatLonGrid.Pixel loc, LatLonGrid curr, LatLonGrid above, LatLonGrid below) {
222
223 // data
224 comment.append(curr.getTypeName()).append(sep);
225 float value = curr.getValue(loc);
226 data.append( value ).append(sep);
227
228 // deltaabove
229 comment.append("d" + curr.getTypeName() + "_up").append(sep);
230 float aboveValue = above.getValue(loc);
231 data.append( diff(aboveValue,value) ).append(sep);
232
233 // deltabelow
234 comment.append("d" + curr.getTypeName() + "_down").append(sep);
235 float belowValue = above.getValue(loc);
236 data.append( diff(value,belowValue) ).append(sep);
237
238 }
239
240 private LatLonGrid getClosestGrid(double ht, LatLonGrid[] grids, double max_diff){
241 LatLonGrid best = null;
242 double best_diff = max_diff;
243 for (LatLonGrid grid : grids){
244 double diff = Math.abs(grid.getLocation().getHeightKms() - ht);
245 if ( diff < best_diff ){
246 best = grid;
247 best_diff = diff;
248 }
249 }
250 return best;
251 }
252
253 private LatLonGrid getAboveGrid(double ht, LatLonGrid[] grids){
254 LatLonGrid best = null;
255 double best_diff = Double.MAX_VALUE;
256 for (LatLonGrid grid : grids){
257 double diff = (grid.getLocation().getHeightKms() - ht);
258 if ( diff > 0.001 && diff < best_diff ){
259 best = grid;
260 best_diff = diff;
261 }
262 }
263 return best;
264 }
265
266 private LatLonGrid getBelowGrid(double ht, LatLonGrid[] grids){
267 LatLonGrid best = null;
268 double best_diff = Double.MAX_VALUE;
269 for (LatLonGrid grid : grids){
270 double diff = (ht - grid.getLocation().getHeightKms());
271 if ( diff > 0.001 && diff < best_diff ){
272 best = grid;
273 best_diff = diff;
274 }
275 }
276 return best;
277 }
278
279 private IcingFile getClosestInTime(IcingFile[] candidates, Date time, long maxTimeDiffMsecs ){
280 if ( candidates.length == 0 ){
281 return null;
282 }
283
284 // find closest in time
285 IcingFile best = candidates[0];
286 long bestTimeDiff = Math.abs(best.time.getTime() - time.getTime());
287 for (IcingFile f : candidates){
288 long timeDiff = Math.abs(f.time.getTime() - time.getTime());
289 if ( timeDiff < bestTimeDiff ){
290 bestTimeDiff = timeDiff;
291 best = f;
292 }
293 }
294
295 // check against threshold
296 if ( bestTimeDiff < maxTimeDiffMsecs ){
297 return best;
298 } else {
299 return null;
300 }
301 }
302
303 public static class IcingFile {
304 public final File file;
305 public final Date time;
306 public IcingFile(File f, Date time) {
307 this.file = f;
308 this.time = time;
309 }
310 public static IcingFile[] getAllFiles(String dir, String pattern){
311 File[] icing_files = new File(dir).listFiles();
312 List<IcingFile> result = new ArrayList<IcingFile>();
313 for (File f : icing_files){
314 if (f.getName().contains(pattern)){
315 IcingFile file = new IcingFile(f, IDLNetcdfGridReader.getTimeFromFileName(f.getName()));
316 result.add(file);
317 }
318 }
319 return result.toArray(new IcingFile[0]);
320 }
321 @Override
322 public String toString(){
323 return file.getName() + " at " + time;
324 }
325 }
326 }