Source code

001/*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2025
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.util;
030
031import java.rmi.RemoteException;
032
033import visad.CommonUnit;
034import visad.Data;
035import visad.FlatField;
036import visad.FunctionType;
037import visad.Linear1DSet;
038import visad.Linear2DSet;
039import visad.MathType;
040import visad.RealTupleType;
041import visad.RealType;
042import visad.Set;
043import visad.SetException;
044import visad.SetType;
045import visad.VisADException;
046
047/**
048 * Utilities for efficiently upscaling and downscaling data referenced on the
049 * GOES-16 2km, 1km and hkm Fixed Grid Frame (FGF).
050 *
051 * <p>Adapted from VisAD's {@link FlatField#resample(Set) resample}, but realizes
052 * efficiencies when target and source domain sets are both rectilinear,
053 * especially avoiding the expansion of the domain samples (getSamples) for
054 * this special case.</p>
055 *
056 * <p>Follows resampling rules outlined in the
057 * <a href="http://www.goes-r.gov/products/ATBDs/baseline/Imagery_v2.0_no_color.pdf">
058 * GOES-16 ABI ATBD for the Cloud and Moisture Imagery Product</a> when
059 * incoming grids are at the their full resolution.
060 * Should be thought of as a grid transfer operation, to be used in conjunction with, 
061 * not necessarily a replacement for, {@code FlatField.resample}.</p>
062 */
063public class GEOSgridUtil {
064    
065    public static RealTupleType LambdaTheta;
066    static {
067        try {
068            LambdaTheta = new RealTupleType(
069                RealType.getRealType("lambda", CommonUnit.radian),
070                RealType.getRealType("theta", CommonUnit.radian));
071        } catch (VisADException e) {
072        }
073    }
074    
075    private static int uniqueID = 0;
076    
077    private static double offsetTolerance = 0.01;
078    
079    private static float accumTolerance = 0.01f;
080    
081    /**
082     * DomainSets of {@code red}, {@code green}, {@code blue} must be (lambda, theta): the intermediate view angle coordinates (radians)
083     * for the Fixed Grid Frame. See {@link #makeGEOSRadiansDomainField}.
084     * TargetSet is taken as that of the first component field ({@code red}).
085     *
086     * @param red
087     * @param green
088     * @param blue
089     *
090     * @return
091     *
092     * @throws VisADException
093     * @throws RemoteException
094     *
095     * @see #makeGEOSRadiansDomainField(FlatField, double, double, double, double)
096     */
097    public static FlatField geosRGB(FlatField red, FlatField green, FlatField blue) throws VisADException, RemoteException {
098        return geosRGB(red, green, blue, (Linear2DSet)red.getDomainSet(), false);
099    }
100    
101    /**
102     *
103     * @param red
104     * @param green
105     * @param blue
106     * @param targetSet
107     *
108     * @return
109     *
110     * @throws VisADException
111     * @throws RemoteException
112     */
113    public static FlatField geosRGB(FlatField red, FlatField green, FlatField blue, Linear2DSet targetSet) throws VisADException, RemoteException {
114        return geosRGB(red, green, blue, targetSet, true);
115    }
116    
117    /**
118     * DomainSets of {@code red}, {@code green}, {@code blue} and {@code targetSet}
119     * must be (lambda, theta): the intermediate view angle coordinates (radians)
120     * for the Fixed Grid Frame.
121     * See {@link #makeGEOSRadiansDomainField(FlatField, double, double, double, double)}.
122     *
123     * @param red
124     * @param green
125     * @param blue
126     * @param targetSet the target, or outgoing domain
127     * @param copy to copy range of resulting RGB FlatField (default=true)
128     *
129     * @return
130     *
131     * @throws VisADException
132     * @throws RemoteException
133     *
134     * @see #makeGEOSRadiansDomainField(FlatField, double, double, double, double)
135     */
136    public static FlatField geosRGB(FlatField red,
137                                    FlatField green,
138                                    FlatField blue,
139                                    Linear2DSet targetSet,
140                                    boolean copy)
141        throws VisADException, RemoteException
142    {
143        
144        Linear2DSet setR = (Linear2DSet)red.getDomainSet();
145        Linear2DSet setG = (Linear2DSet)green.getDomainSet();
146        Linear2DSet setB = (Linear2DSet)blue.getDomainSet();
147        
148        float[] redValues;
149        float[] greenValues;
150        float[] blueValues;
151        
152        RealTupleType newRangeType = new RealTupleType(new RealType[]
153            {RealType.getRealType("redimage_"+uniqueID), RealType.getRealType("greenimage_"+uniqueID), RealType.getRealType("blueimage_"+uniqueID)});
154        
155        FlatField rgb = new FlatField(new FunctionType(((SetType)targetSet.getType()).getDomain(), newRangeType), targetSet);
156        
157        if (targetSet.equals(setR) && targetSet.equals(setB) && targetSet.equals(setG)) {
158            redValues = red.getFloats(false)[0];
159            greenValues = green.getFloats(false)[0];
160            blueValues = blue.getFloats(false)[0];
161        } else {
162            redValues = geosResample(red, targetSet).getFloats(false)[0];
163            greenValues = geosResample(green, targetSet).getFloats(false)[0];
164            blueValues = geosResample(blue, targetSet).getFloats(false)[0];
165        }
166        
167        // For RGB composite, if any NaN -> all NaN
168        for (int k=0; k<redValues.length; k++) {
169            if (Float.isNaN(redValues[k]) || Float.isNaN(greenValues[k]) || Float.isNaN(blueValues[k])) {
170                redValues[k] = Float.NaN;
171                greenValues[k] = Float.NaN;
172                blueValues[k] = Float.NaN;
173            }
174        }
175        
176        rgb.setSamples(new float[][] {redValues, greenValues, blueValues}, copy);
177        
178        uniqueID++;
179        return rgb;
180    }
181    
182    /**
183     * Transforms FlatField with DomainCoordinateSystem
184     * {@code (fgf_x, fgf_y) <-> (Lon,Lat)} based on the Geostationary projection
185     * from fixed grid coordinates to intermediate coordinates view angle
186     * coordinates in radians (lambda, theta).
187     *
188     * @param fltFld
189     *         The incoming FlatField
190     * @param scaleX
191     * @param offsetX
192     *         To transform to fgf_x -> lambda (radians)
193     * @param scaleY
194     * @param offsetY
195     *         To transform to fgf_y -> theta (radians)
196     *
197     * @return FlatField with transformed Domain and original (not copied) Range
198     */
199    public static FlatField makeGEOSRadiansDomainField(FlatField fltFld,
200                                                       double scaleX,
201                                                       double offsetX,
202                                                       double scaleY,
203                                                       double offsetY)
204        throws VisADException, RemoteException
205    {
206        Linear2DSet domainSet = (Linear2DSet)fltFld.getDomainSet();
207        MathType rangeType = ((FunctionType)fltFld.getType()).getRange();
208        float[][] rangeVals = fltFld.getFloats(false);
209        Linear1DSet setX = domainSet.getX();
210        Linear1DSet setY = domainSet.getY();
211        
212        int lenX = setX.getLength();
213        int lenY = setY.getLength();
214        
215        double firstX = setX.getFirst()*scaleX + offsetX;
216        double firstY = setY.getFirst()*scaleY + offsetY;
217        
218        double lastX = setX.getLast()*scaleX + offsetX;
219        double lastY = setY.getLast()*scaleY + offsetY;
220        
221        Linear2DSet dSetRadians = new Linear2DSet(firstX, lastX, lenX, firstY, lastY, lenY);
222        fltFld = new FlatField(new FunctionType(LambdaTheta, rangeType), dSetRadians);
223        fltFld.setSamples(rangeVals, false);
224        return fltFld;
225    }
226    
227    /**
228     * Efficiently upscales or downscales between ABI fields with domainSets on the Fixed Grid Frame.
229     * Can be seen as grid transfer process, to be used in conjunction with, not a replacement for, VisAD resample.
230     *
231     * DomainSets must be (lambda, theta): the intermediate view angle coordinates (radians)
232     * for the Fixed Grid Frame. See {@link #makeGEOSRadiansDomainField}.
233     *
234     * @param fld
235     * @param targetSet
236     *
237     * @return Result field with targetSet domain.
238     *
239     * @throws VisADException
240     * @throws RemoteException
241     *
242     * @see #makeGEOSRadiansDomainField(FlatField, double, double, double, double)
243     */
244    public static FlatField geosResample(FlatField fld, Linear2DSet targetSet)
245        throws VisADException, RemoteException
246    {
247        return geosResample(fld, targetSet, Data.WEIGHTED_AVERAGE);
248    }
249    
250    /**
251     * Efficiently upscales or downscales between ABI fields with domainSets on the Fixed Grid Frame.
252     * Can be seen as grid transfer process, to be used in conjunction with, not a replacement for, VisAD resample.
253     *
254     * DomainSets must be (lambda, theta): the intermediate view angle coordinates (radians)
255     * for the Fixed Grid Frame. See {@link #makeGEOSRadiansDomainField}.
256     *
257     * @param fld
258     * @param targetSet
259     * @param mode VisAD resample mode
260     *
261     * @return Result field with targetSet domain.
262     *
263     * @throws VisADException
264     * @throws RemoteException
265     *
266     * @see #makeGEOSRadiansDomainField(FlatField, double, double, double, double)
267     */
268    public static FlatField geosResample(FlatField fld,
269                                         Linear2DSet targetSet,
270                                         int mode)
271        throws VisADException, RemoteException
272    {
273        double targetStepX = targetSet.getX().getStep();
274        double targetStepY = targetSet.getY().getStep();
275        FunctionType fncType = (FunctionType) fld.getType();
276        RealTupleType setType = ((SetType)targetSet.getType()).getDomain();
277        FunctionType targetType = new FunctionType(setType, fncType.getRange());
278        FlatField target = new FlatField(targetType, targetSet);
279        
280        Linear2DSet setR = null;
281        Set set = fld.getDomainSet();
282        if (set instanceof Linear2DSet) {
283            setR = (Linear2DSet)set;
284        } else {
285            throw new VisADException("FlatField must have Linear2DSet domain, use resample");
286        }
287        
288        /* Just return the incoming field in this case */
289        if (setR.equals(targetSet)) {
290            return fld;
291        }
292        
293        double setStepX = setR.getX().getStep();
294        double setStepY = setR.getY().getStep();
295        double stepRatioX = targetStepX/setStepX;
296        double stepRatioY = targetStepY/setStepY;
297        
298        boolean upsample = (stepRatioX > 1) && (stepRatioY > 1);
299        
300        Linear1DSet setX = ((Linear2DSet)targetSet).getX();
301        Linear1DSet setY = ((Linear2DSet)targetSet).getY();
302        
303        int lenX = setX.getLength();
304        int lenY = setY.getLength();
305        int lenXR = setR.getX().getLength();
306        int lenYR = setR.getY().getLength();
307        
308        boolean ok = upsample;
309        
310        ok = ok && ((Math.abs(setX.getFirst() - setR.getX().getFirst()) % (setR.getX().getStep()/2)) < offsetTolerance &&
311                    (Math.abs(setY.getFirst() - setR.getY().getFirst()) % (setR.getY().getStep()/2)) < offsetTolerance);
312        
313        ok = ok && (Math.abs(stepRatioX - 4.0)*lenX < accumTolerance && Math.abs(stepRatioY - 4.0)*lenY < accumTolerance) ||
314                   (Math.abs(stepRatioX - 2.0)*lenX < accumTolerance && Math.abs(stepRatioY - 2.0)*lenY < accumTolerance);
315        
316        float[][] values = fld.getFloats(false);
317        float[][] targetValues = new float[values.length][targetSet.getLength()];
318        for (int t=0; t< targetValues.length; t++) {
319            java.util.Arrays.fill(targetValues[t], Float.NaN);
320        }
321        
322        float[][] xgrid = valueToGrid(setR.getX(), setX.getSamples(false));
323        float[][] ygrid = valueToGrid(setR.getY(), setY.getSamples(false));
324        
325        if (!upsample || !ok) {
326           interp(ygrid[0], xgrid[0], lenYR, lenXR, values, lenY, lenX, targetValues, mode);
327        } else {
328           int[] xidxs = new int[xgrid[0].length];
329           int[] yidxs = new int[ygrid[0].length];
330           
331           for (int k=0; k<xidxs.length; k++) {
332              float x = xgrid[0][k];
333              xidxs[k] = Float.isNaN(x) ? -1 : (int)Math.floor(x);
334           }
335           for (int k=0; k<yidxs.length; k++) {
336              float y = ygrid[0][k];
337              yidxs[k] = Float.isNaN(y) ? -1 : (int)Math.floor(y);
338           }
339           
340           geosUpsample(yidxs, xidxs, lenYR, lenXR, values, lenY, lenX, targetValues, mode);
341        }
342        
343        target.setSamples(targetValues, false);
344        
345        return target;
346    }
347    
348    static void geosUpsample(int[] yidxs,
349                             int[] xidxs,
350                             int lenY,
351                             int lenX,
352                             float[][] values,
353                             int targetLenY,
354                             int targetLenX,
355                             float[][] targetValues,
356                             int mode)
357    {
358        int tupleDimLen = values.length;
359        for (int j=0; j<targetLenY; j++) {
360            int jR = yidxs[j];
361            if (jR >= 0 && jR < lenY) {
362                for (int i=0; i<targetLenX; i++) {
363                    int iR = xidxs[i];
364                    if (iR >= 0 && iR < lenX) {
365                        int k = j*targetLenX + i;
366                        int kR = jR*lenX + iR;
367                        
368                        if (mode == Data.NEAREST_NEIGHBOR) {
369                            for (int t=0; t<tupleDimLen; t++) {
370                                // Upper right corner (view angle space)
371                                // see PUG
372                                targetValues[t][k] = values[t][kR+lenX];
373                            }
374                        } else {
375                            for (int t=0; t<tupleDimLen; t++) {
376                                if (iR < lenX-1 && jR < lenY-1) { // This check should not be needed for NN
377                                   float val = 0;
378                                   val += values[t][kR];
379                                   val += values[t][kR+1];
380                                   val += values[t][kR+lenX];
381                                   val += values[t][kR+lenX+1];
382                                   targetValues[t][k] = val/4;
383                                }
384                            }
385                        }
386                    }
387                }
388            }
389        }
390    }
391    
392    static void interp(float[] yidxs,
393                       float[] xidxs,
394                       int lenY,
395                       int lenX,
396                       float[][] values,
397                       int targetLenY,
398                       int targetLenX,
399                       float[][] targetValues,
400                       int mode)
401    {
402        int tupleDimLen = values.length;
403        int[] idxs = new int[2];
404        float[] flts = new float[4];
405        
406        for (int j=0; j<targetLenY; j++) {
407            float y = yidxs[j];
408            int jR = Float.isNaN(y) ? -1 : (int) Math.floor(y);
409            if (jR >= 0 && jR < lenY-1) {
410                
411                for (int i=0; i<targetLenX; i++) {
412                    float x = xidxs[i];
413                    int iR = Float.isNaN(x) ? -1 : (int) Math.floor(x);
414                    if (iR >= 0 && iR < lenX-1) {
415                        int k = j*targetLenX + i;
416                        int kR = jR*lenX + iR;
417                        
418                        float dx00 = xidxs[i] - iR;
419                        float dy00 = yidxs[j] - jR;
420                        
421                        float dx01 = xidxs[i] - (iR + 1);
422                        float dy01 = yidxs[j] - jR;
423                        
424                        float dx10 = xidxs[i] - iR;
425                        float dy10 = yidxs[j] - (jR+1);
426                        
427                        float dx11 = xidxs[i] - (iR+1);
428                        float dy11 = yidxs[j] - (jR+1);
429                        
430                        float dst00 = (float) (dx00*dx00 + dy00*dy00);
431                        float dst01 = (float) (dx01*dx01 + dy01*dy01);
432                        float dst10 = (float) (dx10*dx10 + dy10*dy10);
433                        float dst11 = (float) (dx11*dx11 + dy11*dy11);
434                        
435                        if (mode == Data.NEAREST_NEIGHBOR) {
436                            flts[0] = dst00;
437                            flts[1] = dst01;
438                            flts[2] = dst10;
439                            flts[3] = dst11;
440                            minmax(flts, 4, idxs);
441                            
442                            for (int t=0; t<tupleDimLen; t++) {
443                                float val = Float.NaN;
444                                if (idxs[0] == 0) {
445                                    val = values[t][kR];
446                                } else if (idxs[0] == 1) {
447                                    val = values[t][kR+1];
448                                } else if (idxs[0] == 2) {
449                                    val = values[t][kR+lenX];
450                                } else if (idxs[0] == 3) {
451                                    val = values[t][kR+lenX+1];
452                                }
453                                targetValues[t][k] = val;
454                            }
455                        } else {
456                            float w00 = 1/dst00;
457                            float w01 = 1/dst01;
458                            float w10 = 1/dst10;
459                            float w11 = 1/dst11;
460                            float sum = w00 + w01 + w10 + w11;
461                            
462                            if (Float.isInfinite(w00)) {
463                                sum = 1f;
464                                w00 = 1f;
465                                w01 = 0;
466                                w10 = 0;
467                                w11 = 0;             
468                            } 
469                            else if (Float.isInfinite(w01)) {
470                                sum = 1f;
471                                w00 = 0;
472                                w01 = 1f;
473                                w10 = 0;
474                                w11 = 0;              
475                            }
476                            else if (Float.isInfinite(w10)) {
477                                sum = 1f;
478                                w00 = 0;
479                                w01 = 0;
480                                w10 = 1f;
481                                w11 = 0;                        
482                            }
483                            else if (Float.isInfinite(w11)) {
484                                sum = 1f;
485                                w00 = 0;
486                                w01 = 0;
487                                w10 = 0;
488                                w11 = 1f;  
489                            }
490                            
491                            for (int t=0; t<tupleDimLen;t++) {
492                                float val = 0;
493                                val += w00*values[t][kR];
494                                val += w01*values[t][kR+1];
495                                val += w10*values[t][kR+lenX];
496                                val += w11*values[t][kR+lenX+1];
497                                targetValues[t][k] = val/sum;
498                            }
499                        }
500                    }
501                }
502            }
503        }
504    }
505    
506    public static float[] minmax(float[] values, int length, int[] indexes) {
507        float min =  Float.MAX_VALUE;
508        float max = -Float.MAX_VALUE;
509        int minIdx = 0;
510        int maxIdx = 0;
511        for (int k = 0; k < length; k++) {
512            float val = values[k];
513            if ((val == val) && (val < Float.POSITIVE_INFINITY) && (val > Float.NEGATIVE_INFINITY)) {
514                if (val < min) {
515                    min = val;
516                    minIdx = k;
517                }
518                if (val > max) {
519                    max = val;
520                    maxIdx = k;
521                }
522            }
523        }
524        if (indexes != null) {
525            indexes[0] = minIdx;
526            indexes[1] = maxIdx;
527        }
528        return new float[] {min, max};
529    }
530    
531    /**
532     * Adapted from VisAD except NaN is returned unless:
533     * val >= First and val <= Last, First > Last
534     * val <= First and val >= Last, Last > First
535     * i.e., clamp to inside the interval inclusive
536     *
537     * @param set
538     * @param value
539     *
540     * @return
541     *
542     * @throws VisADException
543     */
544    private static float[][] valueToGrid(Linear1DSet set, float[][] value) throws VisADException {
545        if (value.length != 1) {
546            throw new SetException("Linear1DSet.valueToGrid: value dimension" +
547                                   " should be 1, not " + value.length);
548        }
549        double First = set.getFirst();
550        double Last = set.getLast();
551        double Step = set.getStep();
552        double Invstep = 1.0/Step;
553        int Length = set.getLength();
554        int length = value[0].length;
555        float[][] grid = new float[1][length];
556        float[] grid0 = grid[0];
557        float[] value0 = value[0];
558        float l = (float) First;
559        float h = (float) Last;
560        float v;
561        
562        if (h < l) {
563            float temp = l;
564            l = h;
565            h = temp;
566        }
567        for (int i=0; i<length; i++) {
568            v = value0[i];
569            grid0[i] = (float) ((l <= v && v <= h) ? (v - First) * Invstep : Float.NaN);
570        }
571        return grid;
572    }
573}