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.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 float val = 0; 377 val += values[t][kR]; 378 val += values[t][kR+1]; 379 val += values[t][kR+lenX]; 380 val += values[t][kR+lenX+1]; 381 targetValues[t][k] = val/4; 382 } 383 } 384 } 385 } 386 } 387 } 388 } 389 390 static void interp(float[] yidxs, 391 float[] xidxs, 392 int lenY, 393 int lenX, 394 float[][] values, 395 int targetLenY, 396 int targetLenX, 397 float[][] targetValues, 398 int mode) 399 { 400 int tupleDimLen = values.length; 401 int[] idxs = new int[2]; 402 float[] flts = new float[4]; 403 404 for (int j=0; j<targetLenY; j++) { 405 float y = yidxs[j]; 406 int jR = Float.isNaN(y) ? -1 : (int) Math.floor(y); 407 if (jR >= 0 && jR < lenY) { 408 409 for (int i=0; i<targetLenX; i++) { 410 float x = xidxs[i]; 411 int iR = Float.isNaN(x) ? -1 : (int) Math.floor(x); 412 if (iR >= 0 && iR < lenX) { 413 int k = j*targetLenX + i; 414 int kR = jR*lenX + iR; 415 416 float dx00 = xidxs[i] - iR; 417 float dy00 = yidxs[j] - jR; 418 419 float dx01 = xidxs[i] - (iR + 1); 420 float dy01 = yidxs[j] - jR; 421 422 float dx10 = xidxs[i] - iR; 423 float dy10 = yidxs[j] - (jR+1); 424 425 float dx11 = xidxs[i] - (iR+1); 426 float dy11 = yidxs[j] - (jR+1); 427 428 float dst00 = (float) (dx00*dx00 + dy00*dy00); 429 float dst01 = (float) (dx01*dx01 + dy01*dy01); 430 float dst10 = (float) (dx10*dx10 + dy10*dy10); 431 float dst11 = (float) (dx11*dx11 + dy11*dy11); 432 433 float sum = 1/dst00 + 1/dst01 + 1/dst10 + 1/dst11; 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 461 for (int t=0; t<tupleDimLen;t++) { 462 float val = 0; 463 val += w00*values[t][kR]; 464 val += w01*values[t][kR+1]; 465 val += w10*values[t][kR+lenX]; 466 val += w11*values[t][kR+lenX+1]; 467 targetValues[t][k] = val/sum; 468 } 469 } 470 } 471 } 472 } 473 } 474 } 475 476 public static float[] minmax(float[] values, int length, int[] indexes) { 477 float min = Float.MAX_VALUE; 478 float max = -Float.MAX_VALUE; 479 int minIdx = 0; 480 int maxIdx = 0; 481 for (int k = 0; k < length; k++) { 482 float val = values[k]; 483 if ((val == val) && (val < Float.POSITIVE_INFINITY) && (val > Float.NEGATIVE_INFINITY)) { 484 if (val < min) { 485 min = val; 486 minIdx = k; 487 } 488 if (val > max) { 489 max = val; 490 maxIdx = k; 491 } 492 } 493 } 494 if (indexes != null) { 495 indexes[0] = minIdx; 496 indexes[1] = maxIdx; 497 } 498 return new float[] {min, max}; 499 } 500 501 /** 502 * Adapted from VisAD except NaN is returned unless: 503 * val >= First and val <= Last, First > Last 504 * val <= First and val >= Last, Last > First 505 * i.e., clamp to inside the interval inclusive 506 * 507 * @param set 508 * @param value 509 * 510 * @return 511 * 512 * @throws VisADException 513 */ 514 private static float[][] valueToGrid(Linear1DSet set, float[][] value) throws VisADException { 515 if (value.length != 1) { 516 throw new SetException("Linear1DSet.valueToGrid: value dimension" + 517 " should be 1, not " + value.length); 518 } 519 double First = set.getFirst(); 520 double Last = set.getLast(); 521 double Step = set.getStep(); 522 double Invstep = 1.0/Step; 523 int Length = set.getLength(); 524 int length = value[0].length; 525 float[][] grid = new float[1][length]; 526 float[] grid0 = grid[0]; 527 float[] value0 = value[0]; 528 float l = (float) First; 529 float h = (float) Last; 530 float v; 531 532 if (h < l) { 533 float temp = l; 534 l = h; 535 h = temp; 536 } 537 for (int i=0; i<length; i++) { 538 v = value0[i]; 539 grid0[i] = (float) ((l <= v && v <= h) ? (v - First) * Invstep : Float.NaN); 540 } 541 return grid; 542 } 543}