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.control;
030
031import visad.georef.MapProjection;
032import visad.data.hdfeos.LambertAzimuthalEqualArea;
033import visad.RealTupleType;
034import visad.CoordinateSystem;
035import visad.Data;
036import visad.SI;
037import visad.Unit;
038import java.awt.geom.Rectangle2D;
039import visad.VisADException;
040import java.rmi.RemoteException;
041
042
043public class LambertAEA extends MapProjection {
044
045   CoordinateSystem cs;
046   Rectangle2D rect;
047   float earthRadius = 6367470; //- meters
048
049   public LambertAEA(float[][] corners, float lonCenter, float latCenter) throws VisADException {
050      super(RealTupleType.SpatialEarth2DTuple, new Unit[] {SI.meter, SI.meter});
051
052      cs = new LambertAzimuthalEqualArea(RealTupleType.SpatialEarth2DTuple, earthRadius,
053                   lonCenter*Data.DEGREES_TO_RADIANS, latCenter*Data.DEGREES_TO_RADIANS,
054                         0,0);
055
056     float[][] xy = cs.fromReference(corners);
057
058     float min_x = Float.MAX_VALUE;
059     float min_y = Float.MAX_VALUE;
060     float max_x = -Float.MAX_VALUE;
061     float max_y = -Float.MAX_VALUE;
062
063     for (int k=0; k<xy[0].length;k++) {
064       if (xy[0][k] < min_x) min_x = xy[0][k];
065       if (xy[1][k] < min_y) min_y = xy[1][k];
066       if (xy[0][k] > max_x) max_x = xy[0][k];
067       if (xy[1][k] > max_y) max_y = xy[1][k];
068     }
069
070     float del_x = max_x - min_x;
071     float del_y = max_y - min_y;
072
073     boolean forceSquareMapArea = true;
074     if (forceSquareMapArea) {
075       if (del_x < del_y) {
076          del_x = del_y;
077       }
078       else if (del_y < del_x) {
079          del_y = del_x;
080       }
081     }
082
083     rect = new Rectangle2D.Float(-del_x/2, -del_y/2, del_x, del_y);
084   }
085
086   public LambertAEA(float[][] corners) throws VisADException {
087     super(RealTupleType.SpatialEarth2DTuple, new Unit[] {SI.meter, SI.meter});
088
089     boolean spanGM = false;
090     float lonA = corners[0][0];
091     float lonB = corners[0][1];
092     float lonC = corners[0][2];
093     float lonD = corners[0][3];
094     float latA = corners[1][0];
095     float latB = corners[1][1];
096     float latC = corners[1][2];
097     float latD = corners[1][3];
098
099     float diffAD = lonA - lonD;
100     //float diffBC = lonB - lonC;
101
102     if (Math.abs(diffAD) > 180) spanGM = true;
103     //if (Math.abs(diffBC) > 180) spanGM = true;
104
105     float lonCenter;
106
107     if (spanGM) {
108        float[] vals = MultiSpectralControl.minmax(new float[] {lonA, lonD});
109        float wLon = vals[1];
110        float eLon = vals[0];
111        float del = 360f - wLon + eLon;
112        lonCenter = wLon + del/2;
113        if (lonCenter > 360) lonCenter -= 360f;
114     }
115     else {
116        float[] vals = MultiSpectralControl.minmax(new float[] {lonA, lonD});
117        float minLon = vals[0];
118        float maxLon = vals[1];
119        lonCenter = minLon + (maxLon - minLon)/2;
120     }
121
122     float[] vals = MultiSpectralControl.minmax(corners[1]);
123     float minLat = vals[0];
124     float maxLat = vals[1];
125     float latCenter = minLat + (maxLat - minLat)/2;
126
127     cs = new LambertAzimuthalEqualArea(RealTupleType.SpatialEarth2DTuple, earthRadius,
128                   lonCenter*Data.DEGREES_TO_RADIANS, latCenter*Data.DEGREES_TO_RADIANS,
129                         0,0);
130
131     float[][] xy = cs.fromReference(corners);
132
133     float min_x = Float.MAX_VALUE;
134     float min_y = Float.MAX_VALUE;
135     float max_x = Float.MIN_VALUE;
136     float max_y = Float.MIN_VALUE;
137
138     for (int k=0; k<xy[0].length;k++) {
139       if (xy[0][k] < min_x) min_x = xy[0][k];
140       if (xy[1][k] < min_y) min_y = xy[1][k];
141       if (xy[0][k] > max_x) max_x = xy[0][k];
142       if (xy[1][k] > max_y) max_y = xy[1][k];
143     }
144
145     float del_x = max_x - min_x;
146     float del_y = max_y - min_y;
147
148     boolean forceSquareMapArea = true;
149     if (forceSquareMapArea) {
150       if (del_x < del_y) {
151          del_x = del_y;
152       }
153       else if (del_y < del_x) {
154          del_y = del_x;
155       }
156     }
157
158     min_x = -del_x/2;
159     min_y = -del_y/2;
160
161     rect = new Rectangle2D.Float(min_x, min_y, del_x, del_y);
162   }
163
164   public LambertAEA(Rectangle2D ll_rect) throws VisADException {
165     this(ll_rect, true);
166   }
167
168   public LambertAEA(Rectangle2D ll_rect, boolean forceSquareMapArea) throws VisADException {
169     super(RealTupleType.SpatialEarth2DTuple, new Unit[] {SI.meter, SI.meter});
170
171     float minLon = (float) ll_rect.getX();
172     float minLat = (float) ll_rect.getY();
173     float del_lon = (float) ll_rect.getWidth();
174     float del_lat = (float) ll_rect.getHeight();
175     float maxLon = minLon + del_lon;
176     float maxLat = minLat + del_lat;
177
178     float lonDiff = maxLon - minLon;
179     float lonCenter = minLon + (maxLon - minLon)/2;
180     if (lonDiff > 180f) {
181       lonCenter += 180f;
182     }
183     float latCenter = minLat + (maxLat - minLat)/2;
184
185     cs = new LambertAzimuthalEqualArea(RealTupleType.SpatialEarth2DTuple, earthRadius,
186                   lonCenter*Data.DEGREES_TO_RADIANS, latCenter*Data.DEGREES_TO_RADIANS,
187                         0,0);
188
189     float[][] xy = cs.fromReference(new float[][] {{minLon,maxLon,minLon,maxLon}, 
190                                                    {minLat,minLat,maxLat,maxLat}});
191
192
193     float min_x = Float.MAX_VALUE;
194     float min_y = Float.MAX_VALUE;
195     float max_x = Float.MIN_VALUE;
196     float max_y = Float.MIN_VALUE;
197
198     for (int k=0; k<xy[0].length;k++) {
199       if (xy[0][k] < min_x) min_x = xy[0][k];
200       if (xy[1][k] < min_y) min_y = xy[1][k];
201       if (xy[0][k] > max_x) max_x = xy[0][k];
202       if (xy[1][k] > max_y) max_y = xy[1][k];
203     }
204
205     float del_x = max_x - min_x;
206     float del_y = max_y - min_y;
207 
208     if (forceSquareMapArea) {
209       if (del_x < del_y) {
210         del_x = del_y;
211       }
212       else if (del_y < del_x) {
213         del_y = del_x;
214       }
215     }
216
217     min_x = -del_x/2;
218     min_y = -del_y/2;
219  
220     rect = new Rectangle2D.Float(min_x, min_y, del_x, del_y);
221   }
222
223   public Rectangle2D getDefaultMapArea() {
224     return rect;
225   }
226     
227   public float[][] toReference(float[][] values) throws VisADException {
228     return cs.toReference(values);
229   }
230
231   public float[][] fromReference(float[][] values) throws VisADException {
232     return cs.fromReference(values);
233   }
234
235   public double[][] toReference(double[][] values) throws VisADException {
236     return cs.toReference(values);
237   }
238
239   public double[][] fromReference(double[][] values) throws VisADException {
240     return cs.fromReference(values);
241   }
242
243   public boolean equals(Object cs) {
244     if ( cs instanceof LambertAEA ) {
245        LambertAEA that = (LambertAEA) cs;
246        if ( (this.cs.equals(that.cs)) && this.getDefaultMapArea().equals(that.getDefaultMapArea())) {
247           return true;
248        }
249     }
250     return false;
251   }
252
253}