001 /*
002 * $Id: MultiSpectralData.java,v 1.22 2012/02/19 17:35:41 davep Exp $
003 *
004 * This file is part of McIDAS-V
005 *
006 * Copyright 2007-2012
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
031 package edu.wisc.ssec.mcidasv.data.hydra;
032
033 import java.awt.geom.Rectangle2D;
034 import java.rmi.RemoteException;
035 import java.util.ArrayList;
036 import java.util.HashMap;
037 import java.util.Iterator;
038
039 import visad.CoordinateSystem;
040 import visad.FlatField;
041 import visad.FunctionType;
042 import visad.Gridded2DSet;
043 import visad.Linear1DSet;
044 import visad.Linear2DSet;
045 import visad.Real;
046 import visad.RealTuple;
047 import visad.RealTupleType;
048 import visad.RealType;
049 import visad.SampledSet;
050 import visad.Set;
051 import visad.SetType;
052 import visad.VisADException;
053
054 public 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 float lonA = Float.NaN;
362 float lonB = Float.NaN;
363 float lonC = Float.NaN;
364 float lonD = Float.NaN;
365 float latA = Float.NaN;
366 float latB = Float.NaN;
367 float latC = Float.NaN;
368 float latD = Float.NaN;
369 del_x = (stop0 - start0)/20;
370 del_y = (stop1 - start1)/20;
371
372 x = start0;
373 y = start1;
374 try {
375 for (int j=0; j<21; j++) {
376 y = start1+j*del_y;
377 for (int i=0; i<21; i++) {
378 x = start0 + i*del_x;
379 float[][] lonlat = cs.toReference(new float[][] {{x}, {y}});
380 float lon = lonlat[0][0];
381 float lat = lonlat[1][0];
382 if (!Float.isNaN(lon) && !Float.isNaN(lat)) {
383 lonA = lon;
384 latA = lat;
385 break;
386 }
387 }
388 for (int i=0; i<21; i++) {
389 x = stop0 - i*del_x;
390 float[][] lonlat = cs.toReference(new float[][] {{x}, {y}});
391 float lon = lonlat[0][0];
392 float lat = lonlat[1][0];
393 if (!Float.isNaN(lon) && !Float.isNaN(lat)) {
394 lonB = lon;
395 latB = lat;
396 break;
397 }
398 }
399 if (!Float.isNaN(lonA) && !Float.isNaN(lonB)) {
400 break;
401 }
402 }
403
404 for (int j=0; j<21; j++) {
405 y = stop1-j*del_y;
406 for (int i=0; i<21; i++) {
407 x = start0 + i*del_x;
408 float[][] lonlat = cs.toReference(new float[][] {{x}, {y}});
409 float lon = lonlat[0][0];
410 float lat = lonlat[1][0];
411 if (!Float.isNaN(lon) && !Float.isNaN(lat)) {
412 lonC = lon;
413 latC = lat;
414 break;
415 }
416 }
417 for (int i=0; i<21; i++) {
418 x = stop0 - i*del_x;
419 float[][] lonlat = cs.toReference(new float[][] {{x}, {y}});
420 float lon = lonlat[0][0];
421 float lat = lonlat[1][0];
422 if (!Float.isNaN(lon) && !Float.isNaN(lat)) {
423 lonD = lon;
424 latD = lat;
425 break;
426 }
427 }
428 if (!Float.isNaN(lonC) && !Float.isNaN(lonD)) {
429 break;
430 }
431 }
432
433 float[][] corners = {{lonA,lonB,lonC,lonD},{latA,latB,latC,latD}};
434 for (int k=0; k<corners[0].length; k++) {
435 float lon = corners[0][k];
436 float lat = corners[1][k];
437 if (lon < minLon) minLon = lon;
438 if (lat < minLat) minLat = lat;
439 if (lon > maxLon) maxLon = lon;
440 if (lat > maxLat) maxLat = lat;
441 }
442 } catch (Exception e) {
443 }
444 }
445 else if (domainSet instanceof Gridded2DSet) {
446 int[] lens = ((Gridded2DSet)domainSet).getLengths();
447 start0 = 0f;
448 start1 = 0f;
449 stop0 = (float) lens[0];
450 stop1 = (float) lens[1];
451
452 float x, y, del_x, del_y;
453 del_x = (stop0 - start0)/10;
454 del_y = (stop1 - start1)/10;
455 x = start0;
456 y = start1;
457 try {
458 for (int j=0; j<11; j++) {
459 y = start1+j*del_y;
460 for (int i=0; i<11; i++) {
461 x = start0+i*del_x;
462 float[][] lonlat = ((Gridded2DSet)domainSet).gridToValue(new float[][] {{x}, {y}});
463 float lon = lonlat[0][0];
464 float lat = lonlat[1][0];
465 if ((lon > 180 || lon < -180) || (lat > 90 || lat < -90)) continue;
466 if (lon < minLon) minLon = lon;
467 if (lat < minLat) minLat = lat;
468 if (lon > maxLon) maxLon = lon;
469 if (lat > maxLat) maxLat = lat;
470 }
471 }
472 } catch (Exception e) {
473 }
474 }
475
476
477 float del_lon = maxLon - minLon;
478 float del_lat = maxLat - minLat;
479
480 return new Rectangle2D.Float(minLon, minLat, del_lon, del_lat);
481 }
482
483 public float[] radianceToBrightnessTemp(float[] values, float channelValue) {
484 float c1=1.191066E-5f; //- mW/m2/ster/cm^-4
485 float c2=1.438833f; //- K*cm
486 float nu = channelValue; //- nu: wavenumber
487 float B, K, BT;
488
489 int n_values = values.length;
490 float[] new_values = new float[n_values];
491 for (int i=0; i<n_values;i++) {
492 B = values[i];
493 K = (c1*nu*nu*nu)/B;
494 if (K == 0.0) {
495 BT = B;
496 }
497 else {
498 BT = c2*nu/((float) (Math.log((double)((c1*nu*nu*nu)/B)+1.0f)) );
499 }
500 if (BT < 0.01) BT = Float.NaN;
501 new_values[i] = BT;
502 }
503 return new_values;
504 }
505
506 public float[] radianceToBrightnessTemp(float[] values, float channelValue, String platformName, String sensorName)
507 throws Exception {
508 float[] new_values = null;
509
510 if (sensorName == null) {
511 new_values = radianceToBrightnessTemp(values, channelValue);
512 }
513 else if (sensorName == "MODIS") {
514 int channelIndex = spectrumAdapter.getChannelIndexFromWavenumber(channelValue);
515 int band_number = MODIS_L1B_Utility.emissive_indexToBandNumber(channelIndex);
516 new_values = MODIS_L1B_Utility.modis_radiance_to_brightnessTemp(platformName, band_number, values);
517 }
518 return new_values;
519 }
520
521 public float[] radianceToBrightnessTempSpectrum(float[] values, float[] channelValues) {
522 //- Converts radiances [mW/ster/m2/cm^-1] to BT [K]
523 //- Input: nu array of wavenmbers [cm^-1]
524 //- B radiances [mW/ster/m2/cm^-1]
525 //- Output: bt brightness temperature in [K]
526 //- Paolo Antonelli
527 //- Wed Feb 25 16:43:05 CST 1998
528
529 float c1=1.191066E-5f; //- mW/m2/ster/cm^-4
530 float c2=1.438833f; //- K*cm
531
532 float nu; //- wavenumber
533 float B, BT;
534
535 int n_values = values.length;
536 float[] new_values = new float[n_values];
537 for (int i=0; i<n_values; i++) {
538 nu = channelValues[i];
539 B = values[i];
540 BT = c2*nu/((float) (Math.log(((c1*nu*nu*nu)/B)+1.0f)) );
541 new_values[i] = BT;
542 }
543 return new_values;
544 }
545
546
547 public float[] radianceToBrightnessTempSpectrum(float[] values, float[] channelValues,
548 String platformName, String sensorName)
549 throws Exception
550 {
551 float[] new_values = null;
552
553 if (sensorName == null) {
554 new_values = radianceToBrightnessTempSpectrum(values, channelValues);
555 }
556 else if (sensorName == "MODIS") {
557 new_values = new float[values.length];
558 for (int k=0; k<new_values.length; k++) {
559 int channelIndex = spectrumAdapter.getChannelIndexFromWavenumber(channelValues[k]);
560 int band_number = MODIS_L1B_Utility.emissive_indexToBandNumber(channelIndex);
561 float[] tmp = new float[1];
562 tmp[0] = values[k];
563 new_values[k] = (MODIS_L1B_Utility.modis_radiance_to_brightnessTemp(platformName, band_number, tmp))[0];
564 }
565 }
566
567 return new_values;
568 }
569
570 public HashMap getDefaultSubset() {
571 HashMap subset = swathAdapter.getDefaultSubset();
572 double chanIdx=0;
573
574 try {
575 chanIdx = spectrumAdapter.getChannelIndexFromWavenumber(init_wavenumber);
576 }
577 catch (Exception e) {
578 System.out.println("couldn't get chanIdx, using zero");
579 }
580
581 subset.put(SpectrumAdapter.channelIndex_name, new double[] {chanIdx, chanIdx, 1});
582 return subset;
583 }
584
585
586 public SpectrumAdapter getSpectrumAdapter() {
587 return spectrumAdapter;
588 }
589 }