Examples with LocalGeometry org.esa.snap.engine_utilities.eo.LocalGeometry used on opensource projects

Example 1 with LocalGeometry

use of org.esa.snap.engine_utilities.eo.LocalGeometry in project s1tbx by senbox-org.

the class RangeDopplerGeocodingOp method computeTileStack.

/**
 * Called by the framework in order to compute the stack of tiles for the given target bands.
 * <p>The default implementation throws a runtime exception with the message "not implemented".</p>
 *
 * @param targetTiles     The current tiles to be computed for each target band.
 * @param targetRectangle The area in pixel coordinates to be computed (same for all rasters in <code>targetRasters</code>).
 * @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 rasters.
 */
@Override
public void computeTileStack(Map<Band, Tile> targetTiles, Rectangle targetRectangle, ProgressMonitor pm) throws OperatorException {
    try {
        processingStarted = true;
        try {
            if (!isElevationModelAvailable) {
                getElevationModel();
            }
        } catch (Exception e) {
            throw new OperatorException(e);
        }
        if (saveLayoverShadowMask && !isLayoverShadowMaskAvailable) {
            createLayoverShadowMask();
        }
        final int x0 = targetRectangle.x;
        final int y0 = targetRectangle.y;
        final int w = targetRectangle.width;
        final int h = targetRectangle.height;
        // System.out.println("x0 = " + x0 + ", y0 = " + y0 + ", w = " + w + ", h = " + h);
        final TileGeoreferencing tileGeoRef = new TileGeoreferencing(targetProduct, x0 - 1, y0 - 1, w + 2, h + 2);
        double[][] localDEM = new double[h + 2][w + 2];
        if (useAvgSceneHeight) {
            DEMFactory.fillDEM(localDEM, avgSceneHeight);
        } else {
            final boolean valid = DEMFactory.getLocalDEM(dem, demNoDataValue, demResamplingMethod, tileGeoRef, x0, y0, w, h, sourceProduct, nodataValueAtSea, localDEM);
            if (!valid && nodataValueAtSea) {
                for (Band targetBand : targetTiles.keySet()) {
                    ProductData data = targetTiles.get(targetBand).getRawSamples();
                    double nodatavalue = targetBand.getNoDataValue();
                    final int length = data.getNumElems();
                    for (int i = 0; i < length; ++i) {
                        data.setElemDoubleAt(i, nodatavalue);
                    }
                }
                return;
            }
        }
        final Rectangle sourceRectangle = getSourceRectangle(x0, y0, w, h, tileGeoRef, localDEM);
        final GeoPos geoPos = new GeoPos();
        final PositionData posData = new PositionData();
        final int srcMaxRange = sourceImageWidth - 1;
        final int srcMaxAzimuth = sourceImageHeight - 1;
        ProductData demBuffer = null, latBuffer = null, lonBuffer = null, localIncidenceAngleBuffer = null, projectedLocalIncidenceAngleBuffer = null, incidenceAngleFromEllipsoidBuffer = null, layoverShadowMaskBuffer = null;
        final List<TileData> tgtTileList = new ArrayList<>();
        final Set<Band> keySet = targetTiles.keySet();
        for (Band targetBand : keySet) {
            if (targetBand.equals(elevationBand)) {
                demBuffer = targetTiles.get(targetBand).getDataBuffer();
                continue;
            }
            if (targetBand.getName().equals("latitude")) {
                latBuffer = targetTiles.get(targetBand).getDataBuffer();
                continue;
            }
            if (targetBand.getName().equals("longitude")) {
                lonBuffer = targetTiles.get(targetBand).getDataBuffer();
                continue;
            }
            if (targetBand.getName().equals("localIncidenceAngle")) {
                localIncidenceAngleBuffer = targetTiles.get(targetBand).getDataBuffer();
                continue;
            }
            if (targetBand.getName().equals("projectedLocalIncidenceAngle")) {
                projectedLocalIncidenceAngleBuffer = targetTiles.get(targetBand).getDataBuffer();
                continue;
            }
            if (targetBand.getName().equals("incidenceAngleFromEllipsoid")) {
                incidenceAngleFromEllipsoidBuffer = targetTiles.get(targetBand).getDataBuffer();
                continue;
            }
            if (targetBand.getName().equals("layoverShadowMask")) {
                layoverShadowMaskBuffer = targetTiles.get(targetBand).getDataBuffer();
                continue;
            }
            final Band[] srcBands = targetBandNameToSourceBand.get(targetBand.getName());
            Tile sourceTileI = null, sourceTileQ = null;
            if (sourceRectangle != null) {
                sourceTileI = getSourceTile(srcBands[0], sourceRectangle);
                sourceTileQ = srcBands.length > 1 ? getSourceTile(srcBands[1], sourceRectangle) : null;
            }
            final TileData td = new TileData(targetTiles.get(targetBand), srcBands, isPolsar, outputComplex, targetBand.getName(), getBandUnit(targetBand.getName()), absRoot, calibrator, imgResampling, sourceTileI, sourceTileQ);
            td.applyRadiometricNormalization = targetBandApplyRadiometricNormalizationFlag.get(targetBand.getName());
            td.applyRetroCalibration = targetBandApplyRetroCalibrationFlag.get(targetBand.getName());
            tgtTileList.add(td);
        }
        final TileData[] tgtTiles = tgtTileList.toArray(new TileData[0]);
        for (TileData tileData : tgtTiles) {
            if (sourceRectangle != null) {
                try {
                    final Band[] srcBands = targetBandNameToSourceBand.get(tileData.bandName);
                    tileData.imgResamplingRaster.setSourceTiles(getSourceTile(srcBands[0], sourceRectangle), srcBands.length > 1 ? getSourceTile(srcBands[1], sourceRectangle) : null);
                } catch (Exception e) {
                    tileData.imgResamplingRaster.setSourceTiles(null, null);
                }
            } else {
                tileData.imgResamplingRaster.setSourceTiles(null, null);
            }
        }
        final int maxY = y0 + h;
        final int maxX = x0 + w;
        final EarthGravitationalModel96 egm = EarthGravitationalModel96.instance();
        final GeoPos posFirst = targetProduct.getSceneGeoCoding().getGeoPos(new PixelPos(0, 0), null);
        final GeoPos posLast = targetProduct.getSceneGeoCoding().getGeoPos(new PixelPos(0, targetImageHeight), null);
        int diffLat = (int) Math.abs(posFirst.lat - posLast.lat);
        for (int y = y0; y < maxY; y++) {
            final int yy = y - y0 + 1;
            for (int x = x0; x < maxX; x++) {
                final int index = tgtTiles[0].targetTile.getDataBufferIndex(x, y);
                Double alt = localDEM[yy][x - x0 + 1];
                if (alt.equals(demNoDataValue) && !useAvgSceneHeight) {
                    if (nodataValueAtSea) {
                        saveNoDataValueToTarget(index, tgtTiles, demBuffer);
                        continue;
                    }
                }
                tileGeoRef.getGeoPos(x, y, geoPos);
                final double lat = geoPos.lat;
                double lon = geoPos.lon;
                if (lon >= 180.0) {
                    lon -= 360.0;
                }
                if (alt.equals(demNoDataValue) && !nodataValueAtSea) {
                    // get corrected elevation for 0
                    alt = (double) egm.getEGM(lat, lon);
                }
                if (!getPosition(lat, lon, alt, posData)) {
                    saveNoDataValueToTarget(index, tgtTiles, demBuffer);
                    continue;
                }
                if (!SARGeocoding.isValidCell(posData.rangeIndex, posData.azimuthIndex, lat, lon, diffLat, sourceProduct.getSceneGeoCoding(), srcMaxRange, srcMaxAzimuth, posData.sensorPos)) {
                    saveNoDataValueToTarget(index, tgtTiles, demBuffer);
                } else {
                    final double[] localIncidenceAngles = { SARGeocoding.NonValidIncidenceAngle, SARGeocoding.NonValidIncidenceAngle };
                    if (saveLocalIncidenceAngle || saveProjectedLocalIncidenceAngle || saveSigmaNought) {
                        final LocalGeometry localGeometry = new LocalGeometry(x, y, tileGeoRef, posData.earthPoint, posData.sensorPos);
                        SARGeocoding.computeLocalIncidenceAngle(localGeometry, demNoDataValue, saveLocalIncidenceAngle, saveProjectedLocalIncidenceAngle, saveSigmaNought, x0, y0, x, y, localDEM, // in degrees
                        localIncidenceAngles);
                        if (saveLocalIncidenceAngle && localIncidenceAngles[0] != SARGeocoding.NonValidIncidenceAngle) {
                            localIncidenceAngleBuffer.setElemDoubleAt(index, localIncidenceAngles[0]);
                        }
                        if (saveProjectedLocalIncidenceAngle && localIncidenceAngles[1] != SARGeocoding.NonValidIncidenceAngle) {
                            projectedLocalIncidenceAngleBuffer.setElemDoubleAt(index, localIncidenceAngles[1]);
                        }
                    }
                    if (saveDEM) {
                        demBuffer.setElemDoubleAt(index, alt);
                    }
                    if (saveLatLon) {
                        latBuffer.setElemDoubleAt(index, lat);
                        lonBuffer.setElemDoubleAt(index, lon);
                    }
                    if (saveIncidenceAngleFromEllipsoid && incidenceAngle != null) {
                        incidenceAngleFromEllipsoidBuffer.setElemDoubleAt(index, incidenceAngle.getPixelDouble(posData.rangeIndex, posData.azimuthIndex));
                    }
                    if (saveLayoverShadowMask) {
                        layoverShadowMaskBuffer.setElemIntAt(index, layoverShadowMask[(int) (posData.azimuthIndex + 0.5)][(int) (posData.rangeIndex + 0.5)]);
                    }
                    double satelliteHeight = 0;
                    double sceneToEarthCentre = 0;
                    if (saveSigmaNought) {
                        satelliteHeight = Math.sqrt(posData.sensorPos.x * posData.sensorPos.x + posData.sensorPos.y * posData.sensorPos.y + posData.sensorPos.z * posData.sensorPos.z);
                        sceneToEarthCentre = Math.sqrt(posData.earthPoint.x * posData.earthPoint.x + posData.earthPoint.y * posData.earthPoint.y + posData.earthPoint.z * posData.earthPoint.z);
                    }
                    for (TileData tileData : tgtTiles) {
                        int[] subSwathIndex = { INVALID_SUB_SWATH_INDEX };
                        double v = getPixelValue(posData.azimuthIndex, posData.rangeIndex, tileData, subSwathIndex);
                        if (v != tileData.noDataValue && tileData.applyRadiometricNormalization) {
                            if (localIncidenceAngles[1] != SARGeocoding.NonValidIncidenceAngle) {
                                v = calibrator.applyCalibration(v, posData.rangeIndex, posData.azimuthIndex, posData.slantRange, satelliteHeight, sceneToEarthCentre, localIncidenceAngles[1], tileData.bandName, tileData.bandPolar, tileData.bandUnit, subSwathIndex);
                            // use projected incidence angle
                            } else {
                                // v = tileData.noDataValue;
                                saveNoDataValueToTarget(index, tgtTiles, demBuffer);
                                continue;
                            }
                        }
                        tileData.tileDataBuffer.setElemDoubleAt(index, v);
                    }
                    orthoDataProduced = true;
                }
            }
        }
        localDEM = null;
    } catch (Throwable e) {
        // to prevent multiple error messages
        orthoDataProduced = true;
        OperatorUtils.catchOperatorException(getId(), e);
    }
}

