Source code

001/*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2017
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.adt;
030
031import java.io.IOException;
032import java.lang.Math;
033import java.lang.String;
034import java.rmi.RemoteException;
035import java.util.List;
036
037import org.slf4j.Logger;
038import org.slf4j.LoggerFactory;
039
040import ucar.unidata.data.grid.GridUtil;
041import ucar.unidata.util.IOUtil;
042import ucar.unidata.util.StringUtil;
043import visad.FlatField;
044import visad.VisADException;
045
046public class ReadIRImage {
047
048    private static final Logger logger = LoggerFactory.getLogger(ReadIRImage.class);
049    
050    public ReadIRImage() {
051    }
052    
053    public static class ImgCoeffs {
054        String sat_id;
055        int sat_num;
056        int chan;
057        int det;
058        float scal_m;
059        float scal_b;
060        float side;
061        float conv_n;
062        float conv_a;
063        float conv_b;
064        float conv_g;
065
066        public ImgCoeffs(List<String> toks) {
067            sat_id = toks.get(0);
068            sat_num = Integer.parseInt(toks.get(1));
069            chan = Integer.parseInt(toks.get(2));
070            det = Integer.parseInt(toks.get(3));
071            scal_m = Float.parseFloat(toks.get(4));
072            scal_b = Float.parseFloat(toks.get(5));
073            side = Float.parseFloat(toks.get(6));
074            conv_n = Float.parseFloat(toks.get(7));
075            conv_a = Float.parseFloat(toks.get(8));
076            conv_b = Float.parseFloat(toks.get(9));
077            conv_g = Float.parseFloat(toks.get(10));
078        }
079        
080        public String getSat_id() {
081            return sat_id;
082        }
083    }
084    
085    public static void ReadIRDataFile(FlatField satgrid,
086                                      float cenlat,
087                                      float cenlon,
088                                      int SatelliteID,
089                                      int satChannel,
090                                      boolean isTemperature)
091        throws IOException, VisADException, RemoteException
092    {
093        // Retrieve temperatures from image. This to be done in IDV
094        GridUtil.Grid2D g2d = GridUtil.makeGrid2D(satgrid);
095        float[][] temps = null;
096        float[][][] satimage = null;
097        int numx = 200;
098        int numy = 200;
099        float[][] LocalLatitude = new float[200][200];
100        float[][] LocalLongitude = new float[200][200];
101        float[][] LocalTemperature = new float[200][200];
102        
103        // now spatial subset numx by numy
104        GridUtil.Grid2D g2d1 = spatialSubset(g2d, cenlat, cenlon, numx, numy);
105        
106        satimage = g2d1.getvalues();
107        float[][] temp0 = satimage[0];
108        
109        if (isTemperature) {
110            temps = temp0;
111        } else {
112            temps = im_gvtota(numx, numy, temp0, SatelliteID, satChannel);
113        }
114        
115        Data.IRData_NumberRows = numy;
116        Data.IRData_NumberColumns = numx;
117        Data.IRData_CenterLatitude = cenlat;
118        Data.IRData_CenterLongitude = cenlon;
119        
120        LocalTemperature = temps;
121        LocalLatitude = g2d1.getlats();
122        LocalLongitude = g2d1.getlons();
123        
124        for(int XInc = 0; XInc < numx; XInc++) {
125            for(int YInc = 0; YInc < numy; YInc++) {
126                // must flip x/y to y/x for ADT automated routines
127                Data.IRData_Latitude[YInc][XInc] = LocalLatitude[XInc][YInc];
128                Data.IRData_Longitude[YInc][XInc] = LocalLongitude[XInc][YInc];
129                Data.IRData_Temperature[YInc][XInc] = LocalTemperature[XInc][YInc];
130            }
131        }
132        int CenterXPos = Data.IRData_NumberColumns / 2;
133        int CenterYPos = Data.IRData_NumberRows / 2;
134        
135        double LocalValue[] =
136            Functions.distance_angle(
137                Data.IRData_Latitude[CenterYPos][CenterXPos],
138                Data.IRData_Longitude[CenterYPos][CenterXPos],
139                Data.IRData_Latitude[CenterYPos+1][CenterXPos],
140                Data.IRData_Longitude[CenterYPos][CenterXPos],
141                1);
142                
143        Data.IRData_ImageResolution = LocalValue[0];
144        
145        History.IRCurrentRecord.date = Data.IRData_JulianDate;
146        History.IRCurrentRecord.time = Data.IRData_HHMMSSTime;
147        History.IRCurrentRecord.latitude = Data.IRData_CenterLatitude;
148        History.IRCurrentRecord.longitude = Data.IRData_CenterLongitude;
149        History.IRCurrentRecord.sattype = SatelliteID;
150        
151        int RetVal[] = Functions.adt_oceanbasin(Data.IRData_CenterLatitude,
152                                                Data.IRData_CenterLongitude);
153        Env.DomainID = RetVal[1];
154        // System.out.printf("lat=%f lon=%f domainID=%d\n",Data.IRData_CenterLatitude,Data.IRData_CenterLongitude,Env.DomainID);
155    }
156    
157    private static GridUtil.Grid2D spatialSubset(GridUtil.Grid2D g2d,
158                                                 float cenlat,
159                                                 float cenlon,
160                                                 int numx,
161                                                 int numy)
162    {
163        float[][] lats = g2d.getlats();
164        float[][] lons = g2d.getlons();
165        float[][][] values = g2d.getvalues();
166        float[][] slats = new float[numx][numy];
167        float[][] slons = new float[numx][numy];
168        float[][][] svalues = new float[1][numx][numy];
169        
170        int ly = lats[0].length;
171        int ly0 = ly / 2;
172        int lx = lats.length;
173        logger.debug("lenx: {}, leny: {}", lx, ly);
174        int lx0 = lx / 2;
175        int ii = numx / 2, jj = numy / 2;
176        
177        for (int j = 0; j < ly - 1; j++) {
178            if (Float.isNaN(lats[lx0][j])) {
179                continue;
180            }
181            if ((lats[lx0][j] > cenlat) && (lats[lx0][j + 1] < cenlat)) {
182                jj = j;
183            }
184        }
185        for (int i = 0; i < lx - 1; i++) {
186            if (Float.isNaN(lons[i][ly0])) {
187                continue;
188            }
189            if ((lons[i][ly0] < cenlon) && (lons[i + 1][ly0] > cenlon)) {
190                ii = i;
191            }
192        }
193        int startx = ii - (numx / 2 - 1);
194        int starty = jj - (numy / 2 - 1);
195        logger.debug("startx: {}, starty: {}", startx, starty);
196        logger.debug("numx: {}, numy: {}", numx, numy);
197        
198        if (startx < 0) {
199            startx = 0;
200        }
201        if (starty < 0) {
202            starty = 0;
203        }
204        for (int i = 0; i < numx; i++) {
205            for (int j = 0; j < numy; j++) {
206                try {
207                    slats[i][j] = lats[i + startx][j + starty];
208                } catch (ArrayIndexOutOfBoundsException aioobe) {
209                    slats[i][j] = Float.NaN;
210                }
211                try {
212                    slons[i][j] = lons[i + startx][j + starty];
213                } catch (ArrayIndexOutOfBoundsException aioobe) {
214                    slats[i][j] = Float.NaN;
215                }
216                try {
217                    svalues[0][i][j] = values[0][i + startx][j + starty];
218                } catch (ArrayIndexOutOfBoundsException aioobe) {
219                    slats[i][j] = Float.NaN;
220                }
221            }
222        }
223        return new GridUtil.Grid2D(slats, slons, svalues);
224    }
225    
226    /**
227     * im_gvtota
228     *
229     * This subroutine converts GVAR counts to actual temperatures based on the
230     * current image set in IM_SIMG.
231     *
232     * im_gvtota ( int *nvals, unsigned int *gv, float *ta, int *iret )
233     *
234     * Input parameters: *nvals int Number of values to convert *gv int Array of
235     * GVAR count values
236     *
237     * Output parameters: *ta float Array of actual temperatures *iret int
238     * Return value = -1 - could not open table = -2 - could not find match
239     *
240     *
241     * Log: D.W.Plummer/NCEP 02/03 D.W.Plummer/NCEP 06/03 Add coeff G for 2nd
242     * order poly conv T. Piper/SAIC 07/06 Added tmpdbl to eliminate warning
243     */
244    private static float[][] im_gvtota(int nx,
245                                       int ny,
246                                       float[][] gv,
247                                       int imsorc,
248                                       int imtype)
249    {
250        double c1 = 1.191066E-5;
251        double c2 = 1.438833;
252          
253        int ii;
254        int ip;
255        int chan;
256        int found;
257        double Rad;
258        double Teff;
259        double tmpdbl;
260        float[][] ta = new float[nx][ny];
261        String fp = "/ucar/unidata/data/storm/ImgCoeffs.tbl";
262          
263        for (ii = 0; ii < nx; ii++) {
264            for (int jj = 0; jj < ny; jj++) {
265                ta[ii][jj] = Float.NaN;
266            }
267        }
268        
269        // Read in coefficient table if necessary.
270        String s = null;
271        try {
272            s = IOUtil.readContents(fp);
273        } catch (Exception re) {
274            logger.warn("Failed to read coefficient table", re);
275        }
276        
277        int i = 0;
278        ImgCoeffs[] ImageConvInfo = new ImgCoeffs[50];
279        for (String line : StringUtil.split(s, "\n", true, true)) {
280            if (line.startsWith("!")) {
281                continue;
282            }
283            List<String> stoks = StringUtil.split(line, " ", true, true);
284            ImageConvInfo[i] = new ImgCoeffs(stoks);
285            i++;
286        }
287        int nImgRecs = i;
288        found = 0;
289        ii = 0;
290        while ((ii < nImgRecs) && (found == 0)) {
291            int tmp = ImageConvInfo[ii].chan - 1;
292            tmpdbl = (double)(tmp * tmp);
293            chan = G_NINT(tmpdbl);
294            if ((imsorc == ImageConvInfo[ii].sat_num) && (imtype == chan)) {
295                found = 1;
296            } else {
297                ii++;
298            }
299        }
300        
301        if (found == 0) {
302            return null;
303        } else {
304            ip = ii;
305            for (ii = 0; ii < nx; ii++) {
306                for (int jj = 0; jj < ny; jj++) {
307                     // Convert GVAR count (gv) to Scene Radiance
308                    Rad = ((double) gv[ii][jj] - ImageConvInfo[ip].scal_b) /
309                    // -------------------------------------
310                    ImageConvInfo[ip].scal_m;
311                      
312                    Rad = Math.max(Rad, 0.0);
313                    
314                    // Convert Scene Radiance to Effective Temperature
315                    Teff = (c2 * ImageConvInfo[ip].conv_n)
316                                    /
317                                    /*
318                                     * --------------------------------------------------
319                                     * -----
320                                     */
321                                    (Math.log(1.0
322                                                    + (c1 * Math.pow(ImageConvInfo[ip].conv_n,
323                                                                    3.0)) / Rad));
324                    // Convert Effective Temperature to Temperature
325                    ta[ii][jj] = (float) (ImageConvInfo[ip].conv_a
326                                    + ImageConvInfo[ip].conv_b * Teff + ImageConvInfo[ip].conv_g
327                                    * Teff * Teff);
328                }
329            }
330        }
331        return ta;
332    }
333
334    public static int G_NINT(double x) {
335        return (((x) < 0.0F) ? ((((x) - (float) ((int) (x))) <= -.5f) ? (int) ((x) - .5f)
336               : (int) (x))
337               : ((((x) - (float) ((int) (x))) >= .5f) ? (int) ((x) + .5f)
338                                          : (int) (x)));
339      }
340}