Java tutorial
/* * Copyright (C) 2014 by Array Systems Computing Inc. http://www.array.ca * * This program is free software; you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by the Free * Software Foundation; either version 3 of the License, or (at your option) * any later version. * This program is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for * more details. * * You should have received a copy of the GNU General Public License along * with this program; if not, see http://www.gnu.org/licenses/ */ package org.esa.s1tbx.calibration.gpf.calibrators; import com.bc.ceres.core.ProgressMonitor; import org.apache.commons.math3.util.FastMath; import org.esa.s1tbx.calibration.gpf.CalibrationOp; import org.esa.s1tbx.calibration.gpf.support.BaseCalibrator; import org.esa.s1tbx.calibration.gpf.support.Calibrator; import org.esa.s1tbx.calibration.gpf.support.TiePointInterpolator; import org.esa.snap.core.datamodel.Band; import org.esa.snap.core.datamodel.MetadataAttribute; import org.esa.snap.core.datamodel.MetadataElement; import org.esa.snap.core.datamodel.Product; import org.esa.snap.core.datamodel.ProductData; import org.esa.snap.core.datamodel.TiePointGrid; import org.esa.snap.core.gpf.Operator; import org.esa.snap.core.gpf.OperatorException; import org.esa.snap.core.gpf.Tile; import org.esa.snap.core.util.SystemUtils; import org.esa.snap.dataio.envisat.EnvisatAuxReader; import org.esa.snap.engine_utilities.datamodel.AbstractMetadata; import org.esa.snap.engine_utilities.datamodel.DownloadableArchive; import org.esa.snap.engine_utilities.datamodel.OrbitStateVector; import org.esa.snap.engine_utilities.datamodel.Unit; import org.esa.snap.engine_utilities.eo.Constants; import org.esa.snap.engine_utilities.eo.GeoUtils; import org.esa.snap.engine_utilities.gpf.OperatorUtils; import org.esa.snap.engine_utilities.gpf.TileIndex; import org.esa.snap.engine_utilities.util.Maths; import org.esa.snap.engine_utilities.util.Settings; import java.awt.*; import java.io.File; import java.io.IOException; import java.net.URL; import java.nio.file.Path; import java.text.SimpleDateFormat; import java.util.Date; /** * Calibration for ASAR data products. */ public class ASARCalibrator extends BaseCalibrator implements Calibrator { private File externalAuxFile = null; private String auxFile = null; private String productType = null; private String oldXCAFileName = null; // the old XCA file private String newXCAFileName = null; // XCA file for radiometric calibration private String newXCAFilePath = null; // absolute path for XCA file private String[] mdsPolar; // polarizations for the two bands in the product private TiePointGrid incidenceAngle = null; private TiePointGrid slantRangeTime = null; private TiePointGrid latitude = null; private boolean srgrFlag = false; private boolean multilookFlag = false; private boolean antElevCorrFlag = false; private boolean rangeSpreadCompFlag = false; private boolean wideSwathProductFlag = false; private boolean retroCalibrationFlag = false; private boolean applyAntennaPatternCorr = false; private boolean applyRangeSpreadingCorr = false; private double firstLineUTC = 0.0; // in days private double lineTimeInterval = 0.0; // in days private double avgSceneHeight = 0.0; // in m private double rangeSpacing = 0.0; // in m private double azimuthSpacing = 0.0; // in m private double rangeSpreadingCompPower = 0.0; // power in range spreading loss compensation calculation private double halfRangeSpreadingCompPower = 0.0; private double latMax = 0.0; private double delLat = 0.0; private double[] earthRadius = null; // Earth radius for all range lines, in m private final double[] newCalibrationConstant = new double[2]; private double[] oldRefElevationAngle = null; // reference elevation angle for given swath in old aux file, in degree private double[] newRefElevationAngle = null; // reference elevation angle for given swath in new aux file, in degree private float[][] oldAntennaPatternSingleSwath = null; // old antenna pattern gains for single swath product, in dB private float[][] oldAntennaPatternWideSwath = null; // old antenna pattern gains for single swath product, in dB private float[][] newAntennaPatternSingleSwath = null; // new antenna pattern gains for single swath product, in dB private float[][] newAntennaPatternWideSwath = null; // new antenna pattern gains for single swath product, in dB private TiePointInterpolator incidenceTPGInterp = null; private TiePointInterpolator slantRangeTPGInterp = null; private int numMPPRecords; // number of MPP ADSR records private String swath; private OrbitStateVector[] orbitStateVectors = null; private AbstractMetadata.SRGRCoefficientList[] srgrConvParams = null; private static final int numOfGains = 201; // number of antenna pattern gain values for a given swath and // polarization in the aux file // private double refSlantRange = 800000.0; // m // private double halfLightSpeedByRefSlantRange = Constants.halfLightSpeed / refSlantRange; private static final double refSlantRange800km = 800000.0; // m private static final int INVALID_SUB_SWATH_INDEX = -1; public ASARCalibrator() { } /** * Set external auxiliary file. */ @Override public void setExternalAuxFile(File file) { externalAuxFile = file; } /** * Set auxiliary file flag. */ @Override public void setAuxFileFlag(String file) { auxFile = file; } /** * @param op the Calibration Operator * @param srcProduct The source product. * @param tgtProduct The target product. * @param mustPerformRetroCalibration If true, retro-calibration must be performed if it is applicable. * @throws OperatorException The exception. */ @Override public void initialize(final Operator op, final Product srcProduct, final Product tgtProduct, final boolean mustPerformRetroCalibration, final boolean mustUpdateMetadata) throws OperatorException { try { calibrationOp = op; sourceProduct = srcProduct; targetProduct = tgtProduct; absRoot = AbstractMetadata.getAbstractedMetadata(sourceProduct); getProductType(); getSampleType(); getSRGRFlag(); getCalibrationFlags(); getMultilookFlag(); getProductSwath(); setCalibrationFlags(); mdsPolar = OperatorUtils.getProductPolarization(absRoot); getRangeAzimuthSpacing(); numMPPRecords = getNumOfRecordsInMainProcParam(sourceProduct); //??? getTiePointGridData(sourceProduct); checkXCAFileExsitence(mustPerformRetroCalibration); if (srgrFlag) { getSrgrCoeff(); } if (retroCalibrationFlag) { getOldAntennaPattern(); } if (applyAntennaPatternCorr) { getFirstLineTime(); getLineTimeInterval(); getAverageSceneHeight(); getOrbitStateVectors(); getNewAntennaPattern(); computeEarthRadius(); } getNewCalibrationFactor(); setRangeSpreadingLossCompPower(); if (mustUpdateMetadata) { updateTargetProductMetadata(); } } catch (Exception e) { throw new OperatorException(e); } } /** * Get Product ID from MPH. * * @throws OperatorException The exceptions. */ private void getProductType() throws OperatorException { productType = sourceProduct.getProductType(); // product type could be 1P or 1C if (!productType.contains("ASA_IMP_1") && !productType.contains("ASA_IMM_1") && !productType.contains("ASA_APP_1") && !productType.contains("ASA_APM_1") && !productType.contains("ASA_WSM_1") && !productType.contains("ASA_IMG_1") && !productType.contains("ASA_APG_1") && !productType.contains("ASA_IMS_1") && !productType.contains("ASA_APS_1") && !productType.contains("ASA_GM") && !productType.contains("ASA_WSS_1")) { throw new OperatorException(productType + " is not a valid ASAR product type for calibration."); } } /** * Get antenna elevation correction flag and range spreading compensation flag from Metadata. * * @throws Exception The exceptions. */ private void getCalibrationFlags() throws Exception { if (AbstractMetadata.getAttributeBoolean(absRoot, AbstractMetadata.abs_calibration_flag)) { throw new OperatorException("The product has already been calibrated."); } antElevCorrFlag = AbstractMetadata.getAttributeBoolean(absRoot, AbstractMetadata.ant_elev_corr_flag); rangeSpreadCompFlag = AbstractMetadata.getAttributeBoolean(absRoot, AbstractMetadata.range_spread_comp_flag); } /** * Get SRGR flag from the abstracted metadata. * * @throws Exception The exceptions. */ private void getSRGRFlag() throws Exception { srgrFlag = AbstractMetadata.getAttributeBoolean(absRoot, AbstractMetadata.srgr_flag); } /** * Get multilook flag from the abstracted metadata. * * @throws Exception The exceptions. */ private void getMultilookFlag() throws Exception { multilookFlag = AbstractMetadata.getAttributeBoolean(absRoot, AbstractMetadata.multilook_flag); } /** * Get product swath. */ private void getProductSwath() { swath = absRoot.getAttributeString(AbstractMetadata.SWATH); wideSwathProductFlag = swath.contains("WS"); } /** * Set calibration flags. */ private void setCalibrationFlags() { if (antElevCorrFlag) { if (multilookFlag || auxFile != null && auxFile.contains(CalibrationOp.PRODUCT_AUX)) { retroCalibrationFlag = false; //System.out.println("Only constant and incidence angle corrections will be performed for radiometric calibration"); } else { retroCalibrationFlag = true; } } if (auxFile != null && auxFile.contains(CalibrationOp.PRODUCT_AUX)) { applyAntennaPatternCorr = false; } else { applyAntennaPatternCorr = !srgrFlag || retroCalibrationFlag || !antElevCorrFlag; } applyRangeSpreadingCorr = !rangeSpreadCompFlag; } /** * Get SRGR conversion parameters. */ private void getSrgrCoeff() { srgrConvParams = AbstractMetadata.getSRGRCoefficients(absRoot); } /** * Get range spacing from the abstracted metadata. * * @throws Exception The exceptions. */ private void getRangeAzimuthSpacing() throws Exception { rangeSpacing = AbstractMetadata.getAttributeDouble(absRoot, AbstractMetadata.range_spacing); azimuthSpacing = AbstractMetadata.getAttributeDouble(absRoot, AbstractMetadata.azimuth_spacing); } /** * Get XCA file used for the original radiometric calibration. * * @return The complete path to the XCA file. */ private Path getOldXCAFilePath() throws Exception { oldXCAFileName = absRoot.getAttributeString(AbstractMetadata.external_calibration_file); final File localFolder = SystemUtils.getAuxDataPath().resolve("AuxCal").resolve("ENVISAT").toFile(); final URL remotePath = new URL(Settings.getPath("AuxCal.ENVISAT.remotePath")); File xcaFile = new File(localFolder, oldXCAFileName); if (xcaFile.exists()) { // unzipped return xcaFile.toPath(); } else { xcaFile = new File(localFolder, oldXCAFileName + ".zip"); if (xcaFile.exists()) { // zipped return xcaFile.toPath(); } else { final File localFile = new File(localFolder, "ENVISAT_XCA.zip"); final DownloadableArchive archive = new DownloadableArchive(localFile, remotePath); archive.getContentFiles(); } return xcaFile.toPath(); } } /** * Get old antenna pattern gains. */ private void getOldAntennaPattern() throws Exception { final String xcaFilePath = getOldXCAFilePath().toString(); if (wideSwathProductFlag) { oldRefElevationAngle = new double[5]; // reference elevation angles for 5 sub swathes oldAntennaPatternWideSwath = new float[5][numOfGains]; // antenna pattern gain for 5 sub swathes getWideSwathAntennaPatternGainFromAuxData(xcaFilePath, mdsPolar[0], numOfGains, oldRefElevationAngle, oldAntennaPatternWideSwath); } else { oldRefElevationAngle = new double[1]; // reference elevation angle for 1 swath oldAntennaPatternSingleSwath = new float[2][numOfGains]; // antenna pattern gain for 2 bands getSingleSwathAntennaPatternGainFromAuxData(xcaFilePath, swath, mdsPolar, numOfGains, oldRefElevationAngle, oldAntennaPatternSingleSwath); } } /** * Get number of records in Main Processing Params data set. * * @param sourceProduct The source prodict. * @return The number of records. * @throws OperatorException The exceptions. */ static int getNumOfRecordsInMainProcParam(Product sourceProduct) throws OperatorException { final MetadataElement origRoot = AbstractMetadata.getOriginalProductMetadata(sourceProduct); final MetadataElement dsdElem = origRoot.getElement("DSD"); if (dsdElem == null) { throw new OperatorException("DSD not found"); } final MetadataElement dsd = dsdElem.getElement("DSD.3"); if (dsd == null) { throw new OperatorException("DSD not found"); } final MetadataAttribute numRecordsAttr = dsd.getAttribute("num_records"); if (numRecordsAttr == null) { throw new OperatorException("num_records not found"); } int numMPPRecords = numRecordsAttr.getData().getElemInt(); if (numMPPRecords < 1) { throw new OperatorException("Invalid num_records."); } //System.out.println("The number of Main Processing Params records is " + numMPPRecords); return numMPPRecords; } /** * Get incidence angle and slant range time tie point grids. * * @param sourceProduct the source */ private void getTiePointGridData(Product sourceProduct) { slantRangeTime = OperatorUtils.getSlantRangeTime(sourceProduct); incidenceAngle = OperatorUtils.getIncidenceAngle(sourceProduct); latitude = OperatorUtils.getLatitude(sourceProduct); incidenceTPGInterp = new TiePointInterpolator(incidenceAngle); slantRangeTPGInterp = new TiePointInterpolator(slantRangeTime); } /** * Check if old or new XCA file exists. * * @param mustPerformRetroCalibration If true, retro-calibration must be performed if it is applicable. * @throws Exception The exception. */ private void checkXCAFileExsitence(boolean mustPerformRetroCalibration) throws Exception { String oldXCAFilePath = null; if (retroCalibrationFlag) { oldXCAFilePath = getOldXCAFilePath().toString(); if (!isFileExisting(oldXCAFilePath)) { if (mustPerformRetroCalibration) { throw new OperatorException("Cannot find XCA file: " + oldXCAFilePath); } else { retroCalibrationFlag = false; applyAntennaPatternCorr = false; // todo should warn the user that retro-calibration will not be performed } } } if (applyAntennaPatternCorr) { getNewXCAFile(); if (!isFileExisting(newXCAFilePath)) { throw new OperatorException("Cannot find XCA file: " + newXCAFilePath); } if (retroCalibrationFlag && !mustPerformRetroCalibration && newXCAFilePath != null && (oldXCAFilePath.contains(newXCAFilePath) || newXCAFilePath.contains(oldXCAFilePath))) { retroCalibrationFlag = false; applyAntennaPatternCorr = false; } } } /** * Return true if a given file exists, otherwise, false. * * @param filePath The complete path to the given file. * @return Return true if the file exists, false otherwise. */ private boolean isFileExisting(String filePath) { File file = null; final String[] exts = new String[] { "", ".gz", ".zip" }; for (String ext : exts) { file = new File(filePath + ext); if (file.exists()) { return true; } } return false; } private void getNewXCAFile() throws Exception { if (auxFile != null && auxFile.contains(CalibrationOp.EXTERNAL_AUX)) { if (externalAuxFile != null && externalAuxFile.exists()) { if (!externalAuxFile.getName().contains("ASA_XCA")) { throw new OperatorException("Invalid XCA file for ASAR product"); } newXCAFileName = externalAuxFile.getName(); newXCAFilePath = externalAuxFile.getAbsolutePath(); } else { throw new OperatorException("No external auxiliary file is specified."); } } else { final Date startDate = sourceProduct.getStartTime().getAsDate(); final Date endDate = sourceProduct.getEndTime().getAsDate(); final File localFolder = SystemUtils.getAuxDataPath().resolve("AuxCal").resolve("ENVISAT").toFile(); newXCAFileName = findXCAFile(localFolder, startDate, endDate); if (newXCAFileName == null) { final URL remotePath = new URL(Settings.getPath("AuxCal.ENVISAT.remotePath")); final File localFile = new File(localFolder, "ENVISAT_XCA.zip"); final DownloadableArchive archive = new DownloadableArchive(localFile, remotePath); archive.getContentFiles(); newXCAFileName = findXCAFile(localFolder, startDate, endDate); } newXCAFilePath = localFolder.toString() + File.separator + newXCAFileName; } } /** * Find the latest XVA file available. * * @param xcaFileDir The complete path to the XCA file directory. * @param productStartDate The product start date. * @param productEndDate The product end data. * @return The name of the XCA file found. * @throws Exception The exceptions. */ public static String findXCAFile(File xcaFileDir, Date productStartDate, Date productEndDate) throws Exception { final File[] list = xcaFileDir.listFiles(); if (list == null) { return null; } final SimpleDateFormat dateformat = new SimpleDateFormat("yyyyMMdd_HHmmss"); Date latestCreationDate = dateformat.parse("19000101_000000"); String xcaFileName = null; for (File f : list) { final String fileName = f.getName(); if (fileName.length() < 61 || !fileName.substring(0, 10).equals("ASA_XCA_AX")) { continue; } final Date creationDate = dateformat.parse(fileName.substring(14, 29)); final Date validStartDate = dateformat.parse(fileName.substring(30, 45)); final Date validStopDate = dateformat.parse(fileName.substring(46, 61)); if (productStartDate.after(validStartDate) && productEndDate.before(validStopDate) && latestCreationDate.before(creationDate)) { latestCreationDate = creationDate; xcaFileName = fileName; } } return xcaFileName; } /** * Get calibration factor. */ private void getNewCalibrationFactor() { if (newXCAFilePath != null) { getCalibrationFactorFromExternalAuxFile(newXCAFilePath, swath, mdsPolar, productType, newCalibrationConstant); } else { getCalibrationFactorFromMetadata(); } } /** * Get calibration factor from user specified auxiliary data. * * @param auxFilePath The absolute path to the aux file. * @param swath The product swath. * @param mdsPolar The product polarizations. * @param productType The product type. * @param calibrationFactor The calibration factors. */ public static void getCalibrationFactorFromExternalAuxFile(String auxFilePath, String swath, String[] mdsPolar, String productType, double[] calibrationFactor) { final EnvisatAuxReader reader = new EnvisatAuxReader(); try { reader.readProduct(auxFilePath); final int numOfSwaths = 7; String calibrationFactorName; for (int i = 0; i < 2 && mdsPolar[i] != null && mdsPolar[i].length() != 0; i++) { calibrationFactor[i] = 0; if (productType.contains("ASA_IMP_1")) { calibrationFactorName = "ext_cal_im_pri_" + mdsPolar[i]; } else if (productType.contains("ASA_IMM_1")) { calibrationFactorName = "ext_cal_im_med_" + mdsPolar[i]; } else if (productType.contains("ASA_APP_1")) { calibrationFactorName = "ext_cal_ap_pri_" + mdsPolar[i]; } else if (productType.contains("ASA_APM_1")) { calibrationFactorName = "ext_cal_ap_med_" + mdsPolar[i]; } else if (productType.contains("ASA_WS")) { calibrationFactorName = "ext_cal_ws_" + mdsPolar[i]; } else if (productType.contains("ASA_GM1_1")) { calibrationFactorName = "ext_cal_gm_" + mdsPolar[i]; } else if (productType.contains("ASA_IMG_1")) { calibrationFactorName = "ext_cal_im_geo_" + mdsPolar[i]; } else if (productType.contains("ASA_APG_1")) { calibrationFactorName = "ext_cal_ap_geo_" + mdsPolar[i]; } else if (productType.contains("ASA_IMS_1")) { calibrationFactorName = "ext_cal_im_" + mdsPolar[i]; } else if (productType.contains("ASA_APS_1")) { calibrationFactorName = "ext_cal_ap_" + mdsPolar[i]; } else { throw new OperatorException("Invalid ASAR product type."); } final ProductData factorData = reader.getAuxData(calibrationFactorName); final float[] factors = (float[]) factorData.getElems(); if (productType.contains("ASA_WS") || productType.contains("ASA_GM1")) { calibrationFactor[i] = factors[0]; } else { if (factors.length != numOfSwaths) { throw new OperatorException("Incorrect array length for " + calibrationFactorName); } if (swath.contains("S1")) { calibrationFactor[i] = factors[0]; } else if (swath.contains("S2")) { calibrationFactor[i] = factors[1]; } else if (swath.contains("S3")) { calibrationFactor[i] = factors[2]; } else if (swath.contains("S4")) { calibrationFactor[i] = factors[3]; } else if (swath.contains("S5")) { calibrationFactor[i] = factors[4]; } else if (swath.contains("S6")) { calibrationFactor[i] = factors[5]; } else if (swath.contains("S7")) { calibrationFactor[i] = factors[6]; } else { throw new OperatorException("Invalid swath"); } } } } catch (IOException e) { throw new OperatorException(e); } if (Double.compare(calibrationFactor[0], 0.0) == 0 && Double.compare(calibrationFactor[1], 0.0) == 0) { throw new OperatorException("Calibration factors in user provided auxiliary file are zero"); } } /** * Get calibration factors from Metadata for each band in the product. * Here it is assumed that the calibration factor values do not change in case that there are * multiple records in the Main Processing Parameters data set. */ private void getCalibrationFactorFromMetadata() { newCalibrationConstant[0] = absRoot.getAttributeDouble(AbstractMetadata.calibration_factor); if (productType.startsWith("ASA_AP")) { newCalibrationConstant[1] = absRoot.getAttributeDouble(AbstractMetadata.calibration_factor + ".2"); } } /** * Get orbit state vectors from the abstracted metadata. */ private void getOrbitStateVectors() { orbitStateVectors = AbstractMetadata.getOrbitStateVectors(absRoot); } /** * Get first line time from the abstracted metadata (in days). */ private void getFirstLineTime() { firstLineUTC = absRoot.getAttributeUTC(AbstractMetadata.first_line_time).getMJD(); // in days } /** * Get line time interval from the abstracted metadata (in days). */ private void getLineTimeInterval() { lineTimeInterval = absRoot.getAttributeDouble(AbstractMetadata.line_time_interval) / Constants.secondsInDay; // s to day } /** * Get average scene height from abstracted metadata. * * @throws Exception The exceptions. */ private void getAverageSceneHeight() throws Exception { avgSceneHeight = AbstractMetadata.getAttributeDouble(absRoot, AbstractMetadata.avg_scene_height); } /** * Get the new antenna pattern gain from the latest XCA file available. */ private void getNewAntennaPattern() { if (wideSwathProductFlag) { newRefElevationAngle = new double[5]; // reference elevation angles for 5 sub swathes newAntennaPatternWideSwath = new float[5][numOfGains]; // antenna pattern gain for 5 sub swathes getWideSwathAntennaPatternGainFromAuxData(newXCAFilePath, mdsPolar[0], numOfGains, newRefElevationAngle, newAntennaPatternWideSwath); } else { newRefElevationAngle = new double[1]; // reference elevation angle for 1 swath newAntennaPatternSingleSwath = new float[2][numOfGains]; // antenna pattern gain for 2 bands getSingleSwathAntennaPatternGainFromAuxData(newXCAFilePath, swath, mdsPolar, numOfGains, newRefElevationAngle, newAntennaPatternSingleSwath); } } /** * Get reference elevation angle and antenna pattern gain from auxiliary file for single swath product. * * @param fileName The auxiliary data file name * @param swath The swath name. * @param pol The polarizations for 2 bands. * @param numOfGains The number of gains for given swath and polarization (201). * @param refElevAngle The reference elevation angle array. * @param antPatArray The antenna pattern array. * @throws OperatorException The IO exception. */ public static void getSingleSwathAntennaPatternGainFromAuxData(String fileName, String swath, String[] pol, int numOfGains, double[] refElevAngle, float[][] antPatArray) throws OperatorException { final EnvisatAuxReader reader = new EnvisatAuxReader(); try { reader.readProduct(fileName); String swathName; if (swath.contains("S1")) { swathName = "is1"; } else if (swath.contains("S2")) { swathName = "is2"; } else if (swath.contains("S3")) { swathName = "is3_ss2"; } else if (swath.contains("S4")) { swathName = "is4_ss3"; } else if (swath.contains("S5")) { swathName = "is5_ss4"; } else if (swath.contains("S6")) { swathName = "is6_ss5"; } else if (swath.contains("S7")) { swathName = "is7"; } else { throw new OperatorException("Invalid swath"); } final String refElevAngleName = "elev_ang_" + swathName; final ProductData refElevAngleData = reader.getAuxData(refElevAngleName); refElevAngle[0] = (double) refElevAngleData.getElemFloat(); final String patternName = "pattern_" + swathName; final ProductData patternData = reader.getAuxData(patternName); final float[] pattern = ((float[]) patternData.getElems()); if (pattern.length != 804) { throw new OperatorException("Incorret array length for " + patternName); } for (int i = 0; i < 2 && pol[i] != null && pol[i].length() != 0; i++) { if (pol[i].contains("hh")) { System.arraycopy(pattern, 0, antPatArray[i], 0, numOfGains); } else if (pol[i].contains("vv")) { System.arraycopy(pattern, numOfGains, antPatArray[i], 0, numOfGains); } else if (pol[i].contains("hv")) { System.arraycopy(pattern, 2 * numOfGains, antPatArray[i], 0, numOfGains); } else if (pol[i].contains("vh")) { System.arraycopy(pattern, 3 * numOfGains, antPatArray[i], 0, numOfGains); } } } catch (IOException e) { throw new OperatorException(e); } } /** * Get reference elevation angle and antenna pattern gain from auxiliary file for wide swath product. * * @param fileName The auxiliary data file name * @param pol The polarization. * @param numOfGains The number of gains for given swath and polarization (201). * @param refElevAngle The reference elevation angle array. * @param antPatArray The antenna pattern array. * @throws OperatorException The IO exception. */ public static void getWideSwathAntennaPatternGainFromAuxData(String fileName, String pol, int numOfGains, double[] refElevAngle, float[][] antPatArray) throws OperatorException { final EnvisatAuxReader reader = new EnvisatAuxReader(); try { reader.readProduct(fileName); final String[] swathName = { "ss1", "is3_ss2", "is4_ss3", "is5_ss4", "is6_ss5" }; for (int i = 0; i < swathName.length; i++) { // read elevation angles final String refElevAngleName = "elev_ang_" + swathName[i]; final ProductData refElevAngleData = reader.getAuxData(refElevAngleName); refElevAngle[i] = (double) refElevAngleData.getElemFloat(); // read antenna pattern gains final String patternName = "pattern_" + swathName[i]; final ProductData patternData = reader.getAuxData(patternName); final float[] pattern = ((float[]) patternData.getElems()); if (pattern.length != 804) { throw new OperatorException("Incorret array length for " + patternName); } if (pol.contains("hh")) { System.arraycopy(pattern, 0, antPatArray[i], 0, numOfGains); } else if (pol.contains("vv")) { System.arraycopy(pattern, numOfGains, antPatArray[i], 0, numOfGains); } else if (pol.contains("hv")) { System.arraycopy(pattern, 2 * numOfGains, antPatArray[i], 0, numOfGains); } else if (pol.contains("vh")) { System.arraycopy(pattern, 3 * numOfGains, antPatArray[i], 0, numOfGains); } } } catch (IOException e) { throw new OperatorException(e); } } /** * Set power coefficient used in range spreading loss compensation computation for slant range images. */ private void setRangeSpreadingLossCompPower() { rangeSpreadingCompPower = 3.0; if (productType.contains("ASA_APS_1") || productType.contains("ASA_WSS")) { rangeSpreadingCompPower = 4.0; } halfRangeSpreadingCompPower = rangeSpreadingCompPower / 2.0; } /** * Update the metadata in the target product. */ public void updateTargetProductMetadata() { final MetadataElement tgtAbsRoot = AbstractMetadata.getAbstractedMetadata(targetProduct); if (applyAntennaPatternCorr) { AbstractMetadata.setAttribute(tgtAbsRoot, AbstractMetadata.ant_elev_corr_flag, 1); } if (applyRangeSpreadingCorr) { AbstractMetadata.setAttribute(tgtAbsRoot, AbstractMetadata.range_spread_comp_flag, 1); } AbstractMetadata.setAttribute(tgtAbsRoot, AbstractMetadata.abs_calibration_flag, 1); if (newXCAFileName != null) { AbstractMetadata.setAttribute(tgtAbsRoot, AbstractMetadata.external_calibration_file, newXCAFileName); } AbstractMetadata.setAttribute(tgtAbsRoot, AbstractMetadata.calibration_factor, newCalibrationConstant[0]); } /** * Called by the framework in order to compute a tile for the given target band. * <p>The default implementation throws a runtime exception with the message "not implemented".</p> * * @param targetBand The target band. * @param targetTile The current tile associated with the target band to be computed. * @param pm A progress monitor which should be used to determine computation cancelation requests. * @throws OperatorException If an error occurs during computation of the target raster. */ @Override public void computeTile(Band targetBand, Tile targetTile, ProgressMonitor pm) throws OperatorException { final Rectangle targetTileRectangle = targetTile.getRectangle(); final int x0 = targetTileRectangle.x; final int y0 = targetTileRectangle.y; final int w = targetTileRectangle.width; final int h = targetTileRectangle.height; //System.out.println("x0 = " + x0 + ", y0 = " + y0 + ", w = " + w + ", h = " + h); Band sourceBand1; Tile sourceRaster1; ProductData srcData1; ProductData srcData2 = null; final String[] srcBandNames = targetBandNameToSourceBandName.get(targetBand.getName()); if (srcBandNames.length == 1) { sourceBand1 = sourceProduct.getBand(srcBandNames[0]); sourceRaster1 = calibrationOp.getSourceTile(sourceBand1, targetTileRectangle); srcData1 = sourceRaster1.getDataBuffer(); } else { sourceBand1 = sourceProduct.getBand(srcBandNames[0]); final Band sourceBand2 = sourceProduct.getBand(srcBandNames[1]); sourceRaster1 = calibrationOp.getSourceTile(sourceBand1, targetTileRectangle); final Tile sourceRaster2 = calibrationOp.getSourceTile(sourceBand2, targetTileRectangle); srcData1 = sourceRaster1.getDataBuffer(); srcData2 = sourceRaster2.getDataBuffer(); } final Unit.UnitType tgtBandUnit = Unit.getUnitType(targetBand); final Unit.UnitType srcBandUnit = Unit.getUnitType(sourceBand1); // copy band if unit is phase if (tgtBandUnit == Unit.UnitType.PHASE) { targetTile.setRawSamples(sourceRaster1.getRawSamples()); return; } final String pol = OperatorUtils.getBandPolarization(srcBandNames[0], absRoot); int prodBand = 0; if (pol != null && mdsPolar[1] != null && mdsPolar[1].contains(pol)) { prodBand = 1; } final ProductData trgData = targetTile.getDataBuffer(); final TileIndex srcIndex = new TileIndex(sourceRaster1); final TileIndex tgtIndex = new TileIndex(targetTile); final int maxY = y0 + h; final int maxX = x0 + w; final double[] incidenceAnglesArray = new double[w]; final double[] slantRangeTimeArray = new double[w]; double[][] targetTileOldAntPat = null; // old antenna pattern gains for row pixels in a tile, in linear scale double[][] targetTileNewAntPat = null; // new antenna pattern gains for row pixels in a tile, in linear scale double[][] targetTileSlantRange = null; // slant range for pixels in a tile, in m if (applyAntennaPatternCorr) { targetTileNewAntPat = new double[h][w]; targetTileSlantRange = new double[h][w]; if (retroCalibrationFlag) { targetTileOldAntPat = new double[h][w]; } if (wideSwathProductFlag) { computeWideSwathAntennaPatternForCurrentTile(x0, y0, w, h, targetTileOldAntPat, targetTileNewAntPat, targetTileSlantRange, slantRangeTPGInterp); } else { computeSingleSwathAntennaPatternForCurrentTile(x0, y0, w, h, targetTileOldAntPat, targetTileNewAntPat, targetTileSlantRange, prodBand, slantRangeTPGInterp); } } double sigma, dn, dn2, i, q, phaseTerm = 0.0; final double theCalibrationFactor = newCalibrationConstant[prodBand]; int srcIdx, tgtIdx; for (int y = y0, yy = 0; y < maxY; ++y, ++yy) { srcIndex.calculateStride(y); tgtIndex.calculateStride(y); incidenceTPGInterp.getPixels(x0, y, w, 1, incidenceAnglesArray, pm, TiePointInterpolator.InterpMode.QUADRATIC); if (applyRangeSpreadingCorr) { slantRangeTPGInterp.getPixels(x0, y, w, 1, slantRangeTimeArray, pm, TiePointInterpolator.InterpMode.QUADRATIC); } for (int x = x0, xx = 0; x < maxX; ++x, ++xx) { srcIdx = srcIndex.getIndex(x); tgtIdx = tgtIndex.getIndex(x); if (srcBandUnit == Unit.UnitType.AMPLITUDE) { dn = srcData1.getElemDoubleAt(srcIdx); dn2 = dn * dn; } else if (srcBandUnit == Unit.UnitType.INTENSITY) { dn2 = srcData1.getElemDoubleAt(srcIdx); } else if (srcBandUnit == Unit.UnitType.REAL) { i = srcData1.getElemDoubleAt(srcIdx); q = srcData2.getElemDoubleAt(srcIdx); dn2 = i * i + q * q; if (tgtBandUnit == Unit.UnitType.REAL) { phaseTerm = i / Math.sqrt(dn2); } else if (tgtBandUnit == Unit.UnitType.IMAGINARY) { phaseTerm = q / Math.sqrt(dn2); } } else if (srcBandUnit == Unit.UnitType.INTENSITY_DB) { dn2 = FastMath.pow(10, srcData1.getElemDoubleAt(srcIdx) / 10.0); // convert dB to linear scale } else { throw new OperatorException("ASAR Calibration: unhandled unit"); } double calFactor = 1.0; if (retroCalibrationFlag) { // remove old antenna pattern gain calFactor *= targetTileOldAntPat[yy][xx]; // see Andrea's email dated Nov. 11, 2008 } // apply calibration constant and incidence angle corrections calFactor *= FastMath.sin(incidenceAnglesArray[xx] * Constants.DTOR) / theCalibrationFactor; if (applyRangeSpreadingCorr && targetTileSlantRange != null) { // apply range spreading loss compensation calFactor *= FastMath.pow(targetTileSlantRange[yy][xx] / refSlantRange800km, rangeSpreadingCompPower); } if (applyAntennaPatternCorr) { // apply antenna pattern correction calFactor /= targetTileNewAntPat[yy][xx]; // see Andrea's email dated Nov. 11, 2008 } sigma = dn2 * calFactor; if (isComplex && outputImageInComplex) { sigma = Math.sqrt(sigma) * phaseTerm; } if (outputImageScaleInDb) { // convert calibration result to dB if (sigma < underFlowFloat) { sigma = -underFlowFloat; } else { sigma = 10.0 * Math.log10(sigma); } } trgData.setElemDoubleAt(tgtIdx, sigma); } } } /** * Compute antenna pattern for the middle row of the given tile for single swath product. * Here it is assumed that the elevation angles for pixels in the same column are the same. * * @param x0 The x coordinate of the upper left point in the current tile. * @param y0 The y coordinate of the upper left point in the current tile. * @param w The width of the current tile. * @param h The height of the current tile. * @param band The band index. */ private void computeSingleSwathAntennaPatternForCurrentTile(final int x0, final int y0, final int w, final int h, final double[][] targetTileOldAntPat, final double[][] targetTileNewAntPat, final double[][] targetTileSlantRange, final int band, final TiePointInterpolator slantRangeTPGInterp) { final int yMax = y0 + h; for (int y = y0; y < yMax; y++) { final double zeroDopplerTime = firstLineUTC + y * lineTimeInterval; final double satelliteHeight = computeSatelliteHeight(zeroDopplerTime, orbitStateVectors); AbstractMetadata.SRGRCoefficientList srgrConvParam = null; if (srgrFlag) { srgrConvParam = getSRGRCoefficientsForARangeLine(zeroDopplerTime); } final int yy = y - y0; final int xMax = x0 + w; for (int x = x0; x < xMax; x++) { final int xx = x - x0; targetTileSlantRange[yy][xx] = computeSlantRange(x, y, srgrConvParam, slantRangeTPGInterp); // in m final double localEarthRadius = getEarthRadius(x, y); final double theta = computeElevationAngle(targetTileSlantRange[yy][xx], satelliteHeight, avgSceneHeight + localEarthRadius); // in degree /* double alpha = incidenceAngle.getPixelDouble(x, y); // in degree double gamma = Math.asin(targetTileSlantRange[yy][xx]*Math.sin(alpha*MathUtils.DTOR)/satelitteHeight)*MathUtils.RTOD; // in degree double theta = alpha - gamma; // in degree */ targetTileNewAntPat[yy][xx] = computeAntPatGain(theta, newRefElevationAngle[0], newAntennaPatternSingleSwath[band]); if (retroCalibrationFlag) { targetTileOldAntPat[yy][xx] = computeAntPatGain(theta, oldRefElevationAngle[0], oldAntennaPatternSingleSwath[band]); } } } } /** * Compute antenna pattern for the middle row of the given tile for wide swath product. * Here it is assumed that the elevation angles for pixels in the same column are the same. * * @param x0 The x coordinate of the upper left point in the current tile. * @param y0 The y coordinate of the upper left point in the current tile. * @param w The width of the current tile. * @param h The height of the current tile. */ private void computeWideSwathAntennaPatternForCurrentTile(final int x0, final int y0, final int w, final int h, final double[][] targetTileOldAntPat, final double[][] targetTileNewAntPat, final double[][] targetTileSlantRange, final TiePointInterpolator slantRangeTPGInterp) { final int yMax = y0 + h; for (int y = y0; y < yMax; y++) { final double zeroDopplerTime = firstLineUTC + y * lineTimeInterval; final double satelitteHeight = computeSatelliteHeight(zeroDopplerTime, orbitStateVectors); AbstractMetadata.SRGRCoefficientList srgrConvParam = null; if (srgrFlag) { srgrConvParam = getSRGRCoefficientsForARangeLine(zeroDopplerTime); } final int yy = y - y0; final int xMax = x0 + w; for (int x = x0; x < xMax; x++) { final int xx = x - x0; targetTileSlantRange[yy][xx] = computeSlantRange(x, y, srgrConvParam, slantRangeTPGInterp); // in m final double localEarthRadius = getEarthRadius(x, y); final double theta = computeElevationAngle(targetTileSlantRange[yy][xx], satelitteHeight, avgSceneHeight + localEarthRadius); // in degree int subSwathIndex = findSubSwath(theta, newRefElevationAngle); targetTileNewAntPat[yy][xx] = computeAntPatGain(theta, newRefElevationAngle[subSwathIndex], newAntennaPatternWideSwath[subSwathIndex]); if (retroCalibrationFlag) { subSwathIndex = findSubSwath(theta, oldRefElevationAngle); targetTileOldAntPat[yy][xx] = computeAntPatGain(theta, oldRefElevationAngle[subSwathIndex], oldAntennaPatternWideSwath[subSwathIndex]); } } } } /** * Get Earth radius (in m) for given pixel. * * @param x The x coordinate of the given pixel. * @param y The y coordinate of the given pixel. * @return The earth radius. */ private double getEarthRadius(final int x, final int y) { // use linear rather than quadratic interpolation for subset because quadratic interpolation // does not give accurate result in this case int i = (int) ((latMax - latitude.getPixelDouble(x, y)) / delLat + 0.5); if (i < 0) { i = 0; } else if (i >= earthRadius.length) { i = earthRadius.length - 1; } return earthRadius[i]; } /** * Find the sub swath index for given elevation angle. * * @param theta The elevation angle. * @param refElevationAngle The reference elevation array. * @return The sub swath index. */ public static int findSubSwath(double theta, double[] refElevationAngle) { // The method below finds the nearest reference elevation angle to the given elevation angle theta. // The method is equivalent to the one proposed by Romain in his email dated April 28, 2009, in which // middle point of the overlapped area of two adjacent sub swathes is used as boundary of sub swath. int idx = -1; double min = 360.0; for (int i = 0; i < refElevationAngle.length; i++) { double d = Math.abs(theta - refElevationAngle[i]); if (d < min) { min = d; idx = i; } } return idx; } /** * Compute antenna pattern gains for the given elevation angle using linear interpolation. * * @param elevAngle The elevation angle (in degree) of a given pixel. * @param refElevationAngle The reference elevation angle (in degree). * @param antPatArray The antenna pattern array. * @return The antenna pattern gain (in linear scale). */ public static double computeAntPatGain(double elevAngle, double refElevationAngle, float[] antPatArray) { final double delta = 0.05; int k0 = (int) ((elevAngle - refElevationAngle + 5.0) / delta); if (k0 < 0) { k0 = 0; } else if (k0 >= antPatArray.length - 1) { k0 = antPatArray.length - 2; } final double theta0 = refElevationAngle - 5.0 + k0 * delta; final double theta1 = theta0 + delta; final double gain0 = FastMath.pow(10, (double) antPatArray[k0] / 10.0); // convert dB to linear scale final double gain1 = FastMath.pow(10, (double) antPatArray[k0 + 1] / 10.0); final double mu = (elevAngle - theta0) / (theta1 - theta0); return Maths.interpolationLinear(gain0, gain1, mu); } //============================================================================================================ /** * Get the SRGR coefficients for given zero Doppler time. * * @param zeroDopplerTime The zero Doppler time in MJD. * @return The SRGR coefficients. */ private AbstractMetadata.SRGRCoefficientList getSRGRCoefficientsForARangeLine(final double zeroDopplerTime) { if (srgrConvParams.length == 1) { return srgrConvParams[0]; } int idx = 0; for (int i = 0; i < srgrConvParams.length && zeroDopplerTime >= srgrConvParams[i].timeMJD; i++) { idx = i; } if (idx == srgrConvParams.length - 1) { idx--; } AbstractMetadata.SRGRCoefficientList srgrConvParam = new AbstractMetadata.SRGRCoefficientList(); srgrConvParam.timeMJD = zeroDopplerTime; srgrConvParam.ground_range_origin = srgrConvParams[idx].ground_range_origin; srgrConvParam.coefficients = new double[srgrConvParams[idx].coefficients.length]; final double mu = (zeroDopplerTime - srgrConvParams[idx].timeMJD) / (srgrConvParams[idx + 1].timeMJD - srgrConvParams[idx].timeMJD); for (int i = 0; i < srgrConvParam.coefficients.length; i++) { srgrConvParam.coefficients[i] = Maths.interpolationLinear(srgrConvParams[idx].coefficients[i], srgrConvParams[idx + 1].coefficients[i], mu); } return srgrConvParam; } /** * Compute slant range for given pixel. * * @param x The x coordinate of the pixel in the source image. * @param y The y coordinate of the pixel in the source image. * @param srgrConvParam The SRGR coefficients. * @return The slant range (in meters). */ private double computeSlantRange(int x, int y, AbstractMetadata.SRGRCoefficientList srgrConvParam, final TiePointInterpolator slantRangeTPGInterp) { if (srgrFlag) { // for ground detected product, compute slant range from SRGR coefficients return Maths.computePolynomialValue(x * rangeSpacing + srgrConvParam.ground_range_origin, srgrConvParam.coefficients); } else { // for slant range product, compute slant range from slant range time final double time = slantRangeTPGInterp.getPixelDouble(x, y, TiePointInterpolator.InterpMode.QUADRATIC) / Constants.oneBillion; //convert ns to s return time * Constants.halfLightSpeed; // in m } } /** * Compute distance from satelitte to the Earth centre (in meters). * * @param zeroDopplerTime The zero Doppler time (in days). * @param orbitStateVectors The orbit state vectors. * @return The distance. */ public static double computeSatelliteHeight(final double zeroDopplerTime, final OrbitStateVector[] orbitStateVectors) { // todo should use the 3rd state vector as suggested by the doc? int idx = 0; for (int i = 0; i < orbitStateVectors.length && zeroDopplerTime >= orbitStateVectors[i].time_mjd; i++) { idx = i; } final double xPos = orbitStateVectors[idx].x_pos; final double yPos = orbitStateVectors[idx].y_pos; final double zPos = orbitStateVectors[idx].z_pos; return Math.sqrt(xPos * xPos + yPos * yPos + zPos * zPos); } /** * Compute earth radius for all range lines (in m). */ private void computeEarthRadius() { final OperatorUtils.ImageGeoBoundary imageGeoBoundary = OperatorUtils .computeImageGeoBoundary(sourceProduct); latMax = imageGeoBoundary.latMax; final double latMin = imageGeoBoundary.latMin; final double minSpacing = Math.min(rangeSpacing, azimuthSpacing); double minAbsLat; if (latMin * latMax > 0) { minAbsLat = Math.min(Math.abs(latMin), Math.abs(latMax)) * Constants.DTOR; } else { minAbsLat = 0.0; } delLat = minSpacing / Constants.MeanEarthRadius * Constants.RTOD; final double delLon = minSpacing / (Constants.MeanEarthRadius * FastMath.cos(minAbsLat)) * Constants.RTOD; delLat = Math.min(delLat, delLon); final int h = (int) ((latMax - latMin) / delLat) + 1; earthRadius = new double[h + 1]; for (int i = 0; i <= h; i++) { earthRadius[i] = computeEarthRadius((float) (latMax - i * delLat), 0.0f); } } /** * Compute Earth radius (in meters) for given pixel in source image. * * @param lat The latitude of a given pixel in source image. * @param lon The longitude of a given pixel in source image. * @return The Earth radius. */ private static double computeEarthRadius(float lat, float lon) { final double[] xyz = new double[3]; GeoUtils.geo2xyz(lat, lon, 0.0, xyz, GeoUtils.EarthModel.WGS84); return Math.sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2]); } /** * Compute elevation angle (in degree). * * @param slantRange The slant range (in meters). * @param satelliteHeight The distance from satelitte to the Earth centre (in meters). * @param sceneToEarthCentre The distance from the backscatter element to the Earth centre (in meters). * @return The elevation angle. */ public static double computeElevationAngle(final double slantRange, final double satelliteHeight, final double sceneToEarthCentre) { return FastMath .acos((slantRange * slantRange + satelliteHeight * satelliteHeight - sceneToEarthCentre * sceneToEarthCentre) / (2 * slantRange * satelliteHeight)) * Constants.RTOD; } //==================================== pixel calibration used by RD ====================================== /** * Remove the antenna pattern compensation and range spreading loss applied to the pixel. * * @param x The x coordinate of the pixel in the source image. * @param y The y coordinate of the pixel in the source image. * @param v The pixel value. * @param bandPolar The polarization of the source band. * @param bandUnit The source band unit. * @param subSwathIndex The sub swath index for current pixel for wide swath product case. * @return The pixel value with antenna pattern compensation and range spreading loss correction removed. */ public double applyRetroCalibration(int x, int y, double v, String bandPolar, final Unit.UnitType bandUnit, int[] subSwathIndex) { if (!retroCalibrationFlag) { return v; } int bandPolarIdx = 0; if (bandPolar != null && mdsPolar[1] != null && mdsPolar[1].contains(bandPolar)) { bandPolarIdx = 1; } final double zeroDopplerTime = firstLineUTC + y * lineTimeInterval; final double satelitteHeight = computeSatelliteHeight(zeroDopplerTime, orbitStateVectors); AbstractMetadata.SRGRCoefficientList srgrConvParam = null; if (srgrFlag) { srgrConvParam = getSRGRCoefficientsForARangeLine(zeroDopplerTime); } final TiePointInterpolator slantRangeTPGInterp = new TiePointInterpolator(slantRangeTime); final double slantRange = computeSlantRange(x, y, srgrConvParam, slantRangeTPGInterp); // in m final double elevationAngle = computeElevationAngle(slantRange, satelitteHeight, avgSceneHeight + getEarthRadius(x, y)); double gain = 0.0; if (wideSwathProductFlag) { gain = getAntennaPatternGain(elevationAngle, bandPolarIdx, oldRefElevationAngle, oldAntennaPatternWideSwath, true, subSwathIndex); } else { gain = computeAntPatGain(elevationAngle, oldRefElevationAngle[0], oldAntennaPatternSingleSwath[bandPolarIdx]); } if (bandUnit == Unit.UnitType.AMPLITUDE) { return v * Math.sqrt(gain) * FastMath.pow(refSlantRange800km / slantRange, halfRangeSpreadingCompPower); // amplitude } else if (bandUnit == Unit.UnitType.AMPLITUDE_DB) { return 10.0 * Math.log10(FastMath.pow(10, v / 10.0) * Math.sqrt(gain) * FastMath.pow(refSlantRange800km / slantRange, halfRangeSpreadingCompPower)); } else if (bandUnit == Unit.UnitType.INTENSITY || bandUnit == Unit.UnitType.REAL || bandUnit == Unit.UnitType.IMAGINARY) { return v * gain * FastMath.pow(refSlantRange800km / slantRange, rangeSpreadingCompPower); // intensity } else if (bandUnit == Unit.UnitType.INTENSITY_DB) { return 10.0 * Math.log10(FastMath.pow(10, v / 10.0) * gain * FastMath.pow(refSlantRange800km / slantRange, rangeSpreadingCompPower)); } else { throw new OperatorException("Unknown band unit"); } } /** * Get antenna pattern gain value for given elevation angle. * * @param elevationAngle The elevation angle (in degree). * @param bandPolar The source band polarization index. * @param refElevationAngle The reference elevation angles for different swathes or sub swathes. * @param antennaPattern The antenna pattern array. For single swath product, it contains two 201-length arrays * corresponding to the two bands of different polarizations. For wide swath product, it * contains five 201-length arrays with each for a sub swath. * @param compSubSwathIdx The boolean flag indicating if sub swath index should be computed. * @param subSwathIndex The sub swath index for current pixel for wide swath product case. * @return The antenna pattern gain value. */ private static double getAntennaPatternGain(double elevationAngle, int bandPolar, double[] refElevationAngle, float[][] antennaPattern, boolean compSubSwathIdx, int[] subSwathIndex) { if (refElevationAngle.length == 1) { // single swath return computeAntPatGain(elevationAngle, refElevationAngle[0], antennaPattern[bandPolar]); } else { // wide swath if (compSubSwathIdx || subSwathIndex[0] == INVALID_SUB_SWATH_INDEX) { subSwathIndex[0] = findSubSwath(elevationAngle, refElevationAngle); } return computeAntPatGain(elevationAngle, refElevationAngle[subSwathIndex[0]], antennaPattern[subSwathIndex[0]]); } } /** * Apply calibrations to the given point. The following calibrations are included: calibration constant, * antenna pattern compensation, range spreading loss correction and incidence angle correction. * * @param v The pixel value. * @param slantRange The slant range (in m). * @param satelliteHeight The distance from satellite to earth centre (in m). * @param sceneToEarthCentre The distance from the backscattering element position to earth centre (in m). * @param localIncidenceAngle The local incidence angle (in degrees). * @param bandPolar The source band polarization index. * @param bandUnit The source band unit. * @param subSwathIndex The sub swath index for current pixel for wide swath product case. * @return The calibrated pixel value. */ public double applyCalibration(final double v, final double rangeIndex, final double azimuthIndex, final double slantRange, final double satelliteHeight, final double sceneToEarthCentre, final double localIncidenceAngle, final String bandName, final String bandPolar, final Unit.UnitType bandUnit, int[] subSwathIndex) { int bandPolarIdx = 0; if (bandPolar != null && mdsPolar[1] != null && mdsPolar[1].contains(bandPolar)) { bandPolarIdx = 1; } double sigma = 0.0; if (bandUnit == Unit.UnitType.AMPLITUDE) { sigma = v * v; } else if (bandUnit == Unit.UnitType.AMPLITUDE_DB) { sigma = FastMath.pow(10, v / 5.0); // convert dB to linear scale, then square } else if (bandUnit == Unit.UnitType.INTENSITY || bandUnit == Unit.UnitType.REAL || bandUnit == Unit.UnitType.IMAGINARY) { sigma = v; } else if (bandUnit == Unit.UnitType.INTENSITY_DB) { sigma = FastMath.pow(10, v / 10.0); // convert dB to linear scale } else { throw new OperatorException("Unknown band unit"); } sigma *= FastMath.sin(Math.abs(localIncidenceAngle) * Constants.DTOR) / newCalibrationConstant[bandPolarIdx]; if (multilookFlag && antElevCorrFlag) { // calibration constant and incidence angle corrections only return sigma; } if (auxFile == null || !auxFile.contains(CalibrationOp.PRODUCT_AUX)) { sigma *= FastMath.pow(slantRange / refSlantRange800km, rangeSpreadingCompPower); } if (applyAntennaPatternCorr) { final double elevationAngle = computeElevationAngle(slantRange, satelliteHeight, sceneToEarthCentre); // in degrees double gain; if (wideSwathProductFlag) { if (subSwathIndex[0] == INVALID_SUB_SWATH_INDEX) { // Rem(AP+RSL)/ApplyADC Op is used final TiePointInterpolator slantRangeTPGInterp = new TiePointInterpolator(slantRangeTime); computeSubSwathIndex(rangeIndex, azimuthIndex, newRefElevationAngle, subSwathIndex, slantRangeTPGInterp); } gain = getAntennaPatternGain(elevationAngle, bandPolarIdx, newRefElevationAngle, newAntennaPatternWideSwath, false, subSwathIndex); } else { gain = computeAntPatGain(elevationAngle, newRefElevationAngle[0], newAntennaPatternSingleSwath[bandPolarIdx]); } sigma /= gain; } return sigma; } private void computeSubSwathIndex(final double rangeIndex, final double azimuthIndex, final double[] refElevationAngle, int[] subSwathIndex, final TiePointInterpolator slantRangeTPGInterp) { final int x = (int) (rangeIndex + 0.5); final int y = (int) (azimuthIndex + 0.5); final double zeroDopplerTime = firstLineUTC + y * lineTimeInterval; final double satelitteHeight = computeSatelliteHeight(zeroDopplerTime, orbitStateVectors); AbstractMetadata.SRGRCoefficientList srgrConvParam = null; if (srgrFlag) { srgrConvParam = getSRGRCoefficientsForARangeLine(zeroDopplerTime); } final double slantRange = computeSlantRange(x, y, srgrConvParam, slantRangeTPGInterp); // in m final double elevationAngle = computeElevationAngle(slantRange, satelitteHeight, avgSceneHeight + getEarthRadius(x, y)); subSwathIndex[0] = findSubSwath(elevationAngle, refElevationAngle); } public void removeFactorsForCurrentTile(Band targetBand, Tile targetTile, String srcBandName) throws OperatorException { if (!srgrFlag) { return; } final Rectangle targetTileRectangle = targetTile.getRectangle(); final int x0 = targetTileRectangle.x; final int y0 = targetTileRectangle.y; final int w = targetTileRectangle.width; final int h = targetTileRectangle.height; //System.out.println("x0 = " + x0 + ", y0 = " + y0 + ", w = " + w + ", h = " + h); final Band sourceBand = sourceProduct.getBand(srcBandName); final Tile sourceTile = calibrationOp.getSourceTile(sourceBand, targetTileRectangle); final ProductData srcData = sourceTile.getDataBuffer(); final Unit.UnitType bandUnit = Unit.getUnitType(sourceBand); final ProductData trgData = targetTile.getDataBuffer(); final String pol = OperatorUtils.getBandPolarization(targetBand.getName(), absRoot); int prodBand = 0; if (pol != null && mdsPolar[1] != null && mdsPolar[1].contains(pol)) { prodBand = 1; } final double[][] targetTileNewAntPat = new double[h][w]; final double[][] targetTileSlantRange = new double[h][w]; final double[][] targetTileOldAntPat = new double[h][w]; final TiePointInterpolator slantRangeTPGInterp = new TiePointInterpolator(slantRangeTime); if (wideSwathProductFlag) { computeWideSwathAntennaPatternForCurrentTile(x0, y0, w, h, targetTileOldAntPat, targetTileNewAntPat, targetTileSlantRange, slantRangeTPGInterp); } else { computeSingleSwathAntennaPatternForCurrentTile(x0, y0, w, h, targetTileOldAntPat, targetTileNewAntPat, targetTileSlantRange, prodBand, slantRangeTPGInterp); } final int maxY = y0 + h; final int maxX = x0 + w; double gain, slantRange, v; for (int y = y0, yy = 0; y < maxY; ++y, ++yy) { for (int x = x0, xx = 0; x < maxX; ++x, ++xx) { v = srcData.getElemDoubleAt(sourceTile.getDataBufferIndex(x, y)); gain = targetTileOldAntPat[yy][xx]; slantRange = targetTileSlantRange[yy][xx]; if (bandUnit == Unit.UnitType.AMPLITUDE) { v *= Math.sqrt(gain) * FastMath.pow(refSlantRange800km / slantRange, 0.5 * rangeSpreadingCompPower); } else if (bandUnit == Unit.UnitType.AMPLITUDE_DB) { v = FastMath.pow(10, v / 10.0) * Math.sqrt(gain) * FastMath.pow(refSlantRange800km / slantRange, 0.5 * rangeSpreadingCompPower); v = 10.0 * Math.log10(v); } else if (bandUnit == Unit.UnitType.INTENSITY) { v *= gain * FastMath.pow(refSlantRange800km / slantRange, rangeSpreadingCompPower); } else if (bandUnit == Unit.UnitType.INTENSITY_DB) { v = FastMath.pow(10, v / 10.0) * gain * FastMath.pow(refSlantRange800km / slantRange, rangeSpreadingCompPower); v = 10.0 * Math.log10(v); } else { throw new OperatorException("Unknown band unit"); } trgData.setElemDoubleAt(targetTile.getDataBufferIndex(x, y), v); } } } }