Example 2 with LocalGeometry

use of org.esa.snap.engine_utilities.eo.LocalGeometry in project s1tbx by senbox-org.

the class SARSimulationOp method computeTileStack.

/**
 * Called by the framework in order to compute the stack of tiles for the given target bands.
 * <p>The default implementation throws a runtime exception with the message "not implemented".</p>
 *
 * @param targetTiles     The current tiles to be computed for each target band.
 * @param targetRectangle The area in pixel coordinates to be computed (same for all rasters in <code>targetRasters</code>).
 * @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 rasters.
 */
@Override
public void computeTileStack(Map<Band, Tile> targetTiles, Rectangle targetRectangle, ProgressMonitor pm) throws OperatorException {
    final int x0 = targetRectangle.x;
    final int y0 = targetRectangle.y;
    final int w = targetRectangle.width;
    final int h = targetRectangle.height;
    // System.out.println("x0 = " + x0 + ", y0 = " + y0 + ", w = " + w + ", h = " + h);
    OverlapPercentage tileOverlapPercentage = null;
    try {
        if (!isElevationModelAvailable) {
            getElevationModel();
        }
        tileOverlapPercentage = computeTileOverlapPercentage(x0, y0, w, h);
    } catch (Exception e) {
        throw new OperatorException(e);
    }
    final Tile targetTile = targetTiles.get(targetProduct.getBand(SIMULATED_BAND_NAME));
    final ProductData masterBuffer = targetTile.getDataBuffer();
    ProductData demBandBuffer = null;
    ProductData zeroHeightBandBuffer = null;
    ProductData localIncidenceAngleBandBuffer = null;
    ProductData layoverShadowMaskBuffer = null;
    if (saveDEM) {
        demBandBuffer = targetTiles.get(targetProduct.getBand(demBandName)).getDataBuffer();
    }
    if (saveZeroHeightSimulation) {
        zeroHeightBandBuffer = targetTiles.get(targetProduct.getBand(zeroHeightSimulationBandName)).getDataBuffer();
    }
    if (saveLocalIncidenceAngle) {
        localIncidenceAngleBandBuffer = targetTiles.get(targetProduct.getBand(simulatedLocalIncidenceAngleBandName)).getDataBuffer();
    }
    if (saveLayoverShadowMask) {
        layoverShadowMaskBuffer = targetTiles.get(targetProduct.getBand(layoverShadowMaskBandName)).getDataBuffer();
    }
    final int ymin = Math.max(y0 - (int) (h * tileOverlapPercentage.tileOverlapUp), 0);
    final int ymax = Math.min(y0 + h + (int) (h * tileOverlapPercentage.tileOverlapDown), sourceImageHeight);
    final int xmin = Math.max(x0 - (int) (w * tileOverlapPercentage.tileOverlapLeft), 0);
    final int xmax = Math.min(x0 + w + (int) (w * tileOverlapPercentage.tileOverlapRight), sourceImageWidth);
    final SARPosition sarPosition = new SARPosition(firstLineUTC, lastLineUTC, lineTimeInterval, wavelength, rangeSpacing, sourceImageWidth, srgrFlag, nearEdgeSlantRange, nearRangeOnLeft, orbit, srgrConvParams);
    sarPosition.setTileConstraints(x0, y0, w, h);
    final SARPosition.PositionData posData = new SARPosition.PositionData();
    final GeoPos geoPos = new GeoPos();
    double[] slrs = null;
    double[] elev = null;
    double[] azIndex = null;
    double[] rgIndex = null;
    boolean[] savePixel = null;
    try {
        if (reGridMethod) {
            final double[] latLonMinMax = new double[4];
            computeImageGeoBoundary(xmin, xmax, ymin, ymax, latLonMinMax);
            final double latMin = latLonMinMax[0];
            final double latMax = latLonMinMax[1];
            final double lonMin = latLonMinMax[2];
            final double lonMax = latLonMinMax[3];
            final int nLat = (int) ((latMax - latMin) / delLat) + 1;
            final int nLon = (int) ((lonMax - lonMin) / delLon) + 1;
            final double[][] tileDEM = new double[nLat + 1][nLon + 1];
            final double[][] neighbourDEM = new double[3][3];
            Double alt;
            if (saveLayoverShadowMask) {
                slrs = new double[nLon];
                elev = new double[nLon];
                azIndex = new double[nLon];
                rgIndex = new double[nLon];
                savePixel = new boolean[nLon];
            }
            for (int i = 0; i < nLat; i++) {
                final double lat = latMin + i * delLat;
                if (saveLayoverShadowMask) {
                    Arrays.fill(slrs, 0.0);
                    Arrays.fill(elev, 0.0);
                    Arrays.fill(azIndex, 0.0);
                    Arrays.fill(rgIndex, 0.0);
                    Arrays.fill(savePixel, Boolean.FALSE);
                }
                for (int j = 0; j < nLon; j++) {
                    double lon = lonMin + j * delLon;
                    if (lon >= 180.0) {
                        lon -= 360.0;
                    }
                    if (saveZeroHeightSimulation) {
                        alt = 1.0;
                    } else {
                        geoPos.setLocation(lat, lon);
                        alt = dem.getElevation(geoPos);
                        if (alt.equals(demNoDataValue))
                            continue;
                    }
                    tileDEM[i][j] = alt;
                    GeoUtils.geo2xyzWGS84(lat, lon, alt, posData.earthPoint);
                    if (!sarPosition.getPosition(posData))
                        continue;
                    final LocalGeometry localGeometry = new LocalGeometry(lat, lon, delLat, delLon, posData.earthPoint, posData.sensorPos);
                    final double[] localIncidenceAngles = { SARGeocoding.NonValidIncidenceAngle, SARGeocoding.NonValidIncidenceAngle };
                    int r = 0;
                    for (int ii = Math.max(0, i - 1); ii <= i + 1; ++ii) {
                        ii = Math.min(nLat, ii);
                        int c = 0;
                        double neighbourLat = latMin + ii * delLat;
                        for (int jj = Math.max(0, j - 1); jj <= j + 1; ++jj) {
                            jj = Math.min(nLon, jj);
                            neighbourDEM[r][c] = tileDEM[ii][jj];
                            if (neighbourDEM[r][c] == 0) {
                                if (saveZeroHeightSimulation) {
                                    neighbourDEM[r][c] = 1;
                                } else {
                                    geoPos.setLocation(neighbourLat, lonMin + jj * delLon);
                                    neighbourDEM[r][c] = dem.getElevation(geoPos);
                                }
                                tileDEM[ii][jj] = neighbourDEM[r][c];
                            }
                            ++c;
                        }
                        ++r;
                    }
                    SARGeocoding.computeLocalIncidenceAngle(localGeometry, demNoDataValue, false, true, false, 0, 0, 0, 0, neighbourDEM, // in degrees
                    localIncidenceAngles);
                    if (localIncidenceAngles[1] == SARGeocoding.NonValidIncidenceAngle) {
                        continue;
                    }
                    final double v = computeBackscatteredPower(localIncidenceAngles[1]);
                    saveSimulatedData(posData.azimuthIndex, posData.rangeIndex, v, x0, y0, w, h, targetTile, masterBuffer);
                    int idx = 0;
                    if (saveDEM || saveLocalIncidenceAngle)
                        idx = targetTile.getDataBufferIndex((int) posData.rangeIndex, (int) posData.azimuthIndex);
                    if (saveDEM && idx >= 0) {
                        demBandBuffer.setElemDoubleAt(idx, alt);
                    }
                    if (saveZeroHeightSimulation) {
                        saveSimulatedData(posData.azimuthIndex, posData.rangeIndex, 1, x0, y0, w, h, targetTile, zeroHeightBandBuffer);
                    }
                    if (saveLocalIncidenceAngle && idx >= 0) {
                        localIncidenceAngleBandBuffer.setElemDoubleAt(idx, localIncidenceAngles[1]);
                    }
                    if (saveLayoverShadowMask) {
                        int rIndex = (int) posData.rangeIndex;
                        int aIndex = (int) posData.azimuthIndex;
                        if (rIndex >= x0 && rIndex < x0 + w && aIndex >= y0 && aIndex < y0 + h) {
                            azIndex[j] = posData.azimuthIndex;
                            rgIndex[j] = posData.rangeIndex;
                            slrs[j] = posData.slantRange;
                            elev[j] = computeElevationAngle(posData.slantRange, posData.earthPoint, posData.sensorPos);
                            savePixel[j] = true;
                        } else {
                            savePixel[j] = false;
                        }
                    }
                }
                if (saveLayoverShadowMask) {
                    computeLayoverShadow(x0, y0, w, h, savePixel, slrs, elev, azIndex, rgIndex, targetTile, layoverShadowMaskBuffer);
                }
            }
        } else {
            final int widthExt = xmax - xmin;
            final int heightExt = ymax - ymin;
            if (saveLayoverShadowMask) {
                slrs = new double[widthExt];
                elev = new double[widthExt];
                azIndex = new double[widthExt];
                rgIndex = new double[widthExt];
                savePixel = new boolean[widthExt];
            }
            final double[][] localDEM = new double[heightExt + 2][widthExt + 2];
            final TileGeoreferencing tileGeoRef = new TileGeoreferencing(targetProduct, xmin, ymin, widthExt, heightExt);
            if (saveZeroHeightSimulation) {
                for (double[] aLocalDEM : localDEM) {
                    Arrays.fill(aLocalDEM, 1);
                }
            } else {
                final boolean valid = DEMFactory.getLocalDEM(dem, demNoDataValue, demResamplingMethod, tileGeoRef, xmin, ymin, widthExt, heightExt, sourceProduct, true, localDEM);
                if (!valid)
                    return;
            }
            for (int y = ymin; y < ymax; y++) {
                final int yy = y - ymin;
                if (saveLayoverShadowMask) {
                    Arrays.fill(slrs, 0.0);
                    Arrays.fill(elev, 0.0);
                    Arrays.fill(azIndex, 0.0);
                    Arrays.fill(rgIndex, 0.0);
                    Arrays.fill(savePixel, Boolean.FALSE);
                }
                for (int x = xmin; x < xmax; x++) {
                    final int xx = x - xmin;
                    Double alt = localDEM[yy + 1][xx + 1];
                    if (alt.equals(demNoDataValue))
                        continue;
                    tileGeoRef.getGeoPos(x, y, geoPos);
                    if (!geoPos.isValid())
                        continue;
                    double lat = geoPos.lat;
                    double lon = geoPos.lon;
                    if (lon >= 180.0) {
                        lon -= 360.0;
                    }
                    if (orbitMethod) {
                        double[] latlon = jOrbit.lp2ell(new Point(x + 0.5, y + 0.5), meta);
                        lat = latlon[0] * Constants.RTOD;
                        lon = latlon[1] * Constants.RTOD;
                        alt = dem.getElevation(new GeoPos(lat, lon));
                    }
                    GeoUtils.geo2xyzWGS84(lat, lon, alt, posData.earthPoint);
                    if (!sarPosition.getPosition(posData))
                        continue;
                    final LocalGeometry localGeometry = new LocalGeometry(x, y, tileGeoRef, posData.earthPoint, posData.sensorPos);
                    final double[] localIncidenceAngles = { SARGeocoding.NonValidIncidenceAngle, SARGeocoding.NonValidIncidenceAngle };
                    SARGeocoding.computeLocalIncidenceAngle(localGeometry, demNoDataValue, false, true, false, xmin, ymin, x, y, localDEM, // in degrees
                    localIncidenceAngles);
                    if (localIncidenceAngles[1] == SARGeocoding.NonValidIncidenceAngle)
                        continue;
                    final double v = computeBackscatteredPower(localIncidenceAngles[1]);
                    saveSimulatedData(posData.azimuthIndex, posData.rangeIndex, v, x0, y0, w, h, targetTile, masterBuffer);
                    int idx = 0;
                    if (saveDEM || saveLocalIncidenceAngle)
                        idx = targetTile.getDataBufferIndex((int) posData.rangeIndex, (int) posData.azimuthIndex);
                    if (saveDEM && idx >= 0) {
                        demBandBuffer.setElemDoubleAt(idx, alt);
                    }
                    if (saveZeroHeightSimulation) {
                        saveSimulatedData(posData.azimuthIndex, posData.rangeIndex, 1, x0, y0, w, h, targetTile, zeroHeightBandBuffer);
                    }
                    if (saveLocalIncidenceAngle && idx >= 0) {
                        localIncidenceAngleBandBuffer.setElemDoubleAt(idx, localIncidenceAngles[1]);
                    }
                    if (saveLayoverShadowMask) {
                        int rIndex = (int) posData.rangeIndex;
                        int aIndex = (int) posData.azimuthIndex;
                        if (rIndex >= x0 && rIndex < x0 + w && aIndex >= y0 && aIndex < y0 + h) {
                            azIndex[xx] = posData.azimuthIndex;
                            rgIndex[xx] = posData.rangeIndex;
                            slrs[xx] = posData.slantRange;
                            elev[xx] = computeElevationAngle(posData.slantRange, posData.earthPoint, posData.sensorPos);
                            savePixel[xx] = true;
                        } else {
                            savePixel[xx] = false;
                        }
                    }
                }
                if (saveLayoverShadowMask) {
                    computeLayoverShadow(x0, y0, w, h, savePixel, slrs, elev, azIndex, rgIndex, targetTile, layoverShadowMaskBuffer);
                }
            }
        }
    } catch (Throwable e) {
        OperatorUtils.catchOperatorException(getId(), e);
    }
}

