Source code

001/*
002 * $Id: MultiSpectralData.java,v 1.20 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 java.awt.geom.Rectangle2D;
034import java.rmi.RemoteException;
035import java.util.ArrayList;
036import java.util.HashMap;
037import java.util.Iterator;
038
039import visad.CoordinateSystem;
040import visad.FlatField;
041import visad.FunctionType;
042import visad.Gridded2DSet;
043import visad.Linear1DSet;
044import visad.Linear2DSet;
045import visad.Real;
046import visad.RealTuple;
047import visad.RealTupleType;
048import visad.RealType;
049import visad.SampledSet;
050import visad.Set;
051import visad.SetType;
052import visad.VisADException;
053
054public class MultiSpectralData extends MultiDimensionAdapter {
055
056  SwathAdapter swathAdapter = null;
057  SpectrumAdapter spectrumAdapter = null;
058  CoordinateSystem cs = null;
059
060  HashMap spectrumSelect = null;
061  HashMap swathSelect = null;
062
063  String sensorName = null;
064  String platformName = null;
065  String paramName = null;
066  String inputParamName = null;
067  String name = null;
068
069  public float init_wavenumber = 919.50f;
070  public String init_bandName = null;
071
072  float[] dataRange = new float[] {180f, 320f};
073
074  boolean hasBandNames = false;
075  ArrayList<String> bandNameList = null;
076  HashMap<String, Float> bandNameMap = null;
077
078  
079  public MultiSpectralData(SwathAdapter swathAdapter, SpectrumAdapter spectrumAdapter,
080                           String inputParamName, String paramName, String sensorName, String platformName) {
081    this.swathAdapter = swathAdapter;
082    this.spectrumAdapter = spectrumAdapter;
083    this.paramName = paramName;
084    this.inputParamName = inputParamName;
085    this.name = swathAdapter.getArrayName();
086
087    if (spectrumAdapter != null) {
088      this.spectrumSelect = spectrumAdapter.getDefaultSubset();
089      if (spectrumAdapter.hasBandNames()) {
090        hasBandNames = true;
091        bandNameList = spectrumAdapter.getBandNames();
092        bandNameMap = spectrumAdapter.getBandNameMap();
093      }
094      try {
095        setInitialWavenumber(getWavenumberFromChannelIndex(0));
096      } 
097      catch (Exception e) {
098        e.printStackTrace();
099        System.out.println("could not initialize initial wavenumber");
100      }
101    }
102
103    if (swathAdapter != null) {
104      this.swathSelect = swathAdapter.getDefaultSubset();
105      if (spectrumAdapter != null) {
106        spectrumAdapter.setRangeProcessor(swathAdapter.getRangeProcessor());
107      }
108    }
109
110    this.sensorName = sensorName;
111    this.platformName = platformName;
112  }
113
114  public MultiSpectralData(SwathAdapter swathAdapter, SpectrumAdapter spectrumAdapter,
115                           String sensorName, String platformName) {
116    this(swathAdapter, spectrumAdapter, "Radiance", "BrightnessTemp", sensorName, platformName);
117  }
118
119  public MultiSpectralData(SwathAdapter swathAdapter, SpectrumAdapter spectrumAdapter) {
120    this(swathAdapter, spectrumAdapter, null, null);
121  }
122
123  public MultiSpectralData() {
124    this(null, null, null, null);
125  }
126
127  public FlatField getSpectrum(int[] coords) 
128      throws Exception, VisADException, RemoteException {
129    if (coords == null) return null;
130    if (spectrumAdapter == null) return null;
131    spectrumSelect.put(SpectrumAdapter.x_dim_name, new double[] {(double)coords[0], (double)coords[0], 1.0});
132    spectrumSelect.put(SpectrumAdapter.y_dim_name, new double[] {(double)coords[1], (double)coords[1], 1.0});
133
134    FlatField spectrum = spectrumAdapter.getData(spectrumSelect);
135    return convertSpectrum(spectrum, paramName);
136  }
137
138  public FlatField getSpectrum(RealTuple location) 
139      throws Exception, VisADException, RemoteException {
140    if (spectrumAdapter == null) return null;
141    int[] coords = getSwathCoordinates(location, cs);
142    if (coords == null) return null;
143    spectrumSelect.put(SpectrumAdapter.x_dim_name, new double[] {(double)coords[0], (double)coords[0], 1.0});
144    spectrumSelect.put(SpectrumAdapter.y_dim_name, new double[] {(double)coords[1], (double)coords[1], 1.0});
145
146    FlatField spectrum = spectrumAdapter.getData(spectrumSelect);
147    return convertSpectrum(spectrum, paramName);
148  }
149
150  public FlatField getImage(HashMap subset) 
151    throws Exception, VisADException, RemoteException {
152    FlatField image = swathAdapter.getData(subset);
153    cs = ((RealTupleType) ((FunctionType)image.getType()).getDomain()).getCoordinateSystem();
154
155    int channelIndex = (int) ((double[])subset.get(SpectrumAdapter.channelIndex_name))[0];
156    float channel = spectrumAdapter.getWavenumberFromChannelIndex(channelIndex);
157
158    return convertImage(image, channel, paramName);
159  }
160
161  public FlatField getImage(float channel, HashMap subset) 
162      throws Exception, VisADException, RemoteException {
163    if (spectrumAdapter == null) return getImage(subset);
164    int channelIndex = spectrumAdapter.getChannelIndexFromWavenumber(channel);
165    subset.put(SpectrumAdapter.channelIndex_name, new double[] {(double)channelIndex, (double)channelIndex, 1.0});
166    FlatField image = swathAdapter.getData(subset);
167    cs = ((RealTupleType) ((FunctionType)image.getType()).getDomain()).getCoordinateSystem();
168
169    return convertImage(image, channel, paramName);
170  }
171
172  public FlatField getData(Object subset) throws Exception {
173    return getImage((HashMap)subset);
174  }
175
176  public Set makeDomain(Object subset) throws Exception {
177    throw new Exception("makeDomain unimplented");
178  } 
179
180
181  FlatField convertImage(FlatField image, float channel, String param) 
182            throws Exception {
183    FlatField new_image = null;
184    FunctionType f_type = (FunctionType)image.getType();
185    if (param.equals("BrightnessTemp")) { //- convert radiance to BrightnessTemp
186      FunctionType new_type = new FunctionType(f_type.getDomain(), RealType.getRealType("BrightnessTemp"));
187      new_image = new FlatField(new_type, image.getDomainSet());
188      float[][] values = image.getFloats(false);
189      float[] bt_values = values[0];
190      if (inputParamName == "Radiance") {
191        bt_values = radianceToBrightnessTemp(values[0], channel, platformName, sensorName);
192      }
193      new_image.setSamples(new float[][] {bt_values}, false);
194    }
195    else if (param.equals("Reflectance")) {
196      FunctionType new_type = new FunctionType(f_type.getDomain(), RealType.getRealType("Reflectance"));
197      new_image = new FlatField(new_type, image.getDomainSet());
198      new_image.setSamples(image.getFloats(false), false);
199    }
200    else {
201      new_image = image;
202    }
203    return new_image;
204  }
205
206
207  FlatField convertSpectrum(FlatField spectrum, String param) throws Exception {
208    FlatField new_spectrum = null;
209    FunctionType f_type = (FunctionType) spectrum.getType();
210
211    if (param.equals("BrightnessTemp")) {
212      FunctionType new_type = new FunctionType(f_type.getDomain(), RealType.getRealType("BrightnessTemp"));
213      float[][] channels = ((SampledSet)spectrum.getDomainSet()).getSamples(false);
214      float[][] values = spectrum.getFloats(false);
215      float[] bt_values = values[0];
216      if (inputParamName == "Radiance") {
217        bt_values = radianceToBrightnessTempSpectrum(values[0], channels[0], platformName, sensorName);
218      }
219      new_spectrum = new FlatField(new_type, spectrum.getDomainSet());
220      new_spectrum.setSamples(new float[][] {bt_values}, true);
221    }
222    else if (param.equals("Reflectance")) {
223      FunctionType new_type = new FunctionType(f_type.getDomain(), RealType.getRealType("Reflectance"));
224      new_spectrum = new FlatField(new_type, spectrum.getDomainSet());
225      new_spectrum.setSamples(spectrum.getFloats(false), false);
226    }
227    else {
228      new_spectrum = spectrum;
229    }
230    return new_spectrum;
231  }
232
233  protected void setDataRange(float[] range) {
234    dataRange = range;
235  }
236
237  public float[] getDataRange() {
238    return dataRange;
239  }
240
241  public String getParameter() {
242    return paramName;
243  }
244
245  public String getName() {
246    return name;
247  }
248
249  public CoordinateSystem getCoordinateSystem() {
250    return cs;
251  }
252
253  public void setCoordinateSystem(CoordinateSystem cs) {
254    this.cs = cs;
255  }
256
257  public boolean hasBandNames() {
258    return hasBandNames;
259  }
260
261  public ArrayList<String> getBandNames() {
262    return bandNameList;
263  }
264
265  public HashMap<String, Float> getBandNameMap() {
266    return bandNameMap;
267  }
268
269  public String getBandNameFromWaveNumber(float channel) {
270    String bandName = null;
271    Iterator iter = bandNameMap.keySet().iterator();
272    while (iter.hasNext()) {
273       String key = (String) iter.next();
274       float mapVal = ((Float)bandNameMap.get(key)).floatValue();
275       if (channel == mapVal) {
276         bandName = key;
277         break;
278       }
279    }
280    return bandName;
281  }
282
283  public void setInitialWavenumber(float val) {
284    init_wavenumber = val;
285    if (hasBandNames) {
286      init_bandName = getBandNameFromWaveNumber(init_wavenumber);
287    }
288  }
289
290  public int[] getSwathCoordinates(RealTuple location, CoordinateSystem cs) 
291      throws VisADException, RemoteException {
292    if (location == null) return null;
293    if (cs == null) return null;
294    Real[] comps = location.getRealComponents();
295    //- trusted: latitude:0, longitude:1
296    float lon = (float) comps[1].getValue();
297    float lat = (float) comps[0].getValue();
298    if (lon < -180) lon += 360f;
299    if (lon > 180) lon -= 360f;
300    float[][] xy = cs.fromReference(new float[][] {{lon}, {lat}});
301    if ((Float.isNaN(xy[0][0])) || Float.isNaN(xy[1][0])) return null;
302    Set domain = swathAdapter.getSwathDomain();
303    int[] idx = domain.valueToIndex(xy);
304    xy = domain.indexToValue(idx);
305    int[] coords = new int[2];
306    coords[0] = (int) xy[0][0];
307    coords[1] = (int) xy[1][0];
308    if ((coords[0] < 0)||(coords[1] < 0)) return null;
309    return coords;
310  }
311
312  public RealTuple getEarthCoordinates(float[] xy)
313      throws VisADException, RemoteException {
314    float[][] tup = cs.toReference(new float[][] {{xy[0]}, {xy[1]}});
315    return new RealTuple(RealTupleType.SpatialEarth2DTuple, new double[] {(double)tup[0][0], (double)tup[1][0]});
316  }
317
318  public int getChannelIndexFromWavenumber(float channel) throws Exception {
319    return spectrumAdapter.getChannelIndexFromWavenumber(channel);
320  }
321
322  public float getWavenumberFromChannelIndex(int index) throws Exception {
323    return spectrumAdapter.getWavenumberFromChannelIndex(index);
324  }
325
326  public Rectangle2D getLonLatBoundingBox(CoordinateSystem cs) {
327    return null;
328  }
329
330  public Rectangle2D getLonLatBoundingBox(HashMap subset) 
331      throws Exception {
332    Set domainSet = swathAdapter.makeDomain(subset);
333    return getLonLatBoundingBox(domainSet);
334  }
335
336  public static Rectangle2D getLonLatBoundingBox(FlatField field) {
337    Set domainSet = field.getDomainSet();
338    return getLonLatBoundingBox(domainSet);
339  }
340
341  public static Rectangle2D getLonLatBoundingBox(Set domainSet) {
342    CoordinateSystem cs = 
343      ((SetType)domainSet.getType()).getDomain().getCoordinateSystem();
344
345    float start0, stop0, start1, stop1;
346    float minLon = Float.MAX_VALUE;
347    float minLat = Float.MAX_VALUE;
348    float maxLon = -Float.MAX_VALUE;
349    float maxLat = -Float.MAX_VALUE;
350
351
352    if (domainSet instanceof Linear2DSet) {
353      Linear1DSet lset = ((Linear2DSet)domainSet).getLinear1DComponent(0);
354      start0 = (float) lset.getFirst();
355      stop0 = (float) lset.getLast();
356      lset = ((Linear2DSet)domainSet).getLinear1DComponent(1);
357      start1 = (float) lset.getFirst();
358      stop1 = (float) lset.getLast();
359
360      float x, y, del_x, del_y;
361      del_x = (stop0 - start0)/4;
362      del_y = (stop1 - start1)/4;
363      x = start0;
364      y = start1;
365      try {
366        for (int j=0; j<5; j++) {
367          y = start1+j*del_y;
368          for (int i=0; i<5; i++) {
369            x = start0+i*del_x;
370            float[][] lonlat = cs.toReference(new float[][] {{x}, {y}});
371            float lon = lonlat[0][0];
372            float lat = lonlat[1][0];
373            if (lon < minLon) minLon = lon;
374            if (lat < minLat) minLat = lat;
375            if (lon > maxLon) maxLon = lon;
376            if (lat > maxLat) maxLat = lat;
377          }
378        }
379       } catch (Exception e) {
380       }
381    }
382    else if (domainSet instanceof Gridded2DSet) {
383      int[] lens = ((Gridded2DSet)domainSet).getLengths();
384      start0 = 0f;
385      start1 = 0f;
386      stop0 = (float) lens[0];
387      stop1 = (float) lens[1];
388
389      float x, y, del_x, del_y;
390      del_x = (stop0 - start0)/4;
391      del_y = (stop1 - start1)/4;
392      x = start0;
393      y = start1;
394      try {
395        for (int j=0; j<5; j++) {
396          y = start1+j*del_y;
397          for (int i=0; i<5; i++) {
398            x = start0+i*del_x;
399            float[][] lonlat = ((Gridded2DSet)domainSet).gridToValue(new float[][] {{x}, {y}});
400            float lon = lonlat[0][0];
401            float lat = lonlat[1][0];
402            if ((lon > 180 || lon < -180) || (lat > 90 || lat < -90)) continue;
403            if (lon < minLon) minLon = lon;
404            if (lat < minLat) minLat = lat;
405            if (lon > maxLon) maxLon = lon;
406            if (lat > maxLat) maxLat = lat;
407          }
408        }
409      } catch (Exception e) {
410      }
411    }
412    
413
414    float del_lon = maxLon - minLon;
415    float del_lat = maxLat - minLat;
416
417    return new Rectangle2D.Float(minLon, minLat, del_lon, del_lat);
418  }
419
420  public float[] radianceToBrightnessTemp(float[] values, float channelValue) {
421    float c1=1.191066E-5f;           //- mW/m2/ster/cm^-4
422    float c2=1.438833f;              //- K*cm
423    float nu = channelValue;         //- nu: wavenumber
424    float B, K, BT;
425
426    int n_values = values.length;
427    float[] new_values = new float[n_values];
428    for (int i=0; i<n_values;i++) {
429      B = values[i];
430      K = (c1*nu*nu*nu)/B;
431      if (K == 0.0) {
432        BT = B;
433      } 
434      else {
435        BT = c2*nu/((float) (Math.log((double)((c1*nu*nu*nu)/B)+1.0f)) );
436      }
437      if (BT < 0.01) BT = Float.NaN;
438      new_values[i] = BT;
439    }
440    return new_values;
441  }
442
443  public float[] radianceToBrightnessTemp(float[] values, float channelValue, String platformName, String sensorName) 
444     throws Exception {
445    float[] new_values = null;
446
447    if (sensorName == null) {
448      new_values = radianceToBrightnessTemp(values, channelValue);
449    }
450    else if (sensorName == "MODIS") {
451      int channelIndex = spectrumAdapter.getChannelIndexFromWavenumber(channelValue);
452      int band_number = MODIS_L1B_Utility.emissive_indexToBandNumber(channelIndex);
453      new_values = MODIS_L1B_Utility.modis_radiance_to_brightnessTemp(platformName, band_number, values);
454    }
455    return new_values;
456  }
457
458  public float[] radianceToBrightnessTempSpectrum(float[] values, float[] channelValues) {
459    //- Converts radiances [mW/ster/m2/cm^-1] to BT [K]
460    //-  Input: nu  array of wavenmbers [cm^-1]
461    //-          B   radiances [mW/ster/m2/cm^-1]
462    //-  Output: bt brightness temperature in [K]
463    //-   Paolo Antonelli
464    //-   Wed Feb 25 16:43:05 CST 1998
465
466    float c1=1.191066E-5f;           //- mW/m2/ster/cm^-4
467    float c2=1.438833f;              //- K*cm
468
469    float nu;                        //- wavenumber
470    float B, BT;
471
472    int n_values = values.length;
473    float[] new_values = new float[n_values];
474    for (int i=0; i<n_values; i++) {
475      nu = channelValues[i];
476      B = values[i];
477      BT = c2*nu/((float) (Math.log(((c1*nu*nu*nu)/B)+1.0f)) );
478      new_values[i] = BT;
479    }
480    return new_values;
481  }
482
483
484  public float[] radianceToBrightnessTempSpectrum(float[] values, float[] channelValues,
485                                 String platformName, String sensorName) 
486     throws Exception
487  {
488    float[] new_values = null;
489
490    if (sensorName == null) {
491      new_values =  radianceToBrightnessTempSpectrum(values, channelValues);
492    }
493    else if (sensorName == "MODIS") {
494      new_values = new float[values.length];
495      for (int k=0; k<new_values.length; k++) {
496        int channelIndex = spectrumAdapter.getChannelIndexFromWavenumber(channelValues[k]);
497        int band_number = MODIS_L1B_Utility.emissive_indexToBandNumber(channelIndex);
498        float[] tmp = new float[1];
499        tmp[0] = values[k];
500        new_values[k] = (MODIS_L1B_Utility.modis_radiance_to_brightnessTemp(platformName, band_number, tmp))[0];
501      }
502    }
503
504    return new_values;
505  }
506
507  public HashMap getDefaultSubset() {
508    HashMap subset = swathAdapter.getDefaultSubset();
509    double chanIdx=0;
510
511    try {
512       chanIdx = spectrumAdapter.getChannelIndexFromWavenumber(init_wavenumber);
513    }
514    catch (Exception e) {
515      System.out.println("couldn't get chanIdx, using zero");
516    }
517      
518    subset.put(SpectrumAdapter.channelIndex_name, new double[] {chanIdx, chanIdx, 1});
519    return subset;
520  }
521 
522
523  public SpectrumAdapter getSpectrumAdapter() {
524    return spectrumAdapter;
525  }
526}