Source code

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