Source code

001/*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2024
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 https://www.gnu.org/licenses/.
027 */
028
029package edu.wisc.ssec.mcidasv.data.hydra;
030
031import visad.CoordinateSystem;
032import visad.VisADException;
033import visad.RealTupleType;
034import visad.RealType;
035import visad.Linear2DSet;
036import visad.Gridded2DSet;
037import visad.Gridded1DSet;
038import visad.Set;
039import java.awt.geom.Rectangle2D;
040
041import visad.data.hdfeos.LambertAzimuthalEqualArea;
042import visad.Data;
043
044
045public class LongitudeLatitudeCoordinateSystem extends CoordinateSystem {
046
047   private static final long serialVersionUID = 1L;
048
049   Linear2DSet domainSet;
050   Linear2DSet subSet;
051   Gridded2DSet gset;
052   Gridded2DSet projSet;
053   CoordinateSystem projCS;
054   float[][] lonlat;
055
056   Gridded1DSet goodLinesSet;
057
058   //- assumes incoming GriddedSet is (longitude,latitude) with range (-180,+180)
059   boolean neg180pos180 = true;  //false: longitude range (0,+360)
060
061   public LongitudeLatitudeCoordinateSystem(Linear2DSet domainSet, Gridded2DSet gset) throws VisADException {
062     this(domainSet, gset, false);
063   }
064
065   public LongitudeLatitudeCoordinateSystem(Linear2DSet domainSet, Gridded2DSet gset, boolean lonFlag) throws VisADException {
066     super(RealTupleType.SpatialEarth2DTuple, null);
067     this.gset = gset;
068     this.domainSet = domainSet;
069     this.neg180pos180 = lonFlag;
070     int[] lengths = domainSet.getLengths();
071     int[] gset_lengths = gset.getLengths();
072     subSet = new Linear2DSet(0.0, gset_lengths[0]-1, lengths[0],
073                              0.0, gset_lengths[1]-1, lengths[1]);
074
075     lonlat = gset.getSamples(false);
076
077     // Check for missing (complete) scan lines.
078     int[] goodLines = checkForMissingLines(lonlat, gset_lengths[0], gset_lengths[1]);
079     int lenY = gset_lengths[1];
080     if (goodLines != null) {
081        lonlat = removeMissingLines(lonlat, gset_lengths[0], gset_lengths[1], goodLines);
082        lenY = goodLines.length;
083        float[] flts = new float[lenY];
084        for (int k=0; k<lenY; k++) {
085            flts[k] = (float) goodLines[k];
086        }
087        goodLinesSet = new Gridded1DSet(RealType.Generic, new float[][] {flts}, lenY);
088     }
089
090     // Create map projected set of locations to handle problems using Lon/Lat
091     // coordinates around the dateline, prime-meridian or poles.
092     // Use missing removed lonlat array to get center lon,lat
093     int gxc = gset_lengths[0]/2;
094     int gyc = gset_lengths[1]/2;
095     int iic = gyc*gset_lengths[0] + gxc;
096     float[][] cntr = new float[][] {{lonlat[0][iic]}, {lonlat[1][iic]}};
097
098     float earthRadius = 6367470;
099     float lonCenter = cntr[0][0];
100     float latCenter = cntr[1][0];
101     projCS = new LambertAzimuthalEqualArea(RealTupleType.SpatialEarth2DTuple, earthRadius,
102                     lonCenter*Data.DEGREES_TO_RADIANS, latCenter*Data.DEGREES_TO_RADIANS,
103                          0,0);
104
105     float[][] grdLocs = projCS.fromReference(lonlat);
106
107     projSet = new Gridded2DSet(RealTupleType.SpatialCartesian2DTuple,
108                           grdLocs, gset_lengths[0], lenY,
109                                               null, null, null, false, false);
110   }
111
112   public float[][] toReference(float[][] values) throws VisADException {
113     float[][] coords = domainSet.valueToGrid(values);
114     coords = subSet.gridToValue(coords);
115
116     // Use nearest neighbor if grid box straddles dateline because of resulting
117     // interpolation problems.
118
119     int[] lens = gset.getLengths();
120     int lenX = lens[0];
121     float[][] lonlat = gset.getSamples(false);
122
123     float min =  Float.MAX_VALUE;
124     float max = -Float.MAX_VALUE;
125     float lon;
126     int gx, gy, idx;
127
128     for (int i=0; i<coords[0].length; i++) {
129        gx = (int) coords[0][i];
130        gy = (int) coords[1][i];
131        idx = gy*lenX + gx;
132
133        min =  Float.MAX_VALUE;
134        max = -Float.MAX_VALUE;
135
136        // look at grid cell corners
137        lon = lonlat[0][idx];
138        if (lon < min) min = lon;
139        if (lon > max) max = lon;
140
141        if ((gx+1) < lens[0]) {
142           lon = lonlat[0][idx+1];
143           if (lon < min) min = lon;
144           if (lon > max) max = lon;
145        }
146
147        if ((gy+1) < lens[1]) {
148           lon = lonlat[0][idx + lenX];
149           if (lon < min) min = lon;
150           if (lon > max) max = lon;
151        }
152
153        if (((gx+1) < lens[0]) && ((gy+1) < lens[1])) {
154           lon = lonlat[0][idx + lenX + 1];
155           if (lon < min) min = lon;
156           if (lon > max) max = lon;
157        }
158
159        if ((max - min) > 300) { // grid cell probably straddles the dateline so force nearest neighbor
160           coords[0][i] = (float) Math.floor(coords[0][i] + 0.5);
161           coords[1][i] = (float) Math.floor(coords[1][i] + 0.5);
162        }
163     }
164
165     coords = gset.gridToValue(coords);
166     // original set of lon,lat may contain fill values so perform a valid lat range check
167     for (int k=0; k<coords[0].length; k++) {
168        if (Math.abs(coords[1][k]) > 90) {
169           coords[0][k] = Float.NaN;
170           coords[1][k] = Float.NaN;
171        }
172     }
173     return coords;
174   }
175
176   public float[][] fromReference(float[][] values) throws VisADException {
177     if (!neg180pos180) { // force to longitude range (0,360)
178       for (int t=0; t<values[0].length; t++) {
179         if (values[0][t] > 180f) {
180           values[0][t] -= 360f;
181         }
182       }
183     }
184
185     //float[][] grid_vals = gset.valueToGrid(values);
186     // use the projected set
187     values = projCS.fromReference(values);
188     float[][] grid_vals = projSet.valueToGrid(values);
189
190     // return original domain coordinates if missing geo lines were removed
191     if (goodLinesSet != null) {
192        float[][] tmp = goodLinesSet.gridToValue(new float[][] {grid_vals[1]});
193        grid_vals[1] = tmp[0];
194     }
195
196     float[][] coords = subSet.valueToGrid(grid_vals);
197     coords = domainSet.gridToValue(coords);
198     return coords;
199   }
200
201   public double[][] toReference(double[][] values) throws VisADException {
202     float[][] coords = domainSet.valueToGrid(Set.doubleToFloat(values));
203     coords = subSet.gridToValue(coords);
204
205     // Use nearest neighbor if grid box straddles dateline because of resulting
206     // interpolation problems.
207
208     int[] lens = gset.getLengths();
209     int lenX = lens[0];
210
211     float min =  Float.MAX_VALUE;
212     float max = -Float.MAX_VALUE;
213     float lon;
214     int gx, gy, idx;
215
216     for (int i=0; i<coords[0].length; i++) {
217        gx = (int) coords[0][i];
218        gy = (int) coords[1][i];
219        idx = gy*lenX + gx;
220
221        min =  Float.MAX_VALUE;
222        max = -Float.MAX_VALUE;
223
224        lon = lonlat[0][idx];
225        if (lon < min) min = lon;
226        if (lon > max) max = lon;
227
228        if ((gx+1) < lens[0]) {
229           lon = lonlat[0][idx+1];
230           if (lon < min) min = lon;
231           if (lon > max) max = lon;
232        }
233
234        if ((gy+1) < lens[1]) {
235           lon = lonlat[0][idx + lenX];
236           if (lon < min) min = lon;
237           if (lon > max) max = lon;
238        }
239
240        if (((gx+1) < lens[0]) && ((gy+1) < lens[1])) {
241           lon = lonlat[0][idx + lenX + 1];
242           if (lon < min) min = lon;
243           if (lon > max) max = lon;
244        }
245
246        if ((max - min) > 300) { // grid cell probably straddles the dateline so force nearest neighbor
247           coords[0][i] = (float) Math.floor(coords[0][i] + 0.5);
248           coords[1][i] = (float) Math.floor(coords[1][i] + 0.5);
249        }
250     }
251
252     coords = gset.gridToValue(coords);
253     // original set of lon,lat may contain fill values so perform a valid lat range check
254     for (int k=0; k<coords[0].length; k++) {
255        if (Math.abs(coords[1][k]) > 90) {
256           coords[0][k] = Float.NaN;
257           coords[1][k] = Float.NaN;
258        }
259     }
260     return Set.floatToDouble(coords);
261   }
262
263   public double[][] fromReference(double[][] values) throws VisADException {
264     if (!neg180pos180) { // force to longitude range (0,360)
265       for (int t=0; t<values[0].length; t++) {
266         if (values[0][t] > 180.0) {
267           values[0][t] -= 360.0;
268         }
269       }
270     }
271
272     // use the projected set
273     float[][] grid_vals = projSet.valueToGrid(Set.doubleToFloat(values));
274
275     // return original domain coordinates if missing geo lines were removed
276     if (goodLinesSet != null) {
277        float[][] tmp = goodLinesSet.gridToValue(new float[][] {grid_vals[1]});
278        grid_vals[1] = tmp[0];
279     }
280
281     float[][] coords = subSet.valueToGrid(grid_vals);
282     coords = domainSet.gridToValue(coords);
283     return Set.floatToDouble(coords);
284   }
285
286   public Rectangle2D getDefaultMapArea() {
287     float[] lo = domainSet.getLow();
288     float[] hi = domainSet.getHi();
289     return new Rectangle2D.Float(lo[0], lo[1], hi[0] - lo[0], hi[1] - lo[1]);
290   }
291
292   public boolean equals(Object obj) {
293     if (obj instanceof LongitudeLatitudeCoordinateSystem) {
294        LongitudeLatitudeCoordinateSystem llcs = (LongitudeLatitudeCoordinateSystem) obj;
295        if (domainSet.equals(llcs.domainSet) && (gset.equals(llcs.gset))) {
296           return true;
297        }
298     }
299     return false;
300   }
301
302   public int[] checkForMissingLines(float[][] lonlat, int lenX, int lenY) {
303      int[] good_lines = new int[lenY];
304      int cnt = 0;
305      for (int j=0; j<lenY; j++) {
306         int idx = j*lenX + lenX/2; // check scan line center: NaN and valid Lat range check
307         if ((!(Float.isNaN(lonlat[0][idx]) || Float.isNaN(lonlat[1][idx]))) && (Math.abs(lonlat[1][idx]) <= 90.0)) {
308            good_lines[cnt++] = j;
309         }
310      }
311
312      if (cnt == lenY) {
313         return null;
314      }
315      else {
316         int[] tmp = new int[cnt];
317         System.arraycopy(good_lines, 0, tmp, 0, cnt);
318         good_lines = tmp;
319         return good_lines;
320      }
321   }
322
323   public float[][] removeMissingLines(float[][] lonlat, int lenX, int lenY, int[] good_lines) {
324      float[][] noMissing = new float[2][lenX*(good_lines.length)];
325
326      for (int k=0; k < good_lines.length; k++) {
327
328         int idx = good_lines[k]*lenX;
329
330         System.arraycopy(lonlat[0], idx, noMissing[0], k*lenX, lenX);
331         System.arraycopy(lonlat[1], idx, noMissing[1], k*lenX, lenX);
332      }
333
334      return noMissing;
335   }
336
337   public Gridded2DSet getTheGridded2DSet() {
338      return gset;
339   }
340}