Example 3 with LocalGeometry

use of org.esa.snap.engine_utilities.eo.LocalGeometry in project s1tbx by senbox-org.

the class SARSimTerrainCorrectionOp method computeTileStack.

/**
 * Called by the framework in order to compute the stack of tiles for the given target bands.
 * <p>The default implementation throws a runtime exception with the message "not implemented".</p>
 *
 * @param targetTiles     The current tiles to be computed for each target band.
 * @param targetRectangle The area in pixel coordinates to be computed (same for all rasters in <code>targetRasters</code>).
 * @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 rasters.
 */
@Override
public void computeTileStack(Map<Band, Tile> targetTiles, Rectangle targetRectangle, ProgressMonitor pm) throws OperatorException {
    processingStarted = true;
    final int x0 = targetRectangle.x;
    final int y0 = targetRectangle.y;
    final int w = targetRectangle.width;
    final int h = targetRectangle.height;
    final int ymax = y0 + h;
    final int xmax = x0 + w;
    // System.out.println("x0 = " + x0 + ", y0 = " + y0 + ", w = " + w + ", h = " + h);
    final GeoPos geoPos = new GeoPos();
    final PosVector earthPoint = new PosVector();
    final PosVector sensorPos = new PosVector();
    final int srcMaxRange = sourceImageWidth - 1;
    final int srcMaxAzimuth = sourceImageHeight - 1;
    ProductData demBuffer = null;
    ProductData latBuffer = null;
    ProductData lonBuffer = null;
    ProductData localIncidenceAngleBuffer = null;
    ProductData projectedLocalIncidenceAngleBuffer = null;
    ProductData layoverShadowingMasksBuffer = null;
    ProductData incidenceAngleFromEllipsoidBuffer = null;
    final Set<Band> keySet = targetTiles.keySet();
    if (!warpDataAvailable) {
        getWarpData(keySet, targetRectangle);
        outputResidualAndShiftFiles();
    }
    final List<RangeDopplerGeocodingOp.TileData> trgTileList = new ArrayList<>();
    for (Band targetBand : keySet) {
        if (targetBand.getName().equals("elevation")) {
            demBuffer = targetTiles.get(targetBand).getDataBuffer();
            continue;
        }
        if (targetBand.getName().equals("latitude")) {
            latBuffer = targetTiles.get(targetBand).getDataBuffer();
            continue;
        }
        if (targetBand.getName().equals("longitude")) {
            lonBuffer = targetTiles.get(targetBand).getDataBuffer();
            continue;
        }
        if (targetBand.getName().equals("localIncidenceAngle")) {
            localIncidenceAngleBuffer = targetTiles.get(targetBand).getDataBuffer();
            continue;
        }
        if (targetBand.getName().equals("projectedLocalIncidenceAngle")) {
            projectedLocalIncidenceAngleBuffer = targetTiles.get(targetBand).getDataBuffer();
            continue;
        }
        if (targetBand.getName().equals(SARSimulationOp.layoverShadowMaskBandName)) {
            layoverShadowingMasksBuffer = targetTiles.get(targetBand).getDataBuffer();
            continue;
        }
        if (targetBand.getName().equals("incidenceAngleFromEllipsoid")) {
            incidenceAngleFromEllipsoidBuffer = targetTiles.get(targetBand).getDataBuffer();
            continue;
        }
        final String[] srcBandNames = targetBandNameToSourceBandName.get(targetBand.getName());
        final Band[] srcBands = new Band[] { sourceProduct.getBand(srcBandNames[0]), srcBandNames.length > 1 ? sourceProduct.getBand(srcBandNames[1]) : null };
        final RangeDopplerGeocodingOp.TileData td = new RangeDopplerGeocodingOp.TileData(targetTiles.get(targetBand), srcBands, isPolsar, outputComplex, targetBand.getName(), getBandUnit(targetBand.getName()), absRoot, calibrator, imgResampling, null, null);
        td.applyRadiometricNormalization = targetBandapplyRadiometricNormalizationFlag.get(targetBand.getName());
        td.applyRetroCalibration = targetBandApplyRetroCalibrationFlag.get(targetBand.getName());
        trgTileList.add(td);
    }
    final RangeDopplerGeocodingOp.TileData[] trgTiles = trgTileList.toArray(new RangeDopplerGeocodingOp.TileData[trgTileList.size()]);
    final TileGeoreferencing tileGeoRef = new TileGeoreferencing(targetProduct, x0 - 1, y0 - 1, w + 2, h + 2);
    int diffLat = Math.abs(latitude.getPixelInt(0, 0) - latitude.getPixelInt(0, targetImageHeight));
    try {
        final double[][] localDEM = new double[h + 2][w + 2];
        if (useAvgSceneHeight) {
            DEMFactory.fillDEM(localDEM, avgSceneHeight);
        } else {
            final boolean valid = DEMFactory.getLocalDEM(dem, demNoDataValue, demResamplingMethod, tileGeoRef, x0, y0, w, h, sourceProduct, true, localDEM);
            if (!valid) {
                return;
            }
        }
        for (int y = y0; y < ymax; y++) {
            final int yy = y - y0 + 1;
            for (int x = x0; x < xmax; x++) {
                final int index = trgTiles[0].targetTile.getDataBufferIndex(x, y);
                final Double alt = localDEM[yy][x - x0 + 1];
                if (!useAvgSceneHeight && alt.equals(demNoDataValue)) {
                    if (saveDEM) {
                        demBuffer.setElemDoubleAt(index, demNoDataValue);
                    }
                    continue;
                }
                tileGeoRef.getGeoPos(x, y, geoPos);
                if (!geoPos.isValid()) {
                    continue;
                }
                final double lat = geoPos.lat;
                double lon = geoPos.lon;
                if (lon >= 180.0) {
                    lon -= 360.0;
                }
                GeoUtils.geo2xyzWGS84(lat, lon, alt, earthPoint);
                final double zeroDopplerTime = getEarthPointZeroDopplerTime(earthPoint);
                if (Double.compare(zeroDopplerTime, NonValidZeroDopplerTime) == 0) {
                    if (saveDEM) {
                        demBuffer.setElemDoubleAt(index, demNoDataValue);
                    }
                    continue;
                }
                double slantRange = SARGeocoding.computeSlantRange(zeroDopplerTime, orbit, earthPoint, sensorPos);
                double zeroDoppler = zeroDopplerTime;
                if (!skipBistaticCorrection) {
                    // skip bistatic correction for COSMO, TerraSAR-X and RadarSAT-2
                    zeroDoppler = zeroDopplerTime + slantRange / Constants.lightSpeedInMetersPerDay;
                    slantRange = SARGeocoding.computeSlantRange(zeroDoppler, orbit, earthPoint, sensorPos);
                }
                final double azimuthIndex = (zeroDoppler - firstLineUTC) / lineTimeInterval;
                double rangeIndex = SARGeocoding.computeRangeIndex(srgrFlag, sourceImageWidth, firstLineUTC, lastLineUTC, rangeSpacing, zeroDoppler, slantRange, nearEdgeSlantRange, srgrConvParams);
                // temp fix for descending Radarsat2
                if (!nearRangeOnLeft) {
                    rangeIndex = srcMaxRange - rangeIndex;
                }
                if (!SARGeocoding.isValidCell(rangeIndex, azimuthIndex, lat, lon, diffLat, sourceProduct.getSceneGeoCoding(), srcMaxRange, srcMaxAzimuth, sensorPos)) {
                    if (saveDEM) {
                        demBuffer.setElemDoubleAt(index, demNoDataValue);
                    }
                } else {
                    final double[] localIncidenceAngles = { SARGeocoding.NonValidIncidenceAngle, SARGeocoding.NonValidIncidenceAngle };
                    if (saveLocalIncidenceAngle || saveProjectedLocalIncidenceAngle || saveSigmaNought) {
                        final LocalGeometry localGeometry = new LocalGeometry(x, y, tileGeoRef, earthPoint, sensorPos);
                        SARGeocoding.computeLocalIncidenceAngle(localGeometry, demNoDataValue, saveLocalIncidenceAngle, saveProjectedLocalIncidenceAngle, saveSigmaNought, x0, y0, x, y, localDEM, // in degrees
                        localIncidenceAngles);
                        if (saveLocalIncidenceAngle && localIncidenceAngles[0] != SARGeocoding.NonValidIncidenceAngle) {
                            localIncidenceAngleBuffer.setElemDoubleAt(index, localIncidenceAngles[0]);
                        }
                        if (saveProjectedLocalIncidenceAngle && localIncidenceAngles[1] != SARGeocoding.NonValidIncidenceAngle) {
                            projectedLocalIncidenceAngleBuffer.setElemDoubleAt(index, localIncidenceAngles[1]);
                        }
                    }
                    if (saveDEM) {
                        demBuffer.setElemDoubleAt(index, alt);
                    }
                    if (saveLatLon) {
                        latBuffer.setElemDoubleAt(index, lat);
                        lonBuffer.setElemDoubleAt(index, lon);
                    }
                    if (saveLayoverShadowMask) {
                        final Rectangle srcRect = new Rectangle((int) (rangeIndex + 0.5), (int) (azimuthIndex + 0.5), 1, 1);
                        final Tile sourceTile = getSourceTile(maskBand, srcRect);
                        final int m = sourceTile.getDataBuffer().getElemIntAt(sourceTile.getDataBufferIndex((int) (rangeIndex + 0.5), (int) (azimuthIndex + 0.5)));
                        layoverShadowingMasksBuffer.setElemIntAt(index, m);
                    }
                    if (saveIncidenceAngleFromEllipsoid) {
                        incidenceAngleFromEllipsoidBuffer.setElemDoubleAt(index, incidenceAngle.getPixelDouble(rangeIndex, azimuthIndex));
                    }
                    for (RangeDopplerGeocodingOp.TileData tileData : trgTiles) {
                        final Unit.UnitType bandUnit = getBandUnit(tileData.bandName);
                        final String[] srcBandName = targetBandNameToSourceBandName.get(tileData.bandName);
                        final Band srcBand = sourceProduct.getBand(srcBandName[0]);
                        final PixelPos pixelPos = new PixelPos();
                        final WarpData warpData = warpDataMap.get(srcBand);
                        if (!warpData.isValid()) {
                            continue;
                        }
                        warpData.getWarpedCoords(warpPolynomialOrder, rangeIndex, azimuthIndex, pixelPos);
                        if (pixelPos.x < 0.0 || pixelPos.x >= srcMaxRange || pixelPos.y < 0.0 || pixelPos.y >= srcMaxAzimuth) {
                            tileData.tileDataBuffer.setElemDoubleAt(index, tileData.noDataValue);
                            continue;
                        }
                        final int[] subSwathIndex = { INVALID_SUB_SWATH_INDEX };
                        double v = getPixelValue(pixelPos.y, pixelPos.x, tileData, subSwathIndex);
                        if (v != tileData.noDataValue && tileData.applyRadiometricNormalization) {
                            if (localIncidenceAngles[1] != SARGeocoding.NonValidIncidenceAngle) {
                                final double satelliteHeight = Math.sqrt(sensorPos.x * sensorPos.x + sensorPos.y * sensorPos.y + sensorPos.z * sensorPos.z);
                                final double sceneToEarthCentre = Math.sqrt(earthPoint.x * earthPoint.x + earthPoint.y * earthPoint.y + earthPoint.z * earthPoint.z);
                                v = calibrator.applyCalibration(v, rangeIndex, azimuthIndex, slantRange, satelliteHeight, sceneToEarthCentre, localIncidenceAngles[1], tileData.bandName, tileData.bandPolar, bandUnit, // use projected incidence angle
                                subSwathIndex);
                            } else {
                                v = tileData.noDataValue;
                            }
                        }
                        tileData.tileDataBuffer.setElemDoubleAt(index, v);
                    }
                    orthoDataProduced = true;
                }
            }
        }
    } catch (Throwable e) {
        // to prevent multiple error messages
        orthoDataProduced = true;
        OperatorUtils.catchOperatorException(getId(), e);
    }
}