001/*
002 * $Id: ProfileAlongTrack.java,v 1.10 2011/03/24 16:06:33 davep Exp $
003 *
004 * This file is part of McIDAS-V
005 *
006 * Copyright 2007-2011
007 * Space Science and Engineering Center (SSEC)
008 * University of Wisconsin - Madison
009 * 1225 W. Dayton Street, Madison, WI 53706, USA
010 * https://www.ssec.wisc.edu/mcidas
011 * 
012 * All Rights Reserved
013 * 
014 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and
015 * some McIDAS-V source code is based on IDV and VisAD source code.  
016 * 
017 * McIDAS-V is free software; you can redistribute it and/or modify
018 * it under the terms of the GNU Lesser Public License as published by
019 * the Free Software Foundation; either version 3 of the License, or
020 * (at your option) any later version.
021 * 
022 * McIDAS-V is distributed in the hope that it will be useful,
023 * but WITHOUT ANY WARRANTY; without even the implied warranty of
024 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
025 * GNU Lesser Public License for more details.
026 * 
027 * You should have received a copy of the GNU Lesser Public License
028 * along with this program.  If not, see http://www.gnu.org/licenses.
029 */
030
031package edu.wisc.ssec.mcidasv.data.hydra;
032
033import visad.Data;
034import visad.FlatField;
035import visad.Set;
036import visad.Gridded1DSet;
037import visad.CoordinateSystem;
038import visad.RealType;
039import visad.RealTupleType;
040import visad.SetType;
041import visad.Linear1DSet;
042import visad.Linear2DSet;
043import visad.Unit;
044import visad.FunctionType;
045import visad.VisADException;
046import visad.QuickSort;
047import java.rmi.RemoteException;
048
049import java.util.Hashtable;
050import java.util.HashMap;
051import java.util.StringTokenizer;
052
053import java.io.BufferedReader;
054import java.io.FileInputStream;
055import java.io.IOException;
056import java.io.InputStream;
057import java.io.InputStreamReader;
058
059
060public abstract class ProfileAlongTrack extends MultiDimensionAdapter {
061
062      private FunctionType mathtype;
063
064      int TrackLen;
065      int VertLen;
066
067      private float[] vertLocs = null;
068      private float[] trackTimes = null;
069      private float[] trackLongitude = null;
070      private float[] trackLatitude = null;
071
072      public static String longitude_name = "Longitude";
073      public static String latitude_name  = "Latitude";
074      public static String trackDim_name  = "TrackDim";
075      public static String vertDim_name  = "VertDim";
076      public static String array_name = "array_name";
077      public static String profileTime_name = "ProfileTime";
078      public static String profileTime_unit = "ProfileTime_unit";
079      public static String altitude_unit = "altitude_unit";
080      public static String sfcElev_name = "SurfaceElev";
081      public static String range_name = "range_name";
082      public static String scale_name = "scale_name";
083      public static String offset_name = "offset_name";
084      public static String fill_value_name = "fill_value_name";
085      public static String valid_range = "valid_range";
086      public static String ancillary_file_name = "ancillary_file";
087      static String product_name = "product_name";
088      
089      String[] rangeName_s  = null;
090      Class[] arrayType_s = null;
091      Unit[] rangeUnit_s  = new Unit[] {null};
092
093      RealType track  = RealType.getRealType(trackDim_name);
094      RealType vert = RealType.getRealType(vertDim_name);
095      RealType[] domainRealTypes = new RealType[2];
096
097      RealType vertLocType;
098      RealType trackTimeType;
099
100      int track_idx      = -1;
101      int vert_idx       = -1;
102      int range_rank     = -1;
103
104      int track_tup_idx;
105      int vert_tup_idx;
106
107      //-private RangeProcessor rangeProcessor;
108
109      CoordinateSystem cs = null;
110
111      public static HashMap getEmptySubset() {
112        HashMap<String, double[]> subset = new HashMap<String, double[]>();
113        subset.put(trackDim_name, new double[3]);
114        subset.put(vertDim_name, new double[3]);
115        return subset;
116      }
117
118      public static HashMap getEmptyMetadataTable() {
119        HashMap<String, String> metadata = new HashMap<String, String>();
120        metadata.put(array_name, null);
121        metadata.put(trackDim_name, null);
122        metadata.put(vertDim_name, null);
123        metadata.put(longitude_name, null);
124        metadata.put(latitude_name, null);
125        metadata.put(profileTime_name, null);
126        metadata.put(profileTime_unit, null);
127        metadata.put(altitude_unit, null);
128        metadata.put(sfcElev_name, null);
129        metadata.put(scale_name, null);
130        metadata.put(offset_name, null);
131        metadata.put(fill_value_name, null);
132        /*
133        metadata.put(range_name, null);
134        metadata.put(range_unit, null);
135        metadata.put(valid_range, null);
136        */
137        return metadata;
138      }
139
140      public ProfileAlongTrack() {
141
142      }
143
144      public ProfileAlongTrack(MultiDimensionReader reader, HashMap metadata) {
145        super(reader, metadata);
146        this.init();
147      }
148
149
150      private void init() {
151        for (int k=0; k<array_rank;k++) {
152          if ( ((String)metadata.get(trackDim_name)).equals(array_dim_names[k]) ) {
153            track_idx = k;
154          }
155          if ( ((String)metadata.get(vertDim_name)).equals(array_dim_names[k]) ) {
156            vert_idx = k;
157          }
158        }
159
160        int[] lengths = new int[2];
161
162        if (track_idx < vert_idx) {
163          domainRealTypes[0] = vert;
164          domainRealTypes[1] = track;
165          track_tup_idx = 1;
166          vert_tup_idx = 0;
167          lengths[0] = array_dim_lengths[vert_idx];
168          lengths[1] = array_dim_lengths[track_idx];
169        }
170        else {
171          domainRealTypes[0] = track;
172          domainRealTypes[1] = vert;
173          track_tup_idx = 0;
174          vert_tup_idx = 1;
175          lengths[0] = array_dim_lengths[track_idx];
176          lengths[1] = array_dim_lengths[vert_idx];
177        }
178
179        TrackLen = array_dim_lengths[track_idx];
180        VertLen = array_dim_lengths[vert_idx];
181
182        String rangeName = null;
183        if (metadata.get("range_name") != null) {
184          rangeName = (String)metadata.get("range_name");
185        } 
186        else {
187          rangeName = (String)metadata.get("array_name");
188        }
189        rangeType = RealType.getRealType(rangeName, rangeUnit_s[0]);
190
191        try {
192          rangeProcessor = RangeProcessor.createRangeProcessor(reader, metadata);
193        } 
194        catch (Exception e) {
195          System.out.println("RangeProcessor failed to create");
196          e.printStackTrace();
197        }
198
199        try {
200          vertLocs = getVertBinAltitude();
201          trackTimes = getTrackTimes();
202          vertLocType = makeVertLocType();
203          trackTimeType = makeTrackTimeType();
204          trackLongitude = getTrackLongitude();
205          trackLatitude = getTrackLatitude();
206        } 
207        catch (Exception e) {
208          System.out.println(e);
209        }
210
211      }
212
213      public int getTrackLength() {
214        return TrackLen;
215      }
216
217      public int getVertLength() {
218        return VertLen;
219      }
220
221      public int getVertIdx() {
222        return vert_idx;
223      }
224
225      public int getTrackIdx() {
226        return track_idx;
227      }
228
229      public int getVertTupIdx() {
230        return vert_tup_idx;
231      }
232
233      public int getTrackTupIdx() {
234        return track_tup_idx;
235      }
236                                                                                                                                                     
237      public Set makeDomain(Object subset) throws Exception {
238        double[] first = new double[2];
239        double[] last = new double[2];
240        int[] length = new int[2];
241
242        HashMap<String, double[]> domainSubset = new HashMap<String, double[]>();
243        domainSubset.put(trackDim_name, (double[]) ((HashMap)subset).get(trackDim_name));
244        domainSubset.put(vertDim_name, (double[]) ((HashMap)subset).get(vertDim_name));
245
246        for (int kk=0; kk<2; kk++) {
247          RealType rtype = domainRealTypes[kk];
248          double[] coords = (double[]) ((HashMap)subset).get(rtype.getName());
249          first[kk] = coords[0];
250          last[kk] = coords[1];
251          length[kk] = (int) ((last[kk] - first[kk])/coords[2] + 1);
252          last[kk] = first[kk]+coords[2]*(length[kk]-1);
253        }
254        Linear2DSet domainSet = new Linear2DSet(first[0], last[0], length[0], first[1], last[1], length[1]);
255        final Linear1DSet[] lin1DSet_s = new Linear1DSet[] {domainSet.getLinear1DComponent(0),
256                                           domainSet.getLinear1DComponent(1)};
257
258        float[] new_altitudes = new float[length[vert_tup_idx]];
259        float[] new_times = new float[length[track_tup_idx]];
260
261        int track_start = (int) first[track_tup_idx];
262        int vert_start = (int) first[vert_tup_idx];
263        int vert_skip = (int) ((double[]) ((HashMap)subset).get(vertDim_name))[2];
264        int track_skip = (int) ((double[]) ((HashMap)subset).get(trackDim_name))[2];
265        for (int k=0; k<new_altitudes.length; k++) {
266          new_altitudes[k] = vertLocs[vert_start+(k*vert_skip)];
267        }
268
269        for (int k=0; k<new_times.length; k++) {
270          new_times[k] = trackTimes[track_start+(k*track_skip)];
271        }
272
273        final Gridded1DSet alt_set = new Gridded1DSet(vertLocType, new float[][] {new_altitudes}, new_altitudes.length);
274        final Gridded1DSet time_set = new Gridded1DSet(trackTimeType, new float[][] {new_times}, new_times.length);
275        final float vert_offset = (float) first[vert_tup_idx];
276        final float track_offset = (float) first[track_tup_idx];
277
278        RealTupleType reference = new RealTupleType(vertLocType, trackTimeType);
279        
280        CoordinateSystem cs = null;
281
282        try {
283        cs = new CoordinateSystem(reference, null) {
284          public float[][] toReference(float[][] vals) throws VisADException {
285            int[] indexes = lin1DSet_s[0].valueToIndex(new float[][] {vals[0]});
286            for (int k=0; k<vals[0].length;k++) {
287              //-indexes[k] = (int) (vals[vert_tup_idx][k] - vert_offset); ?
288            }
289            float[][] alts = alt_set.indexToValue(indexes);
290
291            indexes = lin1DSet_s[1].valueToIndex(new float[][] {vals[1]});
292            for (int k=0; k<vals[0].length;k++) {
293              //-indexes[k] = (int) (vals[track_tup_idx][k] - track_offset); ?
294            }
295            float[][] times = time_set.indexToValue(indexes);
296
297            return new float[][] {alts[0], times[0]};
298          }
299          public float[][] fromReference(float[][] vals) throws VisADException {
300            int[] indexes = alt_set.valueToIndex(new float[][] {vals[vert_tup_idx]});
301            float[][] vert_coords = lin1DSet_s[vert_tup_idx].indexToValue(indexes);
302            indexes = time_set.valueToIndex(new float[][] {vals[track_tup_idx]});
303            float[][] track_coords = lin1DSet_s[track_tup_idx].indexToValue(indexes);
304            return new float[][] {vert_coords[0], track_coords[0]};
305          }
306          public double[][] toReference(double[][] vals) throws VisADException {
307            return Set.floatToDouble(toReference(Set.doubleToFloat(vals)));
308          }
309          public double[][] fromReference(double[][] vals) throws VisADException {
310            return Set.floatToDouble(fromReference(Set.doubleToFloat(vals)));
311          }
312          public boolean equals(Object obj) {
313            return true;
314          }
315        };
316        }
317        catch (Exception e) {
318        }
319
320        RealTupleType domainTupType = new RealTupleType(domainRealTypes[0], domainRealTypes[1], cs, null);
321        domainSet = new Linear2DSet(domainTupType, first[0], last[0], length[0], first[1], last[1], length[1]);
322
323        return domainSet;
324      }
325
326      public FunctionType getMathType() {
327        return null;
328      }
329
330      public RealType[] getDomainRealTypes() {
331        return domainRealTypes;
332      }
333
334      public HashMap getDefaultSubset() {
335        HashMap subset = ProfileAlongTrack.getEmptySubset();
336
337        double[] coords = (double[])subset.get("TrackDim");
338        coords[0] = 20000.0;
339        coords[1] = (TrackLen - 15000.0) - 1;
340        //-coords[2] = 30.0;
341        coords[2] = 5.0;
342        subset.put("TrackDim", coords);
343
344        coords = (double[])subset.get("VertDim");
345        coords[0] = 98.0;
346        coords[1] = (VertLen) - 1;
347        coords[2] = 2.0;
348        subset.put("VertDim", coords);
349        return subset;
350      }
351
352      public int[] getTrackRangeInsideLonLatRect(double minLat, double maxLat, double minLon, double maxLon) {
353        int nn = 100;
354        int skip = TrackLen/nn;
355        double lon;
356        double lat;
357        
358        int idx = 0;
359        while (idx < TrackLen) {
360          lon = (double)trackLongitude[idx];
361          lat = (double)trackLatitude[idx];
362          if (((lon > minLon) && (lon < maxLon)) && ((lat > minLat)&&(lat < maxLat))) break;
363          idx += skip;
364        }
365        if (idx > TrackLen-1) idx = TrackLen-1;
366        if (idx == TrackLen-1) return new int[] {-1,-1};
367
368        int low_idx = idx;
369        while (low_idx > 0) {
370          lon = (double)trackLongitude[low_idx];
371          lat = (double)trackLatitude[low_idx];
372          if (((lon > minLon) && (lon < maxLon)) && ((lat > minLat)&&(lat < maxLat))) {
373            low_idx -= 1;
374            continue;
375          }
376          else {
377            break;
378          }
379        }
380
381        int hi_idx = idx;
382        while (hi_idx < TrackLen-1) {
383          lon = (double)trackLongitude[hi_idx];
384          lat = (double)trackLatitude[hi_idx];
385          if (((lon > minLon) && (lon < maxLon)) && ((lat > minLat)&&(lat < maxLat))) {
386            hi_idx += 1;
387            continue;
388          }
389          else {
390            break;
391          }
392        }
393        return new int[] {low_idx, hi_idx};
394      }
395
396      public HashMap getSubsetFromLonLatRect(HashMap subset, double minLat, double maxLat, double minLon, double maxLon) {
397        double[] coords = (double[])subset.get("TrackDim");
398        int[] idxs = getTrackRangeInsideLonLatRect(minLat, maxLat, minLon, maxLon);
399        coords[0] = (double) idxs[0];
400        coords[1] = (double) idxs[1];
401        if ((coords[0] == -1) || (coords[1] == -1)) return null;
402        return subset;
403      }
404
405      public HashMap getSubsetFromLonLatRect(HashMap subset, double minLat, double maxLat, double minLon, double maxLon,
406                                             int xStride, int yStride, int zStride) {
407
408        double[] coords = (double[])subset.get("TrackDim"); 
409        int[] idxs = getTrackRangeInsideLonLatRect(minLat, maxLat, minLon, maxLon);
410        coords[0] = (double) idxs[0];
411        coords[1] = (double) idxs[1];
412        if ((coords[0] == -1) || (coords[1] == -1)) return null;
413        if (xStride > 0) {
414          coords[2] = xStride;
415        }
416
417        coords = (double[])subset.get("VertDim");
418        if (yStride > 0) {
419          coords[2] = yStride;
420        }
421        return subset;
422      }
423
424      public HashMap getSubsetFromLonLatRect(double minLat, double maxLat, double minLon, double maxLon) {
425        return getSubsetFromLonLatRect(getDefaultSubset(), minLat, maxLat, minLon, maxLon);
426      }
427
428      public HashMap getSubsetFromLonLatRect(double minLat, double maxLat, double minLon, double maxLon,
429                                             int xStride, int yStride, int zStride) {
430
431        return getSubsetFromLonLatRect(getDefaultSubset(), minLat, maxLat, minLon, maxLon, xStride, yStride, zStride);
432      }
433
434      public abstract float[] getVertBinAltitude() throws Exception;
435      public abstract float[] getTrackTimes() throws Exception;
436      public abstract RealType makeVertLocType() throws Exception;
437      public abstract RealType makeTrackTimeType() throws Exception;
438      public abstract float[] getTrackLongitude() throws Exception;
439      public abstract float[] getTrackLatitude() throws Exception;
440
441      public static float[] medianFilter(float[] A, int lenx, int leny, int window_lenx, int window_leny)
442           throws VisADException {
443        float[] result =  new float[A.length];
444        float[] window =  new float[window_lenx*window_leny];
445        float[] new_window =  new float[window_lenx*window_leny];
446        int[] sort_indexes = new int[window_lenx*window_leny];
447                                                                                                                                    
448        int a_idx;
449        int w_idx;
450                                                                                                                                    
451        int w_lenx = window_lenx/2;
452        int w_leny = window_leny/2;
453                                                                                                                                    
454        int lo;
455        int hi;
456        int ww_jj;
457        int ww_ii;
458        int cnt;
459                                                                                                                                    
460        for (int j=0; j<leny; j++) {
461          for (int i=0; i<lenx; i++) {
462            a_idx = j*lenx + i;
463            cnt = 0;
464            for (int w_j=-w_leny; w_j<w_leny; w_j++) {
465              for (int w_i=-w_lenx; w_i<w_lenx; w_i++) {
466                ww_jj = w_j + j;
467                ww_ii = w_i + i;
468                w_idx = (w_j+w_leny)*window_lenx + (w_i+w_lenx);
469                if ((ww_jj >= 0) && (ww_ii >=0) && (ww_jj < leny) && (ww_ii < lenx)) {
470                  window[cnt] = A[ww_jj*lenx+ww_ii];
471                  cnt++;
472                }
473              }
474            }
475            System.arraycopy(window, 0, new_window, 0, cnt);
476            //-sort_indexes = QuickSort.sort(new_window, sort_indexes);
477            sort_indexes = QuickSort.sort(new_window);
478            result[a_idx] = new_window[cnt/2];
479          }
480        }
481        return result;
482      }
483
484}