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; import com.bc.ceres.core.ProgressMonitor; import org.apache.commons.math3.util.FastMath; import org.esa.s1tbx.calibration.gpf.support.BaseCalibrator; import org.esa.s1tbx.calibration.gpf.support.Calibrator; import org.esa.snap.dataio.envisat.EnvisatAuxReader; import org.esa.snap.datamodel.AbstractMetadata; import org.esa.snap.datamodel.Unit; import org.esa.snap.eo.Constants; import org.esa.snap.framework.datamodel.Band; import org.esa.snap.framework.datamodel.MetadataAttribute; import org.esa.snap.framework.datamodel.MetadataElement; import org.esa.snap.framework.datamodel.Product; import org.esa.snap.framework.datamodel.ProductData; import org.esa.snap.framework.datamodel.RasterDataNode; import org.esa.snap.framework.datamodel.TiePointGrid; import org.esa.snap.framework.gpf.Operator; import org.esa.snap.framework.gpf.OperatorException; import org.esa.snap.framework.gpf.Tile; import org.esa.snap.gpf.OperatorUtils; import org.esa.snap.util.ResourceInstaller; import javax.media.jai.BorderExtender; import javax.media.jai.JAI; import javax.media.jai.PlanarImage; import javax.media.jai.RasterFactory; import javax.media.jai.operator.SubsampleAverageDescriptor; import java.awt.Point; import java.awt.Rectangle; import java.awt.RenderingHints; import java.awt.image.BufferedImage; import java.awt.image.ColorModel; import java.awt.image.DataBuffer; import java.awt.image.DataBufferDouble; import java.awt.image.Raster; import java.awt.image.RenderedImage; import java.awt.image.SampleModel; import java.awt.image.WritableRaster; import java.awt.image.renderable.ParameterBlock; import java.io.BufferedReader; import java.io.File; import java.io.FileInputStream; import java.io.FileNotFoundException; import java.io.IOException; import java.io.InputStreamReader; import java.nio.file.Path; import java.text.ParseException; import java.text.SimpleDateFormat; import java.util.Date; import java.util.Hashtable; import java.util.StringTokenizer; import java.util.TimeZone; /** * The objective of SAR calibration is to provide imagery in which the pixel values can be directly related * to the radar backscatter of the scene. Though uncalibrated SAR imagery is sufficient for qualitative use, * calibrated SAR images are essential to quantitative use of SAR data. Typically processing of SAR data to * produce level 1 images does not include geometric and radiometric corrections and significant radiometric * and geometric bias remains. Corrections to the radiometric and geometric bias will therefore be applied * to SAR images so that the pixel values of SAR images are a true representation of the reflecting surface's * radar backscatter. Calibration is necessary to compare SAR images acquired with different sensors as well * as when acquired from the same sensor but at different times, in different modes, or processed by different * processors. * <p> * This operator applies the following calibrations to ERS-1 and ERS-2 SAR data products generated by different * processors at ESA/ESRIN and at different Processing and Archiving Facilities such as the German PAF(D-PAF), * the Italian PAF (I-PAF) and the United-Kingdom PAF (UK-PAF): * <p> * 1. calibration constant correction * 2. antenna elevation pattern correction * 3. replica pulse power variations correction * 4. analogue to digital converter non-linearity correction */ public final class ERSCalibrator extends BaseCalibrator implements Calibrator { private String pafID; // processing facility identifier private String psID; // processing system identifier private String pvID; // processing version identifier private String extXCAFileName; private int sourceImageWidth; private int sourceImageHeight; private int windowWidth; private int windowHeight; private int blockWidth; private int blockHeight; private boolean applyAntennaPatternCorrection = false; private boolean applyRangeSpreadingLossCorrection = false; private boolean applyReplicaPowerCorrection = false; private boolean applyADCSaturationCorrection = false; private boolean isERS1Mission = false; private boolean isCEOSFormat = false; private boolean isAntPattAvailable = false; private boolean adcHasBeenTestedFlag = false; private boolean antennaPatternCorrectionFlag = false; private boolean rangeSpreadingLossCompFlag = false; private boolean useExtXCAFile = false; private boolean multilookFlag = false; private double rangeSpacing; // m private double azimuthSpacing; // m private double calibrationConstant; private double sceneCentreLatitude; // in degree private double replicaPulseVariationsCorrectionFactor; private double elevationAngle; // elevation angle for given swath IS2 VV, in degree private Date processingTime; private Date acquisitionTime; private Date time19910801; // = 681004800.0; 01-Aug-1991 00:00:00.000, in s private Date time19920401; // = 702086400.0; 01-Apr-1992 00:00:00.000, in s private Date time19920414; // = 703209600.0; 14-Apr-1992 00:00:00.000, in s private Date time19920901; // = 715305600.0; 01-Sep-1992 00:00:00.000, in s private Date time19930408; // = 734227200.0; 08-Apr-1993 00:00:00.000, in s private Date time19930628; // = 741225600.0; 28-Jun-1993 00:00:00.000, in s private Date time19941207; // = 786758400.0; 07-Dec-1994 00:00:00.000, in s private Date time19950317; // = 795398400.0; 17-Mar-1995 00:00:00.000, in s private Date time19950713; // = 805593600.0; 13-Jul-1995 00:00:00.000, in s private Date time19950716; // = 805852800.0; 16-Jul-1995 00:00:00.000, in s private Date time19951016; // = 813801600.0; 16-Oct-1995 00:00:00.000, in s private Date time19970120; // = 853718400.0; 20-Jan-1997 00:00:00.000, in s private Date time19980224; // = 888278400.0; 24-Feb-1998 00:00:00.000, in s private Date time20040904; // = 1094292254.0; 04-Sep-2004 10:04:14.000, in s private Date time20041014; // = 1097764631.0; 14-Oct-2004 14:37:11.000, in s private double[] incidenceAngles = null; // for a complete range line, in radian private double[] lookAngles = null; // for a complete range line, in radian private double[] rangeSpreadingLoss = null; // for a complete range line private double[] antennaPatternCorrFactor = null; // for a range line in current tile, in linear scale private double[] antennaPatternGain = null; // used in ADC, for a range line in current tile, in linear scale private double[][] appendixF1 = null; // ERS-1 SAR ADC Power Loss Correction Look-up Table, in dB private double[][] appendixF2 = null; // ERS-2 SAR ADC Power Loss Correction Look-up Table, in dB private double[][] appendixG1 = null; // initial ERS-1 SAR antenna pattern gain, in dB private double[][] appendixG2 = null; // improved ERS-1 SAR antenna pattern gain, in dB private double[][] appendixG3 = null; // ERS-2 SAR antenna pattern gain, in dB private double[][] appendixH = null; // UK-PAF elevation antenna pattern correction, in dB private float[] antPatForPGS = null; // antenna pattern for swath IS2 VV for PGS product, in dB // parameters used for PGS-ENVISAT calibration private int numMPPRecords; private static final double referenceIncidenceAngle = 23.0 * Constants.DTOR; // radian private static final double relativeLookAngle = 20.355; // degree private static final double aGEM6 = 6378144; // GEM6: equatorial Earth radius in m (for VMP CEOS) private static final double bGEM6 = 6356759; // GEM6: polar Earth radius in m (for VMP CEOS) private static final double aWGS84 = Constants.semiMajorAxis; // WGS 84: equatorial Earth radius in m (for PGS CEOS) private static final double bWGS84 = Constants.semiMinorAxis; // WGS 84: polar Earth radius in m (for PGS CEOS) private static final double referenceSlantRange = 847000; // m private static final double windowDimInRange = 15000.0; // m private static final double windowDimInAzimuth = 5000.0; // m private static final double downSampleBlockSize = 100.0; // m private static final double ers1ApplyADCThreshold = -7.0; // dB private static final double ers2ApplyADCThreshold = -2.0; // dB private static final String D_PAF = "D-PAF"; private static final String I_PAF = "I-PAF"; private static final String UK_PAF = "UK-PAF"; private static final String ESRIN = "ES"; private static final String VMP = "VMP"; /** * Default constructor. The graph processing framework * requires that an operator has a default constructor. */ public ERSCalibrator() { } /** * Set external auxiliary file. */ @Override public void setExternalAuxFile(final File file) throws OperatorException { if (file != null) { throw new OperatorException("No external auxiliary file should be selected for ERS product"); } } /** * Set auxiliary file flag. */ @Override public void setAuxFileFlag(String file) { } /** */ 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; sourceImageWidth = sourceProduct.getSceneRasterWidth(); sourceImageHeight = sourceProduct.getSceneRasterHeight(); //System.out.println("sourceImageWidth = " + sourceImageWidth + ", sourceImageHeight = " + sourceImageHeight); absRoot = AbstractMetadata.getAbstractedMetadata(sourceProduct); origMetadataRoot = AbstractMetadata.getOriginalProductMetadata(sourceProduct); getImportantTimes(); getMissionType(); // abs getSampleType(); // abs getProductType(); // abs getCalibrationFlags(); // abs getMultilookFlag(); // abs getPixelSpacings(); // abs getProcessingTime(); // abs getProcessingSystemID(); // abs getProductAcquisitionTime(); // abs if (isCEOSFormat) { // CEOS getProcessingFacilityIDFromCEOS(); // ERS getProcessingVersionID(); // ERS getSceneCentreLatitude(); // ERS getCalibrationConstantFromCEOS(); // ERS } else { // ENVISAT getProcessingFacilityIDFromENVISAT(); // MPH getNumOfRecordsInMainProcParam(); // DSR.3 getCalibrationConstantFromENVISAT(); // MPP } if (absRoot.getAttribute("retro-calibration performed flag") != null) { // product from removeAntennaPatternOp setCalibrationFlags(); } else { // normal product setCorrectionFlags(); } if (applyADCSaturationCorrection) { prepareForADCCorrection(); } if (applyAntennaPatternCorrection || applyADCSaturationCorrection) { getAntennaPatternFromFile(); } if (applyReplicaPowerCorrection || applyADCSaturationCorrection) { computeReplicaPulseVariationsCorrectionFactor(); } if (isCEOSFormat) { // CEOS computeIncidenceAnglesLookAnglesRangeSpreadingLossForCEOS(); } else { // ENVISAT computeIncidenceAnglesLookAnglesRangeSpreadingLossForENVISAT(); } //computeIncidenceAnglesLookAnglesRangeSpreadingLoss(); // common to CEOS and ENVISAT if (mustUpdateMetadata) { updateTargetProductMetadata(); } } catch (Exception e) { throw new OperatorException(e); } } /** * 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 org.esa.snap.framework.gpf.OperatorException If an error occurs during computation of the target raster. */ @Override public void computeTile(Band targetBand, Tile targetTile, ProgressMonitor pm) throws OperatorException { try { 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 ProductData trgData = targetTile.getDataBuffer(); Band sourceBand1 = null; Band sourceBand2 = null; Tile sourceRaster1 = null; Tile sourceRaster2 = null; ProductData srcData1 = null; ProductData srcData2 = null; final String[] srcBandNames = targetBandNameToSourceBandName.get(targetBand.getName()); if (srcBandNames.length == 1) { sourceBand1 = sourceProduct.getBand(srcBandNames[0]); sourceRaster1 = getSourceTile(sourceBand1, targetTileRectangle); srcData1 = sourceRaster1.getDataBuffer(); } else { sourceBand1 = sourceProduct.getBand(srcBandNames[0]); sourceBand2 = sourceProduct.getBand(srcBandNames[1]); sourceRaster1 = getSourceTile(sourceBand1, targetTileRectangle); sourceRaster2 = getSourceTile(sourceBand2, targetTileRectangle); srcData1 = sourceRaster1.getDataBuffer(); srcData2 = sourceRaster2.getDataBuffer(); } final Unit.UnitType bandUnit = Unit.getUnitType(sourceBand1); // copy band if unit is phase if (bandUnit == Unit.UnitType.PHASE) { targetTile.setRawSamples(sourceRaster1.getRawSamples()); return; } if (applyAntennaPatternCorrection && !isAntPattAvailable) { computeAntennaPatternCorrectionFactors(0, sourceImageWidth); } if (applyADCSaturationCorrection && !adcHasBeenTestedFlag) { testADC(sourceBand1, sourceBand2, bandUnit); } boolean applyADCSaturationCorrectionToCurrentTile = false; if (applyADCSaturationCorrection && h >= blockHeight && w >= blockWidth) { applyADCSaturationCorrectionToCurrentTile = true; } double[][] adcPowerLoss = null; if (applyADCSaturationCorrectionToCurrentTile) { adcPowerLoss = computeADCPowerLossValuesForCurrentTile(sourceBand1, sourceBand2, x0, y0, w, h, bandUnit); } final double k = calibrationConstant * FastMath.sin(referenceIncidenceAngle); double sigma, dn, i, q; int index; int adcJ = 0; for (int x = x0; x < x0 + w; x++) { final double sinIncidenceAngleByK = FastMath.sin(incidenceAngles[x]) / k; if (applyADCSaturationCorrectionToCurrentTile) { adcJ = Math.min(((x - x0) / blockWidth), adcPowerLoss[0].length - 1); } for (int y = y0; y < y0 + h; y++) { index = sourceRaster1.getDataBufferIndex(x, y); if (bandUnit == Unit.UnitType.AMPLITUDE) { dn = srcData1.getElemDoubleAt(index); sigma = dn * dn; } else if (bandUnit == Unit.UnitType.INTENSITY) { sigma = srcData1.getElemDoubleAt(index); } else { // COMPLEX i = srcData1.getElemDoubleAt(index); q = srcData2.getElemDoubleAt(index); sigma = i * i + q * q; } sigma *= sinIncidenceAngleByK; if (applyAntennaPatternCorrection) { sigma *= antennaPatternCorrFactor[x]; } if (applyRangeSpreadingLossCorrection) { sigma *= rangeSpreadingLoss[x]; } if (applyReplicaPowerCorrection) { sigma *= replicaPulseVariationsCorrectionFactor; } if (applyADCSaturationCorrectionToCurrentTile) { final int adcI = Math.min(((y - y0) / blockHeight), adcPowerLoss.length - 1); sigma *= adcPowerLoss[adcI][adcJ]; } if (outputImageScaleInDb) { // convert calibration result to dB if (sigma < underFlowFloat) { sigma = -underFlowFloat; } else { sigma = 10.0 * Math.log10(sigma); } } trgData.setElemDoubleAt(targetTile.getDataBufferIndex(x, y), sigma); } } } catch (Throwable e) { OperatorUtils.catchOperatorException("ERSCalibrator", e); } } private synchronized void testADC(final Band sourceBand1, final Band sourceBand2, final Unit.UnitType bandUnit) { if (adcHasBeenTestedFlag) return; if (!isADCNeeded(sourceBand1, sourceBand2, bandUnit)) { applyADCSaturationCorrection = false; } // if (applyADCSaturationCorrection && antennaPatternCorrectionFlag) { if (antennaPatternCorrectionFlag) { computeAntennaPatternGain(0, sourceImageWidth); } adcHasBeenTestedFlag = true; } /** * Compute some important times (in second). The time zone of GMT is always assumed. */ private void getImportantTimes() { final SimpleDateFormat dateformat = new SimpleDateFormat("dd-MMM-yyyy HH:mm:ss.SSS"); try { time19910801 = dateformat.parse("01-Aug-1991 00:00:00.000"); time19920401 = dateformat.parse("01-Apr-1992 00:00:00.000"); time19920414 = dateformat.parse("14-Apr-1992 00:00:00.000"); time19920901 = dateformat.parse("01-Sep-1992 00:00:00.000"); time19930408 = dateformat.parse("08-Apr-1993 00:00:00.000"); time19930628 = dateformat.parse("28-Jun-1993 00:00:00.000"); time19941207 = dateformat.parse("07-Dec-1994 00:00:00.000"); time19950317 = dateformat.parse("17-Mar-1995 00:00:00.000"); time19950713 = dateformat.parse("13-Jul-1995 00:00:00.000"); time19950716 = dateformat.parse("16-Jul-1995 00:00:00.000"); time19951016 = dateformat.parse("16-Oct-1995 00:00:00.000"); time19970120 = dateformat.parse("20-Jan-1997 00:00:00.000"); time19980224 = dateformat.parse("24-Feb-1998 00:00:00.000"); time20040904 = dateformat.parse("04-Sep-2004 10:04:14.000"); time20041014 = dateformat.parse("14-Oct-2004 14:37:11.000"); } catch (ParseException e) { throw new OperatorException(e); } /* System.out.format("time19910801 = %13.3f%n", time19910801); System.out.format("time19920401 = %13.3f%n", time19920401); System.out.format("time19920414 = %13.3f%n", time19920414); System.out.format("time19920901 = %13.3f%n", time19920901); System.out.format("time19930408 = %13.3f%n", time19930408); System.out.format("time19930628 = %13.3f%n", time19930628); System.out.format("time19941207 = %13.3f%n", time19941207); System.out.format("time19950317 = %13.3f%n", time19950317); System.out.format("time19950713 = %13.3f%n", time19950713); System.out.format("time19950716 = %13.3f%n", time19950716); System.out.format("time19951016 = %13.3f%n", time19951016); System.out.format("time19970120 = %13.3f%n", time19970120); System.out.format("time19980224 = %13.3f%n", time19980224); System.out.format("time20040904 = %13.3f%n", time20040904); System.out.format("time20041014 = %13.3f%n", time20041014); */ } /** * Convert UTC time (dd-MMM-yyyy HH:mm:ss.sss) to seconds. The GMT time zone is always assumed. * * @param utcTime The UTC time string. * @return The time in second. */ private static double convertUTCTimes(final String utcTime) { double timeInSecond = 0.0; SimpleDateFormat simpledateformat; if (utcTime.length() == 23) { simpledateformat = new SimpleDateFormat("dd-MM-yyyy HH:mm:ss.SSS"); } else if (utcTime.length() == 24) { simpledateformat = new SimpleDateFormat("dd-MMM-yyyy HH:mm:ss.SSS"); } else { throw new OperatorException("Incorrect UTC time string"); } simpledateformat.setTimeZone(TimeZone.getTimeZone("UTC")); try { timeInSecond = (double) simpledateformat.parse(utcTime).getTime() / 1000.0; // ms to s } catch (ParseException e) { throw new OperatorException(e); } return timeInSecond; } /** * Get the mission type. */ private void getMissionType() { final String missionType = absRoot.getAttributeString(AbstractMetadata.MISSION); if (!missionType.equals("ERS1") && !missionType.equals("ERS2")) { throw new OperatorException(missionType + " is not a valid mission for ERS Calibration"); } if (missionType.equals("ERS1")) isERS1Mission = true; //System.out.println("Mission type is " + missionType); } /** * Get Product type. */ private void getProductType() { final String productType = absRoot.getAttributeString(AbstractMetadata.PRODUCT_TYPE); if (productType.contains("ERS")) { isCEOSFormat = true; } else if (productType.contains("SAR")) { isCEOSFormat = false; } else { throw new OperatorException("Invalid product type: " + productType); } //System.out.println("Product type is " + productType); } /** * Get the processing facility identifier. */ private void getProcessingFacilityIDFromCEOS() { // Field 81 in PRI Data Set Summary Record final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Scene Parameters"); if (facility == null) { throw new OperatorException("Scene Parameters not found"); } final MetadataAttribute attr = facility.getAttribute("Processing facility identifier"); if (attr == null) { throw new OperatorException("Processing facility identifier not found"); } pafID = attr.getData().getElemString(); /* if (!pafID.contains(ESRIN) && !pafID.contains(D_PAF) && !pafID.contains(I_PAF) && !pafID.contains(UK_PAF)) { throw new OperatorException("Invalid processing facility identifier: " + pafID); } */ //System.out.println("Processing facility identifier is " + pafID); } /** * Get the processing system identifier. */ private void getProcessingSystemID() { psID = absRoot.getAttributeString(AbstractMetadata.ProcessingSystemIdentifier); /* if (isCEOSFormat) { if (!psID.contains("VMP") && !psID.contains("PGS")) { throw new OperatorException("Invalid processing system identifier: " + psID); } } else { // ENVISAT psID = "PGS"; } */ //System.out.println("Processing system identifier is " + psID); } /** * Get the processing version identifier. */ private void getProcessingVersionID() { // Field 83 in PRI Data Set Summary Record final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Scene Parameters"); if (facility == null) { throw new OperatorException("Scene Parameters not found"); } final MetadataAttribute attr = facility.getAttribute("Processing version identifier"); if (attr == null) { throw new OperatorException("Processing version identifier not found"); } pvID = attr.getData().getElemString(); //System.out.println("Processing version identifier is " + pvID); } /** * Get the product acquisition time (in second). */ private void getProductAcquisitionTime() { try { final ProductData.UTC acqTimeUTC = AbstractMetadata .parseUTC(absRoot.getAttributeString(AbstractMetadata.first_line_time)); acquisitionTime = acqTimeUTC.getAsDate(); //System.out.println("The acquisition time is " + acquisitionTime); } catch (Exception e) { throw new OperatorException(e.getMessage()); } } /** * Get the processing time (in second). */ private void getProcessingTime() { try { final ProductData.UTC procTimeUTC = AbstractMetadata .parseUTC(absRoot.getAttributeString(AbstractMetadata.PROC_TIME)); processingTime = procTimeUTC.getAsDate(); //System.out.println("The processing time is " + processingTime); } catch (Exception e) { throw new OperatorException(e.getMessage()); } } /** * Get the range and azimuth spacings (in meter). * * @throws Exception The exceptions. */ private void getPixelSpacings() throws Exception { rangeSpacing = AbstractMetadata.getAttributeDouble(absRoot, AbstractMetadata.range_spacing); azimuthSpacing = AbstractMetadata.getAttributeDouble(absRoot, AbstractMetadata.azimuth_spacing); //System.out.println("Range spacing is " + rangeSpacing); //System.out.println("Azimuth spacing is " + azimuthSpacing); } /** * Get the antenna pattern correction flag and range spreading loss flag. * * @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"); } antennaPatternCorrectionFlag = AbstractMetadata.getAttributeBoolean(absRoot, AbstractMetadata.ant_elev_corr_flag); rangeSpreadingLossCompFlag = AbstractMetadata.getAttributeBoolean(absRoot, AbstractMetadata.range_spread_comp_flag); //System.out.println("Antenna pattern correction flag is " + antennaPatternCorrectionFlag); //System.out.println("Range spreding loss compensation flag is " + rangeSpreadingLossCompFlag); } /** * Get multilook flag from the abstracted metadata. * * @throws Exception The exceptions. */ private void getMultilookFlag() throws Exception { multilookFlag = AbstractMetadata.getAttributeBoolean(absRoot, AbstractMetadata.multilook_flag); } /** * Set correction flags in case source product is from removeAntennaPatternOp. */ private void setCalibrationFlags() { applyAntennaPatternCorrection = true; applyRangeSpreadingLossCorrection = true; applyReplicaPowerCorrection = true; applyADCSaturationCorrection = false; } /** * Set correction flags. */ private void setCorrectionFlags() { applyAntennaPatternCorrection = false; applyRangeSpreadingLossCorrection = false; applyReplicaPowerCorrection = false; applyADCSaturationCorrection = false; if (isERS1Mission) { // ERS-1 if (!isComplex) { // detected if (psID.contains(VMP) && processingTime.compareTo(time19910801) >= 0 && processingTime.compareTo(time19950716) < 0) { applyAntennaPatternCorrection = true; // because the antenna gain applied before was the initial } // one and should be replaced by the improved one. applyReplicaPowerCorrection = true; applyADCSaturationCorrection = true; } else { // complex if (!antennaPatternCorrectionFlag) { applyAntennaPatternCorrection = true; } if (!rangeSpreadingLossCompFlag) { applyRangeSpreadingLossCorrection = true; } applyReplicaPowerCorrection = true; applyADCSaturationCorrection = true; } } else { // ERS-2 if (!isComplex) { // detected applyADCSaturationCorrection = true; } else { // complex if (!antennaPatternCorrectionFlag) { applyAntennaPatternCorrection = true; } if (!rangeSpreadingLossCompFlag) { applyRangeSpreadingLossCorrection = true; } applyADCSaturationCorrection = true; } } if (applyADCSaturationCorrection) { adcHasBeenTestedFlag = false; } } /** * Get the calibration constant. */ private void getCalibrationConstantFromCEOS() { // Field 62 in PRI Facility Related Data Record final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Facility Related"); if (facility == null) { throw new OperatorException("Facility Related not found"); } final MetadataAttribute calibrationConstantAttr = facility.getAttribute("Absolute calibration constant K"); if (calibrationConstantAttr == null) { throw new OperatorException("Absolute calibration constant K not found"); } calibrationConstant = calibrationConstantAttr.getData().getElemFloat(); //System.out.println("Calibration constant is " + calibrationConstant); } /** * Get the scene centre latitude (in degree). */ private void getSceneCentreLatitude() { // Field 13 in PRI Data Set Summary Record (in degree) final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Scene Parameters"); if (facility == null) { throw new OperatorException("Scene Parameters not found"); } final MetadataAttribute attr = facility.getAttribute("scene centre geodetic latitude"); if (attr == null) { throw new OperatorException("Scene centre geodetic latitude not found"); } sceneCentreLatitude = attr.getData().getElemFloat(); //System.out.println("Scene centre geodetic latitude is " + sceneCentreLatitude); } /** * Read ERS-1 or ERS-2 antenna pattern from files. */ private void getAntennaPatternFromFile() { if (isERS1Mission) { // ERS-1 if (psID.contains(VMP)) { if (processingTime.compareTo(time19920901) >= 0 && processingTime.compareTo(time19950716) < 0 && (pafID.contains(ESRIN) || pafID.contains(D_PAF) || pafID.contains(I_PAF) || pafID.contains(UK_PAF))) { getInitialERS1ElevAntPat(); } if (processingTime.compareTo(time19910801) >= 0) { // && (pafID.contains(ESRIN) || pafID.contains(D_PAF) || pafID.contains(I_PAF) || pafID.contains(UK_PAF))) { getImprovedERS1ElevAntPat(); } if (processingTime.compareTo(time19920901) >= 0 && processingTime.compareTo(time19930408) <= 0 && pafID.contains(UK_PAF)) { getUKPAFElevAntPatCor(); } } else { // PGS (CEOS or ENVISAT) // get antenna pattern from ERS1 XCA file final Path moduleBasePath = ResourceInstaller.findModuleCodeBasePath(this.getClass()); final Path xcaFile = moduleBasePath.resolve("org/esa/s1tbx/auxdata/ers/" + "ER1_XCA_AXNXXX20050321_000000_19910101_000000_20100101_000000.zip"); getAntennaPatternGainFromAuxData(xcaFile.toFile()); extXCAFileName = "ER1_XCA_AXNXXX20050321_000000_19910101_000000_20100101_000000.zip"; useExtXCAFile = true; } } else { // ERS-2 if (psID.contains(VMP)) { getERS2ElevAntPat(); } else { // PGS (CEOS or ENVISAT) // get antenna pattern from ERS2 XCA file final Path moduleBasePath = ResourceInstaller.findModuleCodeBasePath(this.getClass()); final Path xcaFile = moduleBasePath.resolve("org/esa/s1tbx/auxdata/ers/" + "ER2_XCA_AXNXXX20050321_000000_19950101_000000_20100101_000000.zip"); getAntennaPatternGainFromAuxData(xcaFile.toFile()); extXCAFileName = "ER2_XCA_AXNXXX20050321_000000_19950101_000000_20100101_000000.zip"; useExtXCAFile = true; } } // Note: For both PGS CEOS and PGS ENVISAT (ERS-1 or ERS-2) use XCA files. See Andrea's email // dated Nov. 5, 2008. } /** * Get the initial ERS-1 SAR antenna pattern (Appendix G1). */ private void getInitialERS1ElevAntPat() { // The following time checking is ignored. See Andrea's email dated Sep. 15, 2008. //if (processingTime.compareTo(time19950716) >= 0) { // throw new OperatorException("Processing date must be earlier than 16-July-1995"); //} final String fileName = "Appendix_G1.txt"; appendixG1 = readFile(getERSAuxFile(fileName).toFile()); } /** * Get the improved ERS-1 SAR antenna pattern (Appendix G2). */ private void getImprovedERS1ElevAntPat() { String fileName = ""; if (processingTime.compareTo(time19950716) >= 0 && pvID.compareTo("6.8") < 0) { fileName = "Appendix_G2_b.txt"; } else if (pvID.compareTo("6.8") >= 0) { fileName = "Appendix_G2_c.txt"; } else { throw new OperatorException("The operator does not support VMP product processed before 1995-07-16"); } appendixG2 = readFile(getERSAuxFile(fileName).toFile()); } /** * Get the UK-PAF ERS-1 antenna pattern (Appendix H). */ private void getUKPAFElevAntPatCor() { String fileName = ""; if (acquisitionTime.compareTo(time19920401) <= 0) { //orbit repeat period: 3 days; fileName = "Appendix_H_1.txt"; } else if (acquisitionTime.compareTo(time19920414) >= 0 && acquisitionTime.compareTo(time19930408) <= 0) { //orbit repeat period: 35 days; fileName = "Appendix_H_2.txt"; } else { throw new OperatorException("Incorrect acquisition date"); } appendixH = readFile(getERSAuxFile(fileName).toFile()); } /** * Get the ERS-2 antenna pattern (Appendix G3). */ private void getERS2ElevAntPat() { String fileName = ""; if (pvID.compareTo("6.8") < 0) { fileName = "Appendix_G3_b.txt"; } else if (pvID.compareTo("6.8") >= 0) { fileName = "Appendix_G3_c.txt"; } appendixG3 = readFile(getERSAuxFile(fileName).toFile()); } private Path getERSAuxFile(final String fileName) { final Path moduleBasePath = ResourceInstaller.findModuleCodeBasePath(this.getClass()); return moduleBasePath.resolve("org/esa/s1tbx/auxdata/ers/" + fileName); } /** * Read antenna pattern data from file and save them in a 2D array. * * @param file The file * @return array The 2D array holding antenna pattern data */ private static double[][] readFile(final File file) { // get reader FileInputStream stream; try { stream = new FileInputStream(file); } catch (FileNotFoundException e) { throw new OperatorException("File not found: " + file); } final BufferedReader reader = new BufferedReader(new InputStreamReader(stream)); // read data from file and save them in 2-D array String line = ""; StringTokenizer st; double[][] array; int rowIdx = 0; try { // get the 1st line if ((line = reader.readLine()) == null) { throw new OperatorException("Empty file: " + file); } st = new StringTokenizer(line); if (st.countTokens() != 2) { throw new OperatorException("Incorrect file format: " + file); } final int numRows = Integer.parseInt(st.nextToken()); final int numCols = Integer.parseInt(st.nextToken()); array = new double[numRows][numCols]; // get the rest numRows lines while ((line = reader.readLine()) != null) { st = new StringTokenizer(line); if (st.countTokens() != numCols) { throw new OperatorException("Incorrect file format: " + file); } for (int j = 0; j < numCols; j++) { array[rowIdx][j] = Double.parseDouble(st.nextToken()); } rowIdx++; } if (numRows != rowIdx) { throw new OperatorException("Incorrect number of lines in file: " + file); } reader.close(); stream.close(); } catch (IOException e) { throw new OperatorException(e); } return array; } /** * Obtain from auxiliary data the elevation angles for given swath and the antenna elevation * pattern gains for the swath and the polarization of the product. * * @param file The auxiliary data file. * @throws OperatorException The exceptions. */ private void getAntennaPatternGainFromAuxData(final File file) throws OperatorException { // Note: ERS-1/2 predate swath-selection, in the Envisat-world all ERS image-data is IS2 VV. // See Andrea and Marcus' email dated Nov. 5, 2008. final EnvisatAuxReader reader = new EnvisatAuxReader(); try { reader.readProduct(file); final ProductData elevAngleData = reader.getAuxData("elev_ang_is2"); elevationAngle = (double) elevAngleData.getElemFloat(); //System.out.println("elevation angle is " + elevationAngle); //System.out.println(); final ProductData patData = reader.getAuxData("pattern_is2"); final float[] pattern = ((float[]) patData.getElems()); if (pattern.length != 804) { throw new OperatorException("Incorrect array length for pattern_is2"); } final int numOfGains = 201; antPatForPGS = new float[numOfGains]; System.arraycopy(pattern, numOfGains, antPatForPGS, 0, numOfGains); // polarization = VV /* for (float val : antPatForPGS) { System.out.print(val + ", "); } System.out.println(); */ } catch (IOException e) { throw new OperatorException(e); } } /** * Get chirp average density image. * * @return The chirp average density image. */ private static double getChirpAverageDensityImage() { // Byte location 3449 in "ESA reserved" record in Facility Data Record, PCS type // need to read chirp average density image (see Appendix D3 on p27) return 267.20; } /** * Compute replica pulse variations correction factor. */ private void computeReplicaPulseVariationsCorrectionFactor() { final double replicaPulsePower; if (isCEOSFormat) { replicaPulsePower = getReplicaPulsePowerForCEOS(); } else { // ENVISAT replicaPulsePower = getReplicaPulsePowerForENVISAT(); } if (Double.compare(replicaPulsePower, -9999999.9999999) == 0 || Double.compare(replicaPulsePower, 0.0) == 0) { replicaPulseVariationsCorrectionFactor = 1.0; } else { if (isERS1Mission) { replicaPulseVariationsCorrectionFactor = replicaPulsePower / 205229.0; } else { replicaPulseVariationsCorrectionFactor = replicaPulsePower / 156000.0; } } /* if (isERS1Mission) { if (pafID.contains(D_PAF) || pafID.contains(I_PAF) || pafID.contains(UK_PAF)) { if (Double.compare(replicaPulsePower, -9999999.9999999) == 0) { if (pafID.contains(D_PAF)) { replicaPulseVariationsCorrectionFactor = getChirpAverageDensityImage() / 267.20; } else { throw new OperatorException("Replica pulse power is not available"); } } else { replicaPulseVariationsCorrectionFactor = replicaPulsePower / 205229.0; } } else if (pafID.contains(ESRIN)) { replicaPulseVariationsCorrectionFactor = getChirpAverageDensityImage() / 267.20; } } else { // ERS-2 if (pafID.contains(D_PAF) || pafID.contains(I_PAF) || pafID.contains(UK_PAF) || pafID.contains(ESRIN)) { if (Double.compare(replicaPulsePower, -9999999.9999999) == 0) { throw new OperatorException("Replica pulse power is not available"); } replicaPulseVariationsCorrectionFactor = replicaPulsePower / 156000.0; } } */ } /** * Get the replica pulse power for CEOS product. * * @return The replica pulse power */ private double getReplicaPulsePowerForCEOS() { // Field 55 in PRI Facility Related Data Record final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Facility Related"); if (facility == null) { throw new OperatorException("Facility Related not found"); } final MetadataAttribute attr = facility.getAttribute("Replica pulse power"); if (attr == null) { throw new OperatorException("Replica pulse power"); } return attr.getData().getElemFloat(); } /** * Compute incidence angles (in radian), look angles (in radian) and range spreading loss * for pixels in a complete range line. */ private void computeIncidenceAnglesLookAnglesRangeSpreadingLossForCEOS() { incidenceAngles = new double[sourceImageWidth]; lookAngles = new double[sourceImageWidth]; rangeSpreadingLoss = new double[sourceImageWidth]; double a; double b; if (psID.contains(VMP)) { a = aGEM6; b = bGEM6; } else { // PGS a = aWGS84; b = bWGS84; } final double lambda = sceneCentreLatitude * Constants.DTOR; final double alpha1 = getIncidenceAngleAtFirstRangePixel() * Constants.DTOR; final double cos2 = FastMath.pow(FastMath.cos(lambda), 2.0); final double sin2 = FastMath.pow(FastMath.sin(lambda), 2.0); final double e2 = FastMath.pow(b / a, 2.0); final double rt = a * Math.sqrt((cos2 + e2 * e2 * sin2) / (cos2 + e2 * sin2)); final double rt2 = rt * rt; final double deltaPsi = rangeSpacing / rt; // in radian final double r1 = Constants.halfLightSpeed * getSlantRangeTimeToFirstRangePixel(); double psi = 0.0; double alpha = 0.0; double ri = 0.0; if (!pafID.contains(UK_PAF) || processingTime.compareTo(time19930408) >= 0) { // Method 1 in Appendix B1 final double rtPlusH = Math.sqrt(rt2 + r1 * r1 + 2.0 * rt * r1 * FastMath.cos(alpha1)); final double rtPlusH2 = rtPlusH * rtPlusH; final double theta1 = FastMath.acos((r1 + rt * FastMath.cos(alpha1)) / rtPlusH); final double psi1 = alpha1 - theta1; for (int i = 0; i < sourceImageWidth; i++) { if (!isComplex) { psi = psi1 + i * deltaPsi; ri = Math.sqrt(rt2 + rtPlusH2 - 2.0 * rt * rtPlusH * FastMath.cos(psi)); } else { // see Andrea's email dataed June 9, 2010 ri = r1 + i * rangeSpacing; } alpha = FastMath.acos((rtPlusH2 - ri * ri - rt2) / (2.0 * ri * rt)); incidenceAngles[i] = alpha; lookAngles[i] = FastMath.acos((ri + rt * FastMath.cos(alpha)) / rtPlusH); rangeSpreadingLoss[i] = FastMath.pow(ri / referenceSlantRange, 3.0); } } else { // For UK-PAF products processed prior to 8th April 1993 // Method 2 in Appendix B2 final double del = getTimeIntervalBetweenDataPoints(); final double firstLineTime = getZeroDopplerAzimuthTimeOfFirstAzimuthPixel(); final double middleLineTime = getZeroDopplerAzimuthTimeOfCentreAzimuthPixel(); final int k = (int) ((middleLineTime - firstLineTime) / del + 0.5); final double[] positionVector = new double[3]; getPositionVector(k, positionVector); final double x = positionVector[0]; final double y = positionVector[1]; final double z = positionVector[2]; final double rtPlusH = Math.sqrt(x * x + y * y + z * z); final double rtPlusH2 = rtPlusH * rtPlusH; final double theta1 = FastMath.asin(FastMath.sin(alpha1) * rt / rtPlusH); final double psi1 = alpha1 - theta1; for (int i = 0; i < sourceImageWidth; i++) { if (!isComplex) { psi = psi1 + FastMath.asin(i * deltaPsi); ri = Math.sqrt(rt2 + rtPlusH2 - 2.0 * rt * rtPlusH * FastMath.cos(psi)); } else { // see Andrea's email dataed June 9, 2010 ri = r1 + i * rangeSpacing; } alpha = FastMath.acos((rtPlusH2 - ri * ri - rt2) / (2.0 * ri * rt)); incidenceAngles[i] = alpha; lookAngles[i] = FastMath.asin(FastMath.sin(alpha) * rt / rtPlusH); rangeSpreadingLoss[i] = FastMath.pow(ri / referenceSlantRange, 3.0); } } /* for (int i = 0; i < sourceImageWidth; i=i+20) { System.out.print(incidenceAngles[i] + ", "); } System.out.println(); for (int i = 0; i < sourceImageWidth; i=i+20) { System.out.print(lookAngles[i] + ", "); } System.out.println(); for (int i = 0; i < sourceImageWidth; i=i+20) { System.out.print(rangeSpreadingLoss[i] + ", "); } System.out.println(); */ } /** * Get the incidence angle at the first range pixel (in degree). * * @return The incidence angle. */ private double getIncidenceAngleAtFirstRangePixel() { // Field 56 in PRI Facility Related Data Record (in degree) final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Facility Related"); if (facility == null) { throw new OperatorException("Facility Related not found"); } final MetadataAttribute attr = facility.getAttribute("Incidence angle at first range pixel"); if (attr == null) { throw new OperatorException("Incidence angle at first range pixel not found"); } //System.out.println("Incidence angle at first range pixel is " + incidenceAngleAtFirstRangePixel); return attr.getData().getElemFloat(); } /** * Get slant range time to the first range pixel (in second). * * @return The slant range time. */ private double getSlantRangeTimeToFirstRangePixel() { // Field 126/1 in PRI Data Set Summary Record (in millisec) final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Scene Parameters"); if (facility == null) { throw new OperatorException("Scene Parameters not found"); } final MetadataAttribute attr = facility.getAttribute("Zero-doppler range time of first range pixel"); if (attr == null) { throw new OperatorException("Zero-doppler range time of first range pixel not found"); } //System.out.println("Zero-doppler range time of first range pixel is " + zeroDopplerTimeOfFirstRangePixel); return attr.getData().getElemFloat() / 1000.0; // millisec to s } /** * Get the time interval between data points (in second). * * @return The time interval. */ private double getTimeIntervalBetweenDataPoints() { // Field 20 in PRI Platform Position Data Set Record (in s) final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Platform Position"); if (facility == null) { throw new OperatorException("Platform Position not found"); } final MetadataAttribute attr = facility.getAttribute("Time interval between data points"); if (attr == null) { throw new OperatorException("Time interval between data points not found"); } return attr.getData().getElemDouble(); } /** * Get the zero doppler azimuth time of the first azimuth pixel (in second). * * @return The zero doppler time. */ private double getZeroDopplerAzimuthTimeOfFirstAzimuthPixel() { // Field 126/5 in PRI Data Set Summary Record (in UTC) final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Scene Parameters"); if (facility == null) { throw new OperatorException("Scene Parameters not found"); } final MetadataAttribute attr = facility.getAttribute("Zero-doppler azimuth time of first azimuth pixel"); if (attr == null) { throw new OperatorException("Zero-doppler azimuth time of first azimuth pixel not found"); } return convertUTCTimes(attr.getData().getElemString()); // in s } /** * Get the zero doppler azimuth time of the centre azimuth pixel (in second). * * @return The zero doppler time. */ private double getZeroDopplerAzimuthTimeOfCentreAzimuthPixel() { // Field 126/5 in PRI Data Set Summary Record (in UTC) final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Scene Parameters"); if (facility == null) { throw new OperatorException("Scene Parameters not found"); } final MetadataAttribute attr = facility.getAttribute("Zero-doppler azimuth time of centre azimuth pixel"); if (attr == null) { throw new OperatorException("Zero-doppler azimuth time of centre azimuth pixel not found"); } return convertUTCTimes(attr.getData().getElemString()); // in s } /** * Get the position vector for the kth data point. * * @param k The data point index * @param array The length-3 1D array holding position coordinates (x,y,z) */ private void getPositionVector(final int k, final double[] array) { if (k > getNumOfDataPoints()) { throw new OperatorException("Invalid data point index"); } // Fields 29 to (29 + number of data points) in PRI Platform Position Data Set Record final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Platform Position"); if (facility == null) { throw new OperatorException("Platform Position not found"); } final MetadataAttribute xAttr = facility.getAttribute("Position vector X " + k); if (xAttr == null) { throw new OperatorException("Position vector X " + k + " not found"); } final MetadataAttribute yAttr = facility.getAttribute("Position vector Y " + k); if (yAttr == null) { throw new OperatorException("Position vector X " + k + " not found"); } final MetadataAttribute zAttr = facility.getAttribute("Position vector Z " + k); if (zAttr == null) { throw new OperatorException("Position vector Y " + k + " not found"); } array[0] = xAttr.getData().getElemDouble(); array[1] = yAttr.getData().getElemDouble(); array[2] = zAttr.getData().getElemDouble(); } /** * Get the number of data points. * * @return The number of data points. */ private int getNumOfDataPoints() { // Field 14 in PRI Platform Position Data Set Record final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Platform Position"); if (facility == null) { throw new OperatorException("Platform Position not found"); } final MetadataAttribute attr = facility.getAttribute("Number of data points"); if (attr == null) { throw new OperatorException("Number of data points not found"); } return attr.getData().getElemInt(); } /** * Update the metadata in the target product. */ private void updateTargetProductMetadata() { final MetadataElement abs = AbstractMetadata.getAbstractedMetadata(targetProduct); if (applyAntennaPatternCorrection) { AbstractMetadata.setAttribute(abs, AbstractMetadata.ant_elev_corr_flag, 1); } if (applyRangeSpreadingLossCorrection) { AbstractMetadata.setAttribute(abs, AbstractMetadata.range_spread_comp_flag, 1); } if (applyReplicaPowerCorrection) { AbstractMetadata.setAttribute(abs, AbstractMetadata.replica_power_corr_flag, 1); } AbstractMetadata.setAttribute(abs, AbstractMetadata.abs_calibration_flag, 1); if (useExtXCAFile) { AbstractMetadata.setAttribute(abs, AbstractMetadata.external_calibration_file, extXCAFileName); } } /** * Compute the antenna pattern correction facotrs for a range line in current tile. * * @param x0 The x coordinate for the pixel at the upper left corner * @param w The width of the tile */ private synchronized void computeAntennaPatternCorrectionFactors(final int x0, final int w) { if (isAntPattAvailable) return; antennaPatternCorrFactor = new double[w]; if (psID.contains(VMP)) { computeAntennaPatternCorrectionFactorsForVMPProduct(x0, w); } else { // PGS (CEOS or ENVISAT) computeAntennaPatternCorrectionFactorsForPGSProduct(x0, w); } isAntPattAvailable = true; } /** * Compute antenna pattern correction factors for VMP product for a range line in current tile. * * @param x0 The x coordinate for the pixel at the upper left corner. * @param w The width of the tile. */ private void computeAntennaPatternCorrectionFactorsForVMPProduct(final int x0, final int w) { // This function implements Appendix C. if (processingTime.compareTo(time19950716) >= 0) { for (int x = x0; x < x0 + w; x++) { antennaPatternCorrFactor[x - x0] = 1.0; } return; } double theta = 0.0; if (isERS1Mission) { // ERS-1 if (pafID.contains(ESRIN) || pafID.contains(D_PAF) || pafID.contains(I_PAF)) { if (processingTime.compareTo(time19910801) >= 0 && processingTime.compareTo(time19920901) < 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternCorrFactor[x - x0] = 1.0 / g2Im(theta); } } else if (processingTime.compareTo(time19920901) >= 0 && processingTime.compareTo(time19950716) < 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternCorrFactor[x - x0] = g2Init(theta) / g2Im(theta); } } /* else if (processingTime.compareTo(time19950716) >= 0) { for (int x = x0; x < x0 + w; x++) { antennaPatternCorrFactor[x - x0] = 1.0; } } */ } else if (pafID.contains(UK_PAF)) { if (processingTime.compareTo(time19910801) >= 0 && processingTime.compareTo(time19920901) < 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternCorrFactor[x - x0] = 1.0 / g2Im(theta); } } else if (processingTime.compareTo(time19920901) >= 0 && processingTime.compareTo(time19930408) < 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternCorrFactor[x - x0] = ec(theta) * g2Init(theta) / g2Im(theta); } } else if (processingTime.compareTo(time19930408) >= 0 && processingTime.compareTo(time19950716) < 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternCorrFactor[x - x0] = g2Init(theta) / g2Im(theta); } } /* else if (processingTime.compareTo(time19950716) >= 0) { for (int x = x0; x < x0 + w; x++) { antennaPatternCorrFactor[x - x0] = 1.0; } } */ } } else { // ERS-2 //if (processingTime.compareTo(time19951016) >= 0) { for (int x = x0; x < x0 + w; x++) { antennaPatternCorrFactor[x - x0] = 1.0; } //} } } private void computeAntennaPatternCorrectionFactorsForPGSProduct(final int x0, final int w) { final double[] antennaPatternGain = new double[w]; getPGSAntennaPatternGainForCurrentTile(x0, w, antennaPatternGain); for (int x = x0; x < x0 + w; x++) { antennaPatternCorrFactor[x - x0] = 1.0 / antennaPatternGain[x - x0]; // see Andrea's email dated Nov. 11, 2008 } } private void getPGSAntennaPatternGainForCurrentTile(final int x0, final int w, final double[] array) { final double delta = 0.05; for (int x = x0; x < x0 + w; x++) { final double theta = lookAngles[x] * Constants.RTOD; // in degree final int k = (int) ((theta - elevationAngle + 5.0) / delta); final double theta1 = elevationAngle - 5.0 + k * delta; final double theta2 = theta1 + delta; final double gain1 = FastMath.pow(10.0, (double) antPatForPGS[k] / 10.0); // convert dB to linear scale final double gain2 = FastMath.pow(10.0, (double) antPatForPGS[k + 1] / 10.0); array[x - x0] = ((theta2 - theta) * gain1 + (theta - theta1) * gain2) / (theta2 - theta1); /* System.out.println("Reference elevation angle is " + elevationAngle); System.out.println("Pixel elevation angle is " + theta); System.out.println("theta1 = " + theta1); System.out.println("theta2 = " + theta2); System.out.println("gain1 = " + pattern[k]); System.out.println("gain2 = " + pattern[k+1]); System.out.println("gain = " + targetTileNewAntPat[x - x0]); */ } } /** * Compute the initial ERS-1 antenna pattern (Appendix G1) for given look angle. * * @param lookAngle The look angle (in degree) * @return The antenna pattern gain (in linear scale) */ private double g2Init(final double lookAngle) { return getAntennaPatternGain(lookAngle, appendixG1); } /** * Compute the improved ERS-1 antenna pattern (Appendix G2(a)(b)(c)) for given look angle. * * @param lookAngle The look angle (in degree) * @return The antenna pattern gain (in linear scale) */ private double g2Im(final double lookAngle) { return getAntennaPatternGain(lookAngle, appendixG2); } /** * Compute the ERS-2 antenna pattern (Appendix G3(a)(b)(c)) for given look angle. * * @param lookAngle The look angle (in degree) * @return The antenna pattern gain (in linear scale) */ private double g2ERS2(double lookAngle) { return getAntennaPatternGain(lookAngle, appendixG3); } /** * Compute initial or improved ERS-1 antenna pattern (Appendix G1 or G2) for given look angle. * * @param lookAngle The look angle (in degree) * @param array 2-D array holding antenna pattern data given in Appendix G1 or G2 * @return The antenna pattern gain (in linear scale) */ private static double getAntennaPatternGain(final double lookAngle, final double[][] array) { final int numRows = array.length; final int numCols = array[0].length; if (numCols != 2) { throw new OperatorException("Incorrect array dimension"); } final double boreSightAngle = lookAngle - relativeLookAngle; int row1 = 0; int row2 = 0; if (boreSightAngle < array[0][0]) { row1 = 0; row2 = 1; } else if (boreSightAngle > array[numRows - 1][0]) { row1 = numRows - 2; row2 = numRows - 1; } else { for (int i = 1; i < numRows; i++) { if (boreSightAngle < array[i][0]) { row1 = i - 1; row2 = i; break; } } } final double delTheta1 = array[row1][0]; final double delTheta2 = array[row2][0]; final double gain1 = array[row1][1]; final double gain2 = array[row2][1]; final double lambda = (boreSightAngle - delTheta1) / (delTheta2 - delTheta1); double gain = ((1 - lambda) * gain1 + lambda * gain2); gain = FastMath.pow(10, gain / 10.0); // dB to linear scale return gain; } /** * Compute the UK-PAF ERS-1 antenna pattern (Appendix H) for given look angle. * * @param lookAngle The look angle (in degree) * @return The antenna pattern gain (in linear scale) */ private double ec(final double lookAngle) { final int numRows = appendixH.length; final int numCols = appendixH[0].length; if (numRows < 2 || numCols < 2) { throw new OperatorException("Not enough antenna pattern data"); } final double boreSightAngle = lookAngle - relativeLookAngle; int row1 = 0; int row2 = 0; if (sceneCentreLatitude < appendixH[1][0]) { row1 = 1; row2 = 2; } else if (sceneCentreLatitude > appendixH[numRows - 1][0]) { row1 = numRows - 2; row2 = numRows - 1; } else { for (int i = 2; i < numRows; i++) { if (sceneCentreLatitude < appendixH[i][0]) { row1 = i - 1; row2 = i; break; } } } int col1 = 0; int col2 = 0; if (boreSightAngle < appendixH[0][1]) { col1 = 1; col2 = 2; } else if (boreSightAngle > appendixH[numCols - 1][0]) { col1 = numCols - 2; col2 = numCols - 1; } else { for (int j = 2; j < numCols; j++) { if (boreSightAngle < appendixH[0][j]) { col1 = j - 1; col2 = j; break; } } } final double lat1 = appendixH[row1][0]; final double lat2 = appendixH[row2][0]; final double delTheta1 = appendixH[0][col1]; final double delTheta2 = appendixH[0][col2]; if (Double.compare(lat1, lat2) == 0 || Double.compare(delTheta1, delTheta2) == 0) { throw new OperatorException("Incorrect latitude or look angle data"); } final double gain11 = appendixH[row1][col1]; final double gain12 = appendixH[row1][col2]; final double gain21 = appendixH[row2][col1]; final double gain22 = appendixH[row2][col2]; final double lambda1 = (sceneCentreLatitude - lat1) / (lat2 - lat1); final double lambda2 = (boreSightAngle - delTheta1) / (delTheta2 - delTheta1); double gain = (1 - lambda2) * ((1 - lambda1) * gain11 + lambda1 * gain21) + lambda2 * ((1 - lambda1) * gain12 + lambda1 * gain22); gain = FastMath.pow(10, gain / 10); // dB to linear scale return gain; } //================================================= ADC ====================================================== private void prepareForADCCorrection() { getADCPowerLossCorrLUT(); windowWidth = (int) (windowDimInRange / rangeSpacing); // 1200 pixels windowHeight = (int) (windowDimInAzimuth / azimuthSpacing); // 400 pixels blockWidth = (int) (downSampleBlockSize / rangeSpacing); // 8 pixels blockHeight = (int) (downSampleBlockSize / azimuthSpacing); // 8 pixels } /** * Get ADC Power Loss Correction Look-up Table (Appendix F1 or F2). */ private void getADCPowerLossCorrLUT() { if (isERS1Mission) { appendixF1 = readFile(getERSAuxFile("Appendix_F1.txt").toFile()); } else { appendixF2 = readFile(getERSAuxFile("Appendix_F2.txt").toFile()); } } private boolean isADCNeeded(final Band sourceBand1, final Band sourceBand2, final Unit.UnitType bandUnit) { final int w = Math.min(windowWidth, sourceImageWidth); final int h = Math.min(windowHeight, sourceImageHeight); final int x0 = (sourceImageWidth - w) / 2; final int y0 = (sourceImageHeight - h) / 2; final Rectangle sourceTileRectangle = new Rectangle(x0, y0, w, h); final Tile sourceRaster1 = getSourceTile(sourceBand1, sourceTileRectangle); Tile sourceRaster2 = null; if (sourceBand2 != null) { sourceRaster2 = getSourceTile(sourceBand2, sourceTileRectangle); } final ProductData srcData1 = sourceRaster1.getDataBuffer(); ProductData srcData2 = null; if (sourceRaster2 != null) srcData2 = sourceRaster2.getDataBuffer(); int index; double sigma = 0.0; for (int y = y0; y < y0 + h; y++) { for (int x = x0; x < x0 + w; x++) { index = sourceRaster1.getDataBufferIndex(x, y); if (bandUnit == Unit.UnitType.AMPLITUDE) { final double dn = srcData1.getElemDoubleAt(index); sigma += dn * dn; } else if (bandUnit == Unit.UnitType.INTENSITY) { sigma += srcData1.getElemDoubleAt(index); } else { // COMPLEX final double i = srcData1.getElemDoubleAt(index); final double q = srcData2.getElemDoubleAt(index); sigma += i * i + q * q; } } } sigma /= w * h * calibrationConstant; if (sigma < underFlowFloat) { sigma = Math.min(ers1ApplyADCThreshold, ers2ApplyADCThreshold); } else { sigma = 10.0 * Math.log10(sigma); } if (isERS1Mission && sigma <= ers1ApplyADCThreshold) { return false; } return !(!isERS1Mission && sigma <= ers2ApplyADCThreshold); } /** * Compute antenna pattern gain values for a range line. They are used to remove the antenna * pattern gain correction applied before. * * @param x0 The x coordinate for the pixel at the upper left corner * @param w The width of the tile */ private void computeAntennaPatternGain(final int x0, final int w) { antennaPatternGain = new double[w]; if (psID.contains(VMP)) { computeAntennaPatternGainForVMPProduct(x0, w); } else { // PGS (XEOS or ENVISAT) computeAntennaPatternGainForPGSProduct(x0, w); } } private void computeAntennaPatternGainForVMPProduct(final int x0, final int w) { // This function implements Appendix E. double theta = 0.0; if (isERS1Mission) { if (processingTime.compareTo(time19950716) >= 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternGain[x - x0] = g2Im(theta); } return; } if (pafID.contains(ESRIN) || pafID.contains(D_PAF) || pafID.contains(I_PAF)) { if (processingTime.compareTo(time19910801) >= 0 && processingTime.compareTo(time19920901) < 0) { for (int x = x0; x < x0 + w; x++) { antennaPatternGain[x - x0] = 1.0; } } else if (processingTime.compareTo(time19920901) >= 0 && processingTime.compareTo(time19950716) < 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternGain[x - x0] = g2Init(theta); } } } else if (pafID.contains(UK_PAF)) { if (processingTime.compareTo(time19910801) >= 0 && processingTime.compareTo(time19920901) < 0) { for (int x = x0; x < x0 + w; x++) { antennaPatternGain[x - x0] = 1.0; } } else if (processingTime.compareTo(time19920901) >= 0 && processingTime.compareTo(time19930408) < 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternGain[x - x0] = ec(theta) * g2Init(theta); } } else if (processingTime.compareTo(time19930408) >= 0 && processingTime.compareTo(time19950716) < 0) { for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternGain[x - x0] = g2Init(theta); } } } } else { // ERS-2 for (int x = x0; x < x0 + w; x++) { theta = lookAngles[x] * Constants.RTOD; // in degree antennaPatternGain[x - x0] = g2ERS2(theta); } } } private void computeAntennaPatternGainForPGSProduct(final int x0, final int w) { final double[] antennaPatternGainArray = new double[w]; getPGSAntennaPatternGainForCurrentTile(x0, w, antennaPatternGainArray); System.arraycopy(antennaPatternGainArray, x0 - x0, antennaPatternGain, x0 - x0, x0 + w - x0); } private double[][] computeADCPowerLossValuesForCurrentTile(final Band sourceBand1, final Band sourceBand2, final int tx0, final int ty0, final int tw, final int th, final Unit.UnitType bandUnit) { // 1. Get source tile rectangle TileDescriptionFlags tileDescriptionFlags = new TileDescriptionFlags(); Rectangle sourceTileRectangle = getSourceTileRectangle(tx0, ty0, tw, th, tileDescriptionFlags); // 2. Compute intensity image RenderedImage intendityImage = getIntensityImage(sourceBand1, sourceBand2, sourceTileRectangle, bandUnit); //System.out.println("intendityImage width: " + intendityImage.getWidth()); //System.out.println("intendityImage height: " + intendityImage.getHeight()); //outputRealImage(intendityImage, 0, 999); // 3. Average intensity image with an 8x8 block. RenderedImage downSampledImage = downSampleImage(intendityImage); //System.out.println("downSampledImage width: " + downSampledImage.getWidth()); //System.out.println("downSampledImage height: " + downSampledImage.getHeight()); //outputRealImage(downSampledImage, 0, 999); // 4. Removing original corrections applied: range spreading loss, antenna pattern and replica pulse power. RenderedImage rawImage = removeFactorsApplied(downSampledImage, sourceTileRectangle); //System.out.println("rawImage width: " + rawImage.getWidth()); //System.out.println("rawImage height: " + rawImage.getHeight()); //outputRealImage(rawImage, 0, 999); // 5. Smooth the raw image with a (1200/8)x(400/8) = 150x50 window. RenderedImage smoothedImage = smoothImage(rawImage); //System.out.println("smoothedImage width: " + smoothedImage.getWidth()); //System.out.println("smoothedImage height: " + smoothedImage.getHeight()); //outputRealImage(smoothedImage, 25000, 25999); // 6. Squaring the pixel value and dividing it by calibration constant K. RenderedImage squaredImage = getSquaredImage(smoothedImage); //System.out.println("squaredImage width: " + squaredImage.getWidth()); //System.out.println("squaredImage height: " + squaredImage.getHeight()); //outputRealImage(squaredImage, 25000, 25999); // 7. Generating ADC compensation file using look-up table in Appendix F (F1 or F2) and interpolation. return computeADCPowerLossValue(squaredImage, tileDescriptionFlags); } private Rectangle getSourceTileRectangle(final int tx0, final int ty0, final int tw, final int th, TileDescriptionFlags tileDescriptionFlags) { // The tile height should have window height more pixels than the target tile height final int halfWindowHeight = windowHeight / 2; final int halfWindowWidth = windowWidth / 2; int sx0 = tx0; int sy0 = ty0; int sw = tw; int sh = th; tileDescriptionFlags.adcSourceTileTopExtFlag = false; tileDescriptionFlags.adcSourceTileBottomExtFlag = false; tileDescriptionFlags.adcSourceTileLeftExtFlag = false; tileDescriptionFlags.adcSourceTileRightExtFlag = false; if (ty0 >= halfWindowHeight) { tileDescriptionFlags.adcSourceTileTopExtFlag = true; sy0 = ty0 - halfWindowHeight; sh += halfWindowHeight; } if (ty0 + th + halfWindowHeight <= sourceImageHeight) { tileDescriptionFlags.adcSourceTileBottomExtFlag = true; sh += halfWindowHeight; } if (tx0 >= halfWindowWidth) { tileDescriptionFlags.adcSourceTileLeftExtFlag = true; sx0 = tx0 - halfWindowWidth; sw += halfWindowWidth; } if (tx0 + tw + halfWindowWidth <= sourceImageWidth) { tileDescriptionFlags.adcSourceTileRightExtFlag = true; sw += halfWindowWidth; } return new Rectangle(sx0, sy0, sw, sh); } private RenderedImage getIntensityImage(final Band sourceBand1, final Band sourceBand2, final Rectangle sourceTileRectangle, final Unit.UnitType bandUnit) { final int sx0 = sourceTileRectangle.x; final int sy0 = sourceTileRectangle.y; final int sw = sourceTileRectangle.width; final int sh = sourceTileRectangle.height; final double[] array = new double[sw * sh]; final Tile sourceRaster1 = getSourceTile(sourceBand1, sourceTileRectangle); Tile sourceRaster2 = null; if (sourceBand2 != null) { sourceRaster2 = getSourceTile(sourceBand2, sourceTileRectangle); } final ProductData srcData1 = sourceRaster1.getDataBuffer(); ProductData srcData2 = null; if (sourceRaster2 != null) srcData2 = sourceRaster2.getDataBuffer(); int index; double sigma, dn, i, q; int k = 0; final int maxY = sy0 + sh; final int maxX = sx0 + sw; for (int y = sy0; y < maxY; ++y) { for (int x = sx0; x < maxX; ++x) { index = sourceRaster1.getDataBufferIndex(x, y); if (bandUnit == Unit.UnitType.AMPLITUDE) { dn = srcData1.getElemDoubleAt(index); sigma = dn * dn; } else if (bandUnit == Unit.UnitType.INTENSITY) { sigma = srcData1.getElemDoubleAt(index); } else { // COMPLEX i = srcData1.getElemDoubleAt(index); q = srcData2.getElemDoubleAt(index); sigma = i * i + q * q; } array[k++] = sigma; } } return createRenderedImage(array, sw, sh); } private static RenderedImage createRenderedImage(final double[] array, final int width, final int height) { // create rendered image with demension being width by height final SampleModel sampleModel = RasterFactory.createBandedSampleModel(DataBuffer.TYPE_DOUBLE, width, height, 1); final ColorModel colourModel = PlanarImage.createColorModel(sampleModel); final DataBufferDouble dataBuffer = new DataBufferDouble(array, array.length); final WritableRaster raster = RasterFactory.createWritableRaster(sampleModel, dataBuffer, new Point(0, 0)); return new BufferedImage(colourModel, raster, false, new Hashtable()); } private RenderedImage downSampleImage(RenderedImage intendityImage) { final double scaleX = 1.0 / blockWidth; // 1/8 = 0.125 final double scaleY = 1.0 / blockHeight; // 1/8 = 0.125 return SubsampleAverageDescriptor.create(intendityImage, scaleX, scaleY, null); } private RenderedImage removeFactorsApplied(final RenderedImage downSampledImage, final Rectangle sourceTileRectangle) { final int sx0 = sourceTileRectangle.x; final int w = downSampledImage.getWidth(); final int h = downSampledImage.getHeight(); final double[] array = new double[h * w]; //final Raster data = downSampledImage.getData(); final Raster data = downSampledImage.getData(new Rectangle(0, 0, w, h)); double sigma; int k = 0; for (int y = 0; y < h; y++) { for (int x = 0; x < w; x++) { sigma = data.getSampleDouble(x, y, 0); if (antennaPatternCorrectionFlag) { sigma *= antennaPatternGain[sx0 + x * blockWidth]; } if (rangeSpreadingLossCompFlag) { sigma /= rangeSpreadingLoss[sx0 + x * blockWidth]; } if (!isERS1Mission) { sigma /= replicaPulseVariationsCorrectionFactor; } array[k++] = Math.sqrt(sigma); } } return createRenderedImage(array, w, h); } private RenderedImage smoothImage(final RenderedImage rawImage) { final int slidingWindowWidth = windowWidth / blockWidth; final int slidingWindowHeight = windowHeight / blockHeight; final BorderExtender extender = BorderExtender.createInstance(BorderExtender.BORDER_REFLECT); final RenderingHints hints = new RenderingHints(JAI.KEY_BORDER_EXTENDER, extender); final ParameterBlock pb = new ParameterBlock(); pb.addSource(rawImage); pb.add(slidingWindowWidth); pb.add(slidingWindowHeight); pb.add(slidingWindowWidth / 2); pb.add(slidingWindowHeight / 2); return JAI.create("boxfilter", pb, hints); } private RenderedImage getSquaredImage(final RenderedImage smoothedImage) { final ParameterBlock pb1 = new ParameterBlock(); pb1.addSource(smoothedImage); pb1.addSource(smoothedImage); final RenderedImage squaredImage = JAI.create("multiply", pb1); final double[] v = { calibrationConstant }; final ParameterBlock pb2 = new ParameterBlock(); pb2.addSource(squaredImage); pb2.add(v); return JAI.create("dividebyconst", pb2, null); } private double[][] computeADCPowerLossValue(final RenderedImage squaredImage, final TileDescriptionFlags tileDescriptionFlags) { final int delH = (windowHeight / 2) / blockHeight; final int delW = (windowWidth / 2) / blockWidth; int x0 = 0; int y0 = 0; int w = squaredImage.getWidth(); int h = squaredImage.getHeight(); if (tileDescriptionFlags.adcSourceTileTopExtFlag) { y0 = delH; h -= delH; } if (tileDescriptionFlags.adcSourceTileBottomExtFlag) { h -= delH; } if (h <= 0) h = 1; if (tileDescriptionFlags.adcSourceTileLeftExtFlag) { x0 = delW; w -= delW; } if (tileDescriptionFlags.adcSourceTileRightExtFlag) { w -= delW; } if (w <= 0) w = 1; double[][] adcPowerLoss = new double[h][w]; Raster data = squaredImage.getData(); double dn; for (int y = y0; y < y0 + h; y++) { for (int x = x0; x < x0 + w; x++) { dn = data.getSampleDouble(x, y, 0); if (isERS1Mission) { adcPowerLoss[y - y0][x - x0] = getPowerLossValue(dn, appendixF1); } else { adcPowerLoss[y - y0][x - x0] = getPowerLossValue(dn, appendixF2); } } } /* for (int x = 0; x < w; x++) { System.out.print(adcPowerLoss[25][x] + ","); } System.out.println(); */ return adcPowerLoss; } // This function is for debugging only. private static void outputRealImage(final RenderedImage I, final int startIdx, final int endIdx) { final Raster data = I.getData(); final double[] real = data.getSamples(0, 0, I.getWidth(), I.getHeight(), 0, (double[]) null); for (int i = startIdx; i <= endIdx; i++) { System.out.print(real[i] + ","); } System.out.println(); } /** * Compute ADC power loss value for given pixel value using LUT in Appendix F1 or F2. * * @param dn The pixel value * @param array 2-D array holding ADC power loss value data given in Appendix F1 or F2 * @return The ADC power loss value (in linear scale) */ private static double getPowerLossValue(final double dn, final double[][] array) { final int numRows = array.length; final int numCols = array[0].length; if (numCols != 2) { throw new OperatorException("Incorrect array dimension"); } double dnInDb; if (dn < underFlowFloat) { return -underFlowFloat; } else { dnInDb = 10.0 * Math.log10(dn); } int row1 = 0; int row2 = 0; if (dnInDb < array[0][0]) { row1 = 0; row2 = 1; } else if (dnInDb > array[numRows - 1][0]) { row1 = numRows - 2; row2 = numRows - 1; } else { for (int i = 1; i < numRows; i++) { if (dnInDb < array[i][0]) { row1 = i - 1; row2 = i; break; } } } final double intensityK1 = array[row1][0]; final double intensityK2 = array[row2][0]; final double loss1 = array[row1][1]; final double loss2 = array[row2][1]; final double lambda = (dnInDb - intensityK1) / (intensityK2 - intensityK1); double loss = (1 - lambda) * loss1 + lambda * loss2; loss = FastMath.pow(10.0, loss / 10.0); // dB to linear scale return loss; } // ================================== PGS ENVISAT ============================================================ /** * Get the processing facility identifier. */ private void getProcessingFacilityIDFromENVISAT() { MetadataElement mph = origMetadataRoot.getElement("MPH"); if (mph == null) { throw new OperatorException("MPH not found"); } MetadataAttribute attr = mph.getAttribute("proc_center"); if (attr == null) { throw new OperatorException("proc_center not found"); } pafID = attr.getData().getElemString(); /* if (!pafID.contains(ESRIN) && !pafID.contains(D_PAF) && !pafID.contains(I_PAF) && !pafID.contains(UK_PAF)) { throw new OperatorException("Invalid processing facility identifier: " + pafID); } */ //System.out.println("Processing facility identifier is " + pafID); } /** * Get number of records in Main Processing Params data set. */ private void getNumOfRecordsInMainProcParam() { MetadataElement dsd = origMetadataRoot.getElement("DSD").getElement("DSD.3"); if (dsd == null) { throw new OperatorException("DSD not found"); } MetadataAttribute numRecordsAttr = dsd.getAttribute("num_records"); if (numRecordsAttr == null) { throw new OperatorException("num_records not found"); } 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); } /** * Get calibration factors from Metadata for each band in the product. */ private void getCalibrationConstantFromENVISAT() { MetadataElement ads; if (numMPPRecords == 1) { ads = origMetadataRoot.getElement("MAIN_PROCESSING_PARAMS_ADS"); } else { ads = origMetadataRoot.getElement("MAIN_PROCESSING_PARAMS_ADS") .getElement("MAIN_PROCESSING_PARAMS_ADS.1"); } if (ads == null) { throw new OperatorException("MAIN_PROCESSING_PARAMS_ADS not found"); } MetadataAttribute calibrationFactorsAttr = ads .getAttribute("ASAR_Main_ADSR.sd/calibration_factors.1.ext_cal_fact"); if (calibrationFactorsAttr == null) { throw new OperatorException("calibration_factors.1.ext_cal_fact not found"); } calibrationConstant = (double) calibrationFactorsAttr.getData().getElemFloat(); /* calibrationFactor[0] = (double) calibrationFactorsAttr.getData().getElemFloat(); calibrationFactorsAttr = ads.getAttribute("ASAR_Main_ADSR.sd/calibration_factors.2.ext_cal_fact"); if (calibrationFactorsAttr == null) { throw new OperatorException("calibration_factors.2.ext_cal_fact not found"); } calibrationFactor[1] = (double) calibrationFactorsAttr.getData().getElemFloat(); if (Double.compare(calibrationFactor[0], 0.0) == 0 && Double.compare(calibrationFactor[1], 0.0) == 0) { throw new OperatorException("Calibration factors in metadata are zero"); } */ //System.out.println("calibration factor for band 1 is " + calibrationFactor[0]); //System.out.println("calibration factor for band 2 is " + calibrationFactor[1]); } /** * Compute incidence angles (in radian), look angles (in radian) and range spreading loss * for pixels in a complete range line. */ private void computeIncidenceAnglesLookAnglesRangeSpreadingLossForENVISAT() { // private void computeIncidenceAnglesLookAnglesRangeSpreadingLoss() { incidenceAngles = new double[sourceImageWidth]; lookAngles = new double[sourceImageWidth]; rangeSpreadingLoss = new double[sourceImageWidth]; final TiePointGrid incidenceAngleTiePointGrid = OperatorUtils.getIncidenceAngle(sourceProduct); final TiePointGrid slantRangeTimeTiePointGrid = OperatorUtils.getSlantRangeTime(sourceProduct); final double rSat = getSatelliteToEarthCenterDistanceForENVISAT(); /* double rSat; if (isCEOSFormat) { rSat = getSatelliteToEarthCenterDistanceForCEOS(); } else { rSat = getSatelliteToEarthCenterDistanceForENVISAT(); } */ final int y = sourceImageHeight / 2; for (int x = 0; x < sourceImageWidth; x++) { final double alpha = incidenceAngleTiePointGrid.getPixelDouble(x + 0.5, y + 0.5) * Constants.DTOR; // in radian final double time = slantRangeTimeTiePointGrid.getPixelDouble(x + 0.5, y + 0.5) / Constants.oneBillion; //convert ns to s final double r = time * Constants.halfLightSpeed; // in m final double theta = alpha - FastMath.asin(FastMath.sin(alpha) * r / rSat); // in radian incidenceAngles[x] = alpha; lookAngles[x] = theta; rangeSpreadingLoss[x] = FastMath.pow(r / referenceSlantRange, 3.0); } /* for (int i = 0; i < sourceImageWidth; i=i+20) { System.out.print(incidenceAngles[i] + ", "); } System.out.println(); for (int i = 0; i < sourceImageWidth; i=i+20) { System.out.print(lookAngles[i] + ", "); } System.out.println(); for (int i = 0; i < sourceImageWidth; i=i+20) { System.out.print(rangeSpreadingLoss[i] + ", "); } System.out.println(); */ } /** * Compute distance from satellite to the Earth center using satellite corrodinate in Metadata. * * @return The distance. */ private double getSatelliteToEarthCenterDistanceForENVISAT() { final MetadataElement mppAds = origMetadataRoot.getElement("MAIN_PROCESSING_PARAMS_ADS"); if (mppAds == null) { throw new OperatorException("MAIN_PROCESSING_PARAMS_ADS not found"); } MetadataElement ads; if (numMPPRecords == 1) { ads = mppAds; } else { ads = mppAds.getElement("MAIN_PROCESSING_PARAMS_ADS." + 1); } final MetadataAttribute xPositionAttr = ads.getAttribute("ASAR_Main_ADSR.sd/orbit_state_vectors.3.x_pos_1"); if (xPositionAttr == null) { throw new OperatorException("x_pos_1 not found"); } final float x_pos = xPositionAttr.getData().getElemInt() / 100.0f; // divide 100 to convert unit from 10^-2 m to m //System.out.println("x position is " + x_pos); final MetadataAttribute yPositionAttr = ads.getAttribute("ASAR_Main_ADSR.sd/orbit_state_vectors.3.y_pos_1"); if (yPositionAttr == null) { throw new OperatorException("y_pos_1 not found"); } final float y_pos = yPositionAttr.getData().getElemInt() / 100.0f; // divide 100 to convert unit from 10^-2 m to m //System.out.println("y position is " + y_pos); final MetadataAttribute zPositionAttr = ads.getAttribute("ASAR_Main_ADSR.sd/orbit_state_vectors.3.z_pos_1"); if (zPositionAttr == null) { throw new OperatorException("z_pos_1 not found"); } final float z_pos = zPositionAttr.getData().getElemInt() / 100.0f; // divide 100 to convert unit from 10^-2 m to m //System.out.println("z position is " + z_pos); final double rSat = Math.sqrt(x_pos * x_pos + y_pos * y_pos + z_pos * z_pos); // in m if (Double.compare(rSat, 0.0) == 0) { throw new OperatorException("x, y and z positions in orbit_state_vectors are all zeros"); } return rSat; } /** * Get the replica pulse power for PGS ENVISAT product. * * @return The replica pulse power. */ private double getReplicaPulsePowerForENVISAT() { // Field 9 in CHIRP_PARAMS_ADS final MetadataElement dsd = origMetadataRoot.getElement("CHIRP_PARAMS_ADS"); if (dsd == null) { throw new OperatorException("CHIRP_PARAMS_ADS not found"); } final MetadataAttribute attr = dsd.getAttribute("chirp_power"); if (attr == null) { throw new OperatorException("chirp_power not found"); } double replicaPulsePower = attr.getData().getElemFloat(); // in dB replicaPulsePower = FastMath.pow(10.0, replicaPulsePower / 10.0); // convert to linear scale //System.out.println("Replica pulse power is " + replicaPulsePower); return replicaPulsePower; } private double getSatelliteToEarthCenterDistanceForCEOS() { // Field 99 - 101 in PRI Facility Related Data Record final MetadataElement facility = origMetadataRoot.getElement("Leader").getElement("Facility Related"); if (facility == null) { throw new OperatorException("Facility Related not found"); } final MetadataAttribute xAttr = facility.getAttribute("Input state vector - Position vector X"); if (xAttr == null) { throw new OperatorException("Input state vector - Position vector X"); } final double x = xAttr.getData().getElemDouble(); final MetadataAttribute yAttr = facility.getAttribute("Input state vector - Position vector Y"); if (yAttr == null) { throw new OperatorException("Input state vector - Position vector Y"); } final double y = yAttr.getData().getElemDouble(); final MetadataAttribute zAttr = facility.getAttribute("Input state vector - Position vector Z"); if (zAttr == null) { throw new OperatorException("Input state vector - Position vector Z"); } final double z = zAttr.getData().getElemDouble(); return Math.sqrt(x * x + y * y + z * z); } /** * Gets a {@link Tile} for a given band and rectangle. * * @param rasterDataNode the raster data node of a data product, * e.g. a {@link org.esa.snap.framework.datamodel.Band Band} or * {@link org.esa.snap.framework.datamodel.TiePointGrid TiePointGrid}. * @param rectangle the raster rectangle in pixel coordinates * @return a tile. * @throws OperatorException if the tile request cannot be processed */ public Tile getSourceTile(RasterDataNode rasterDataNode, Rectangle rectangle) throws OperatorException { return calibrationOp.getSourceTile(rasterDataNode, rectangle); } //==================================== pixel calibration used by RD ====================================== public double applyRetroCalibration(int x, int y, double v, String bandPolar, final Unit.UnitType bandUnit, int[] subSwathIndex) { return v; } 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 bandPolar, final Unit.UnitType bandUnit, final int[] subSwathIndex) { // For both detectec and slant range products, // 1) local incidence angle (Remember that for ERS the correction is sin(theta_loc)/sin(theta_ref)) // 2) antenna pattern // 3) range spreading loss // 4) replica pulse power (for both ERS-1 and ERS-2 since for ERS2 it has been removed in the pre-calibration step) // 5) calibration constant 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"); } if (multilookFlag && antennaPatternCorrectionFlag) { // calibration constant and incidence angle corrections only return FastMath.sin(Math.abs(localIncidenceAngle) * Constants.DTOR) / FastMath.sin(referenceIncidenceAngle) / calibrationConstant; } sigma *= FastMath.sin(Math.abs(localIncidenceAngle) * Constants.DTOR) / FastMath.sin(referenceIncidenceAngle); sigma /= getNewAntennaPatternGainSquare((int) rangeIndex); sigma *= rangeSpreadingLoss[(int) rangeIndex]; sigma *= replicaPulseVariationsCorrectionFactor; sigma /= calibrationConstant; return sigma; } /** * Get the new antenna pattern gain square for a given pixel. * * @param rangeIndex The x coordinate for the pixel * @return The antenna pattern gain square. */ private double getNewAntennaPatternGainSquare(final int rangeIndex) { if (psID.contains(VMP)) { return getNewAntennaPatternGainSquareForVMPProduct(rangeIndex); } else { // PGS (CEOS or ENVISAT) return getNewAntennaPatternGainSquareForPGSProduct(rangeIndex); } } /** * Get the new antenna pattern gain for a given pixel for VMP product. * * @param rangeIndex The x coordinate for the pixel * @return The antenna pattern gain square. */ private double getNewAntennaPatternGainSquareForVMPProduct(final int rangeIndex) { if (isERS1Mission) { // ERS-1 return g2Im(lookAngles[rangeIndex] * Constants.RTOD); } else { return g2ERS2(lookAngles[rangeIndex] * Constants.RTOD); } } private double getNewAntennaPatternGainSquareForPGSProduct(final int rangeIndex) { final double delta = 0.05; final double theta = lookAngles[rangeIndex] * Constants.RTOD; // in degree final int k = (int) ((theta - elevationAngle + 5.0) / delta); final double theta1 = elevationAngle - 5.0 + k * delta; final double theta2 = theta1 + delta; final double gain1 = FastMath.pow(10.0, (double) antPatForPGS[k] / 10.0); // convert dB to linear scale final double gain2 = FastMath.pow(10.0, (double) antPatForPGS[k + 1] / 10.0); final double gain = ((theta2 - theta) * gain1 + (theta - theta1) * gain2) / (theta2 - theta1); return gain; // see Andrea's email dated Nov. 11, 2008 } public void removeFactorsForCurrentTile(Band targetBand, Tile targetTile, String srcBandName) throws OperatorException { // For ground range product, // a) remove antenna pattern gain // b) remove range spreading loss corrections // c) remove replica pulse variation (ERS-2 only) // d) multiply calibration constant (done in applyCalibration) // e) apply ADC power loss correction // // For slant range complex product, // c) remove replica pulse variation (ERS-2 only) // d) multiply calibration constant (done in applyCalibration) // e) apply ADC power loss correction final Rectangle targetTileRectangle = targetTile.getRectangle(); final int tx0 = targetTileRectangle.x; final int ty0 = targetTileRectangle.y; final int tw = targetTileRectangle.width; final int th = targetTileRectangle.height; final ProductData trgData = targetTile.getDataBuffer(); //System.out.println("RetroOp: tx0 = " + tx0 + ", ty0 = " + ty0 + ", tw = " + tw + ", th = " + th); final Band sourceBand1 = sourceProduct.getBand(srcBandName); final Tile sourceTile = getSourceTile(sourceBand1, targetTileRectangle); final ProductData srcData = sourceTile.getDataBuffer(); final String[] srcBandNames = { targetBand.getName() }; Band sourceBand2 = null; if (srcBandNames.length > 1) { sourceBand2 = sourceProduct.getBand(srcBandNames[1]); } final Unit.UnitType bandUnit = Unit.getUnitType(sourceBand1); if (applyADCSaturationCorrection && !adcHasBeenTestedFlag) { testADC(sourceBand1, sourceBand2, bandUnit); } boolean applyADCSaturationCorrectionToCurrentTile = false; if (applyADCSaturationCorrection && th >= blockHeight && tw >= blockWidth) { applyADCSaturationCorrectionToCurrentTile = true; } double[][] adcPowerLoss = null; if (applyADCSaturationCorrectionToCurrentTile) { adcPowerLoss = computeADCPowerLossValuesForCurrentTile(sourceBand1, sourceBand2, tx0, ty0, tw, th, bandUnit); } double sigma = 0.0; int adcJ = 0; for (int x = tx0; x < tx0 + tw; x++) { double antennaPatternByRangeSpreadingLoss = 0.0; if (!isComplex) { antennaPatternByRangeSpreadingLoss = antennaPatternGain[x] / rangeSpreadingLoss[x]; } if (applyADCSaturationCorrectionToCurrentTile) { adcJ = Math.min(((x - tx0) / blockWidth), adcPowerLoss[0].length - 1); } for (int y = ty0; y < ty0 + th; y++) { final int srcIndex = sourceTile.getDataBufferIndex(x, y); if (bandUnit == Unit.UnitType.AMPLITUDE) { final double dn = srcData.getElemDoubleAt(srcIndex); sigma = dn * dn; } else if (bandUnit == Unit.UnitType.AMPLITUDE_DB) { sigma = FastMath.pow(10, srcData.getElemDoubleAt(srcIndex) / 5.0); } else if (bandUnit == Unit.UnitType.INTENSITY) { sigma = srcData.getElemDoubleAt(srcIndex); } else if (bandUnit == Unit.UnitType.INTENSITY_DB) { sigma = FastMath.pow(10, srcData.getElemDoubleAt(srcIndex) / 10.0); } else { throw new OperatorException("ERSCalibrator: Unknown band unit"); } if (!isComplex) { // ground range sigma *= antennaPatternByRangeSpreadingLoss; } if (!isERS1Mission) { sigma /= replicaPulseVariationsCorrectionFactor; } if (applyADCSaturationCorrectionToCurrentTile) { final int adcI = Math.min(((y - ty0) / blockHeight), adcPowerLoss.length - 1); sigma *= adcPowerLoss[adcI][adcJ]; } if (bandUnit == Unit.UnitType.AMPLITUDE) { trgData.setElemDoubleAt(srcIndex, Math.sqrt(sigma)); } else if (bandUnit == Unit.UnitType.AMPLITUDE_DB) { trgData.setElemDoubleAt(srcIndex, 5.0 * Math.log10(sigma)); } else if (bandUnit == Unit.UnitType.INTENSITY) { trgData.setElemDoubleAt(srcIndex, sigma); } else if (bandUnit == Unit.UnitType.INTENSITY_DB) { trgData.setElemDoubleAt(srcIndex, 10.0 * Math.log10(sigma)); } } } } public static class TileDescriptionFlags { public boolean adcSourceTileTopExtFlag; public boolean adcSourceTileBottomExtFlag; public boolean adcSourceTileLeftExtFlag; public boolean adcSourceTileRightExtFlag; public TileDescriptionFlags() { } } }