package org.n52.math;
import static java.lang.Math.*;
import java.net.URL;
import Jama.Matrix;
/**
* Supplier class for azimuth, zenith angle, distance and pixel size algorithms,
* to calculate with satellite triangulation.
* SatelliteTriangulation holds the geometric properties (3D) of the triangle ETS:
*
* E | earth-center = (0.0.0),
* |
T | observed terrain-pixel position,
* |
S | satellite position;
* |
* All positions are metric, cartesian and geo-centric.
* T depends on assumptions about ellipsoid and image georeference.
* terrTriangulation contains the geometry depending merely on T.
* S is for the moment geo-stationary (latit = 0, longit = constant).
* The purpose is to produce (output) parameters for (a.o) solar reflection:
*
* - {@linkplain #getZenithAngleAlgorithm zenith angle},
* angle between local vertical and satellite direction;
*
- {@linkplain #getAzimuthAlgorithm azimuth},
* angle between local North and ground-projected satellite direction;
*
- {@linkplain #getDistanceAlgorithm distance},
* distance from terrain point to satellite (meters);
*
- {@linkplain #getPixelSizeAlgorithm pixel size},
* apparent area per pixel at location in square km ;
*
* @author Jan Hendrikse
*/
public class SatelliteTriangulation {
//protected final double EPS10 = 0.0000000001;
private double lat = Double.NaN;
private double lon = Double.NaN;
// private double phi;
// private double lam;
private final double satPsi;
private final double satLam;
private final double pixelSizeAtNadir;
private final double satOrbitRadius;
private final double satHeight;
private Vector3D satellitePosition = null;
private TerrainTriangulation terrTriangulation;
private double zenithAngle;
private double azimuth;
private double distance;
private double pixelSize;
private boolean zenithAngleCalculated;
private boolean azimuthCalculated;
private boolean distanceCalculated;
private boolean pixelSizeCalculated;
private final static String[] paramNames = {
"lat", "lon"
};
private final static String[] paramDescr = {
"latitude in degrees at terrain location",
"longitude in degrees at terrain location"
};
/**
* constructor of SatelliteTriangulation using terrain and satellite -related
* geometric input to allow mainly angle computations
* satHeight is a final member depending on satOrbitRadius and
* the eq radius a of the assumed ellipsoid, for the time being
* @param terrTri holds terrainlocation specific vectors
* @param satLat geocentric latitude (angle with equator-plane)
* @param satLon (subSatellite) longitude
* @param pixNadir
* @param satOrbitRadius
*/
//TODO Eventually satOrbitRadius can be replaced by (time-dependent)
// orbit parameters (possibly partly entered by the user )
public SatelliteTriangulation(TerrainTriangulation terrTri,
double satLat, double satLon,
double pixNadir, double satOrbitRadius) {
this.terrTriangulation = terrTri;
this.satPsi = satLat * PI / 180;
this.satLam = satLon * PI / 180;
this.pixelSizeAtNadir = pixNadir;
this.satOrbitRadius = satOrbitRadius;
this.satHeight = satOrbitRadius - terrTri.getEquatorRadius();
clear();
}
/**
* @return vector from EarthCenter to Satellite center
* defined by polar coordinates: posRadius (length of vec)
* satLam (angle between vec and prime meridian plane)
* satPsi (angle between vec and equator plane), i.e. geocentric,
* not geodetic (phi) latitude
*/
private Vector3D getSatellitePositionVector() {
if (null == satellitePosition) {
Vector3D vec = new Vector3D();
double posRadius = satOrbitRadius;
vec.set(0, posRadius*cos(satLam)*cos(satPsi));
vec.set(1, posRadius*sin(satLam)*cos(satPsi));
vec.set(2, posRadius*sin(satPsi));
satellitePosition = vec;
}
return satellitePosition;
}
private Vector3D getTerrainToSatelliteVector(Vector3D vecTerrainPosition,
Vector3D vecSatellitePosition) {
return vecSatellitePosition.minus(vecTerrainPosition);
}
private double getZenithToSatelliteAngle(Vector3D vecTerrVerticalUnitVector,
Vector3D vecTerrToSatellite) {
double norm = vecTerrToSatellite.normF();
Vector3D U2 = vecTerrToSatellite.times(1/norm);
return acos(vecTerrVerticalUnitVector.dotProduct(U2));
}
private double getDistanceToSatellite(Vector3D vecTerrToSatellite) {
double norm = vecTerrToSatellite.normF();
return norm;
}
private void check(double lat, double lon) {
if ((this.lat != lat ) || (this.lon != lon)) {
clear();
this.lat = lat;
this.lon = lon;
// this.phi = lat * PI / 180;
// this.lam = lon * PI / 180;
}
return;
}
private void clear() {
zenithAngleCalculated = false;
azimuthCalculated = false;
distanceCalculated = false;
pixelSizeCalculated = false;
}
private double getZenithAngle(double lat, double lon) {
check(lat,lon);
if (!zenithAngleCalculated) {
Vector3D unitV1 = terrTriangulation.getUnitVerticalVector(lat,lon);
Vector3D posVecTerr =
terrTriangulation.getTerrainPositionVector(lat,lon);
Vector3D posVecSat = getSatellitePositionVector();
Vector3D V2 = getTerrainToSatelliteVector(posVecTerr, posVecSat);
Vector3D unitV2 = V2.times(1/V2.normF()); //reduce to unit length
zenithAngleCalculated = true;
zenithAngle = acos(unitV1.dotProduct(unitV2)) * 180.0 / PI;
}
return zenithAngle;
}
private double getAzimuth(double lat, double lon) {
check(lat,lon);
if (!azimuthCalculated) {
Vector3D unitNorth = terrTriangulation.getUnitNorthVector(lat,lon);
Vector3D unitEast = terrTriangulation.getUnitEastVector(lat,lon);
Vector3D posVecTerr = terrTriangulation.getTerrainPositionVector(lat,lon);
Vector3D posVecSat = getSatellitePositionVector();
Vector3D V2 = getTerrainToSatelliteVector(posVecTerr, posVecSat);
Vector3D unitV2 = V2.times(1/V2.normF()); //reduce to unit length
///direction-cosines of PS vector with local east- and north-axes:
double nort = unitNorth.dotProduct(unitV2);
double east = unitEast.dotProduct(unitV2);
azimuth = atan2(east,nort) * 180.0 / PI;
if (azimuth < 0)
azimuth += 360;
azimuthCalculated = true;
}
return azimuth;
}
private double getDistance(double lat, double lon) {
check(lat,lon);
if (!distanceCalculated) {
Vector3D posVecTerr = terrTriangulation.getTerrainPositionVector(lat,lon);
Vector3D posVecSat =
getSatellitePositionVector();
Vector3D V2 = getTerrainToSatelliteVector(posVecTerr, posVecSat);
distanceCalculated = true;
distance = V2.length();
}
return distance;
}
private double getPixelSize(double lat, double lon) {
check(lat,lon);
if (!pixelSizeCalculated) {
Vector3D unitV1 = terrTriangulation.getUnitVerticalVector(lat,lon);
Vector3D posVecTerr = terrTriangulation.getTerrainPositionVector(lat,lon);
Vector3D posVecSat = getSatellitePositionVector();
Vector3D V2 = getTerrainToSatelliteVector(posVecTerr, posVecSat);
double dist = V2.length();
Vector3D unitV2 = V2.times(1/V2.normF()); //reduce to unit length
double cosZenithAngle = unitV1.dotProduct(unitV2);
if (cosZenithAngle < 0.0001)
pixelSize = Double.NaN;
else // in square km
pixelSize = pixelSizeAtNadir*dist/cosZenithAngle/satHeight;
}
return pixelSize;
}
/**
* returns the ZenithAngle Algorithm; it produces the zenith angle,
* angle (in degrees) between local vertical and
* satellite direction, given local lat and lon (in degrees)
* @return the Algorithm
*/
public Algorithm getZenithAngleAlgorithm() {
return new AbstractAlgorithm("Zenith Angle",
"angle between local vertical and satellite direction",
paramNames, paramDescr) {
public double calculate(double[] params) {
return getZenithAngle(params[0], params[1]);
}
};
}
/**
* returns the Azimuth Algorithm; it produces the azimuth,
* angle (in degrees) between local North and projected
* satellite direction, given local lat and lon (in degrees)
* @return the Algorithm
*/
public Algorithm getAzimuthAlgorithm() {
return new AbstractAlgorithm("Azimuth",
"angle between local North and projected " +
"satellite direction",
paramNames, paramDescr) {
public double calculate(double[] params) {
return getAzimuth(params[0], params[1]);
}
};
}
/**
* returns the Distance Algorithm; it produces the distance,
* between local position and satellite position, in meters,
* given local lat and lon (in degrees)
* @return the Algorithm
*/
public Algorithm getDistanceAlgorithm() {
return new AbstractAlgorithm("Distance",
"distance from terrain point to satellite",
paramNames, paramDescr) {
public double calculate(double[] params) {
return getDistance(params[0], params[1]);
}
};
}
/**
* returns the Pixel Size Algorithm; it produces the pixel size,
* in square km (area covered by one pixel)at a local position,
* given local lat and lon (in degrees) and nadir pixel size
* @return the Algorithm
*/
public Algorithm getPixelSizeAlgorithm() {
return new AbstractAlgorithm("Pixel Size",
"area per pixel at location in square km",
paramNames, paramDescr) {
public double calculate(double[] params) {
return getPixelSize(params[0], params[1]);
}
};
}
}