/*
* geoProj.c --
*
* This file defines functions that convert back and forth between
* geopoints and map-points. See the user documentation for more
* information.
*
* Copyright (c) 2004 Gordon D. Carrie. All rights reserved.
*
* Licensed under the Open Software License version 2.1
*
* Please address questions and feedback to user0@tkgeomap.org
*
* @(#) $Id: geoProj.c,v 1.9 2007/06/27 18:38:56 tkgeomap Exp $
*
********************************************
*
*/
/*
* Formulas and notation for most projections are from:
* Snyder, John P.
* Map Projections used by the U.S. Geological Survey.
* (Geological Survey bulletin ; 1532)
* United States Government Printing Office, Washington: 1982.
*/
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <assert.h>
#include <tcl.h>
#include "geography.h"
#include "geoProj.h"
/*
* Local functions.
*/
static MapPt latLonToCylEqDist _ANSI_ARGS_((GeoPt geoPt,
struct GeoProj *projPtr));
static GeoPt cylEqDistToLatLon _ANSI_ARGS_((MapPt mapPt,
struct GeoProj *projPtr));
static GeoProjInfo cylEqDistInfo _ANSI_ARGS_((struct GeoProj *projPtr));
static MapPt latLonToCylEqArea _ANSI_ARGS_((GeoPt geoPt,
struct GeoProj *projPtr));
static GeoPt cylEqAreaToLatLon _ANSI_ARGS_((MapPt mapPt,
struct GeoProj *projPtr));
static MapPt latLonToMercator _ANSI_ARGS_((GeoPt geoPt,
struct GeoProj *projPtr));
static GeoPt mercatorToLatLon _ANSI_ARGS_((MapPt mapPt,
struct GeoProj *projPtr));
static MapPt latLonToLambertConfConic _ANSI_ARGS_((GeoPt geoPt,
struct GeoProj *projPtr));
static GeoPt lambrtConfConicToLatLon _ANSI_ARGS_((MapPt mapPt,
struct GeoProj *projPtr));
static GeoProjInfo lambrtConfConicInfo _ANSI_ARGS_((
struct GeoProj *projPtr));
static MapPt latLonToLambrtEqArea _ANSI_ARGS_((GeoPt geoPt,
struct GeoProj *projPtr));
static GeoPt lambrtEqAreaToLatLon _ANSI_ARGS_((MapPt mapPt,
struct GeoProj *projPtr));
static MapPt latLonToStereographic _ANSI_ARGS_((GeoPt geoPt,
struct GeoProj *projPtr));
static GeoPt stereographicToLatLon _ANSI_ARGS_((MapPt mapPt,
struct GeoProj *projPtr));
static MapPt latLonToOrthographic _ANSI_ARGS_((GeoPt geoPt,
struct GeoProj *projPtr));
static GeoPt orthographicToLatLon _ANSI_ARGS_((MapPt mapPt,
struct GeoProj *projPtr));
static GeoProjInfo refLonProjInfo _ANSI_ARGS_((struct GeoProj *projPtr));
static GeoProjInfo refPtProjInfo _ANSI_ARGS_((struct GeoProj *projPtr));
/*
*----------------------------------------------------------------------
*
* GeoProjAlloc --
*
* This procedure allocates a new GeoProj structure.
*
* Results:
* See the user documentation.
*
* Side effects:
* A GeoProj structure is allocated and initialized.
* The memory should eventually be freed with a call to GeoProjDestroy.
*
*----------------------------------------------------------------------
*/
GeoProj
GeoProjAlloc()
{
struct GeoProj *projPtr;
projPtr = (GeoProj )CKALLOC(sizeof(*projPtr));
GeoProjInit(projPtr);
return projPtr;
}
/*
*----------------------------------------------------------------------
*
* GeoProjInit --
*
* This procedure initializes a GeoProj structure.
*
* Results:
* None.
*
* Side effects:
* Modifies the member of a GeoProj structure. Previous contents,
* assumed to be garbage, are overwritten.
*
*----------------------------------------------------------------------
*/
void
GeoProjInit(projPtr)
struct GeoProj *projPtr;
{
projPtr->params = NULL;
SetCylEqDist(projPtr, 0, 0);
GeoProjSetRotation(projPtr, AngleFmDeg(0.0));
}
/*
*----------------------------------------------------------------------
*
* GeoProjFree --
*
* This procedure frees the internal allocations in a GeoProj structure.
*
* Results:
* None.
*
* Side effects:
* See the user documentation.
*
*----------------------------------------------------------------------
*/
void
GeoProjFree(projPtr)
struct GeoProj *projPtr;
{
if( projPtr ) {
CKFREE((char *)projPtr->params);
}
}
/*
*----------------------------------------------------------------------
*
* GeoProjDestroy --
*
* This procedure frees internal storage of a GeoProj structure and
* frees the structure.
*
* Results:
* None.
*
* Side effects:
* See the user documentation.
*
*----------------------------------------------------------------------
*/
void
GeoProjDestroy(projPtr)
struct GeoProj *projPtr;
{
GeoProjFree(projPtr);
CKFREE((char *)projPtr);
}
/*
*----------------------------------------------------------------------
*
* ProjToLatLon --
*
* This procedure converts a map-point to a geopoint.
*
* Results:
* See the user documentation.
*
* Side effects:
* See the user documentation.
*
*----------------------------------------------------------------------
*/
GeoPt
ProjToLatLon(mapPt, projPtr)
MapPt mapPt;
struct GeoProj *projPtr;
{
if ( MapPtIsNowhere(mapPt) ) {
return GeoPtNowhere();
}
if (projPtr->rotation != 0) {
double abs, ord;
abs = mapPt.abs * projPtr->cosr - mapPt.ord * projPtr->sinr;
ord = mapPt.abs * projPtr->sinr + mapPt.ord * projPtr->cosr;
mapPt.abs = abs;
mapPt.ord = ord;
}
return projPtr->projToLatLonProc(mapPt, projPtr);
}
/*
*----------------------------------------------------------------------
*
* LatLonToProj --
*
* This procedure converts a geopoint to a map-point.
*
* Results:
* See the user documentation.
*
* Side effects:
* See the user documentation.
*
*----------------------------------------------------------------------
*/
MapPt
LatLonToProj(geoPt, projPtr)
GeoPt geoPt;
struct GeoProj *projPtr;
{
MapPt mapPt;
if ( GeoPtIsNowhere(geoPt) ) {
return MapPtNowhere();
}
mapPt = projPtr->latLonToProjProc(geoPt, projPtr);
if (projPtr->rotation != 0) {
double abs, ord;
abs = mapPt.abs * projPtr->cosr + mapPt.ord * projPtr->sinr;
ord = mapPt.ord * projPtr->cosr - mapPt.abs * projPtr->sinr;
mapPt.abs = abs;
mapPt.ord = ord;
}
return mapPt;
}
/*
*----------------------------------------------------------------------
*
* GeoProjGetInfo --
*
* This is a member access function.
*
* Results:
* See the user documentation.
*
* Side effects:
* See the user documentation.
*
*----------------------------------------------------------------------
*/
GeoProjInfo
GeoProjGetInfo(projPtr)
struct GeoProj *projPtr;
{
return projPtr->infoProc(projPtr);
}
/*
*----------------------------------------------------------------------
*
* GeoProjGetType --
*
* This is a member access function.
*
* Results:
* See the user documentation.
*
* Side effects:
* See the user documentation.
*
*----------------------------------------------------------------------
*/
enum ProjType
GeoProjGetType(projPtr)
struct GeoProj *projPtr;
{
return projPtr->type;
}
/*
*----------------------------------------------------------------------
*
* GeoProjDescriptor --
*
* This is a member access function.
*
* Results:
* See the user documentation.
*
* Side effects:
* See the user documentation.
*
*----------------------------------------------------------------------
*/
char *
GeoProjDescriptor(projPtr)
struct GeoProj *projPtr;
{
return projPtr->descriptor;
}
/*
*----------------------------------------------------------------------
*
* GeoProjSetRotation --
*
* This function determines how much a map is rotated.
*
* Results:
* None.
*
* Side effects:
* The rotation angle an associated members are set in the structure.
*
*----------------------------------------------------------------------
*/
void
GeoProjSetRotation(projPtr, angle)
struct GeoProj *projPtr; /* Projection to modify */
Angle angle; /* Projection angle */
{
projPtr->rotation = angle;
projPtr->cosr = ICos(angle);
projPtr->sinr = ISin(angle);
}
/*
*----------------------------------------------------------------------
*
* SetCylEqDist --
*
* This procedure sets a projection to Cylindrical Equidistant.
*
* Results:
* None.
*
* Side effects:
* The members of a GeoProj structure are modified.
*
*----------------------------------------------------------------------
*/
void
SetCylEqDist(projPtr, refLat, refLon)
struct GeoProj *projPtr; /* Projection to modify */
Angle refLat; /* Reference latitude */
Angle refLon; /* Reference longitude */
{
GeoProjCylEqDistParams *params;
projPtr->type = CylEqDist;
params = (GeoProjCylEqDistParams *)CKALLOC(sizeof(GeoProjCylEqDistParams));
params->refLat = DomainLat(refLat);
params->cosRLat = ICos(refLat);
params->refLon = GwchLon(refLon);
if (projPtr->params) {
CKFREE(projPtr->params);
}
projPtr->params = params;
projPtr->latLonToProjProc = latLonToCylEqDist;
projPtr->projToLatLonProc = cylEqDistToLatLon;
sprintf(projPtr->descriptor, "CylEqDist {%9.3f %-9.3f}",
AngleToDeg(params->refLat), AngleToDeg(params->refLon));
projPtr->infoProc = cylEqDistInfo;
}
/*
*----------------------------------------------------------------------
*
* cylEqDistInfo --
*
* This is a value for an infoProc member in a GeoProj structure.
* It is usually called by the GeoProjGetInfo procedure
*
* Results:
* Projection information is copied into a GeoProjInfo structure.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoProjInfo
cylEqDistInfo(projPtr)
struct GeoProj *projPtr;
{
GeoProjInfo geoProjInfo;
GeoProjCylEqDistParams *params = projPtr->params;
geoProjInfo.type = projPtr->type;
geoProjInfo.prm.refPt.lat = params->refLat;
geoProjInfo.prm.refPt.lon = params->refLon;
geoProjInfo.rotation = projPtr->rotation;
return geoProjInfo;
}
/*
*----------------------------------------------------------------------
*
* latLonToCylEqDist --
*
* This is a value for a latLonToProjProc member of a GeoProj structure.
* It is usually called by the LatLonToProj procedure.
*
* Results:
* A map-point corresponding to a geopoint.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static MapPt
latLonToCylEqDist(geoPt, projPtr)
GeoPt geoPt;
struct GeoProj *projPtr;
{
MapPt mapPt;
double lat, lon;
GeoProjCylEqDistParams *params = projPtr->params;
GeoPtGetDeg(DomainLonPt(geoPt, params->refLon), &lat, &lon);
mapPt.ord = MPERDEG * lat;
mapPt.abs = MPERDEG * params->cosRLat * lon;
return mapPt;
}
/*
*----------------------------------------------------------------------
*
* cylEqDistToLatLon --
*
* This is a value for a projToLatLonProc member of a GeoProj structure.
* It is usually called by the ProjToLatLon procedure.
*
* Results:
* A geopoint corresponding to a map-point.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
cylEqDistToLatLon(mapPt, projPtr)
MapPt mapPt;
struct GeoProj *projPtr;
{
GeoPt geoPt;
GeoProjCylEqDistParams *params = projPtr->params;
geoPt = GeoPtFmDeg(DEGPERM * mapPt.ord,
DEGPERM * mapPt.abs / params->cosRLat);
return DomainLonPt(geoPt, params->refLon);
}
/*
*----------------------------------------------------------------------
*
* SetMercator --
*
* This procedure sets a projection to Mercator.
*
* Results:
* None.
*
* Side effects:
* The members of a GeoProj structure are modified.
*
*----------------------------------------------------------------------
*/
void
SetMercator(projPtr, refLon)
struct GeoProj *projPtr; /* Projection to modify */
Angle refLon; /* Reference longitude */
{
Angle *refLonPtr;
refLonPtr = (Angle *)CKALLOC(sizeof(Angle));
projPtr->type = Mercator;
*refLonPtr = GwchLon(refLon);
if (projPtr->params) {
CKFREE(projPtr->params);
}
projPtr->params = refLonPtr;
projPtr->latLonToProjProc = latLonToMercator;
projPtr->projToLatLonProc = mercatorToLatLon;
sprintf(projPtr->descriptor, "Mercator %-9.3f", AngleToDeg(*refLonPtr));
projPtr->infoProc = refLonProjInfo;
}
/*
*----------------------------------------------------------------------
*
* latLonToMercator --
*
* This is a value for a latLonToProjProc member of a GeoProj structure.
* It is usually called by the LatLonToProj procedure.
*
* Results:
* A map-point corresponding to a geopoint.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static MapPt
latLonToMercator(geoPt, projPtr)
GeoPt geoPt;
struct GeoProj *projPtr;
{
Angle *refLon = projPtr->params;
double lat, lon;
MapPt mapPt;
Angle d90 = AngleFmDeg(90.0);
/*
* Undefined at pole
*/
if (LatCmp(geoPt.lat, d90) == Eq || LatCmp(geoPt.lat, -d90) == Eq) {
return MapPtNowhere();
}
GeoPtGetRad(DomainLonPt(geoPt, *refLon), &lat, &lon);
mapPt.ord = REarth() * log(tan(M_PI_4 + 0.5 * lat));
mapPt.abs = REarth() * lon;
return mapPt;
}
/*
*----------------------------------------------------------------------
*
* mercatorToLatLon --
*
* This is a value for a projToLatLonProc member of a GeoProj structure.
* It is usually called by the ProjToLatLon procedure.
*
* Results:
* A geopoint corresponding to a map-point.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
mercatorToLatLon(mapPt, projPtr)
MapPt mapPt;
struct GeoProj *projPtr;
{
double lat, lon; /* Latitude and longitude in radians */
Angle *refLon = projPtr->params;
lat = 2.0 * (atan(exp(mapPt.ord / REarth())) - M_PI_4);
lon = mapPt.abs / REarth();
return DomainLonPt(GeoPtFmRad(lat, lon), *refLon);
}
/*
*----------------------------------------------------------------------
*
* SetCylEqArea --
*
* This procedure sets a projection to Cylindrical Equal Area.
*
* Results:
* None.
*
* Side effects:
* The members of a GeoProj structure are modified.
*
*----------------------------------------------------------------------
*/
void
SetCylEqArea(projPtr, refLon)
struct GeoProj *projPtr; /* Projection to modify */
Angle refLon; /* Reference longitude */
{
Angle *refLonPtr;
refLonPtr = (Angle *)CKALLOC(sizeof(Angle));
projPtr->type = CylEqArea;
*refLonPtr = GwchLon(refLon);
projPtr->latLonToProjProc = latLonToCylEqArea;
projPtr->projToLatLonProc = cylEqAreaToLatLon;
sprintf(projPtr->descriptor, "CylEqArea %-9.3f", AngleToDeg(*refLonPtr));
projPtr->infoProc = refLonProjInfo;
}
/*
*----------------------------------------------------------------------
*
* latLonToCylEqArea --
*
* This is a value for a latLonToProjProc member of a GeoProj structure.
* It is usually called by the LatLonToProj procedure.
*
* Results:
* A map-point corresponding to a geopoint.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static MapPt
latLonToCylEqArea(geoPt, projPtr)
GeoPt geoPt;
struct GeoProj *projPtr;
{
MapPt mapPt;
Angle *refLon = projPtr->params;
double lat, lon;
GeoPtGetRad(DomainLonPt(geoPt, *refLon), &lat, &lon);
mapPt.ord = REarth() * sin(lat);
mapPt.abs = REarth() * lon;
return mapPt;
}
/*
*----------------------------------------------------------------------
*
* cylEqAreaToLatLon --
*
* This is a value for a projToLatLonProc member of a GeoProj structure.
* It is usually called by the ProjToLatLon procedure.
*
* Results:
* A geopoint corresponding to a map-point.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
cylEqAreaToLatLon(mapPt, projPtr)
MapPt mapPt;
struct GeoProj *projPtr;
{
Angle *refLon = projPtr->params;
double r, lat, lon;
r = mapPt.ord / REarth();
if (fabs(r) > 1.0) {
return GeoPtNowhere();
}
lat = asin(r);
lon = mapPt.abs / REarth();
return DomainLonPt(GeoPtFmRad(lat, lon), *refLon);
}
/*
*----------------------------------------------------------------------
*
* SetLambertConfConic --
*
* This procedure sets a projection to Lambert Conformal Conic.
*
* Results:
* None.
*
* Side effects:
* The members of a GeoProj structure are modified.
*
*----------------------------------------------------------------------
*/
void
SetLambertConfConic(projPtr, refLat, refLon)
struct GeoProj *projPtr; /* Projection to modify */
Angle refLat; /* Reference latitude */
Angle refLon; /* Reference longitude */
{
/*
* Notation from Snyder, p. 105
*/
double phi0, n, tan_n, RF;
GeoProjLambertConfConicParams *params;
params = (GeoProjLambertConfConicParams *)
CKALLOC(sizeof(GeoProjLambertConfConicParams));
refLat = DomainLat(refLat);
if (AngleCmp(refLat, 0) == 0) {
SetMercator(projPtr, refLat);
return;
}
phi0 = AngleToRad(refLat);
n = sin(phi0);
tan_n = pow(tan(M_PI_4 + 0.5 * phi0), n);
RF = REarth() * cos(phi0) * tan_n / n;
projPtr->type = LambertConfConic;
params->refLat = refLat;
params->refLon = GwchLon(refLon);
params->n = n;
params->RF = RF;
params->rho0 = REarth() / tan(phi0);
if (projPtr->params) {
CKFREE(projPtr->params);
}
projPtr->params = params;
projPtr->projToLatLonProc = lambrtConfConicToLatLon;
projPtr->latLonToProjProc = latLonToLambertConfConic;
sprintf(projPtr->descriptor, "LambertConfConic {%9.3f %-9.3f}",
AngleToDeg(params->refLat), AngleToDeg(params->refLon));
projPtr->infoProc = lambrtConfConicInfo;
}
/*
*----------------------------------------------------------------------
*
* lambrtConfConicInfo --
*
* This is a value for an infoProc member in a GeoProj structure.
* It is usually called by the GeoProjGetInfo procedure
*
* Results:
* Projection information is copied into a GeoProjInfo structure.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoProjInfo
lambrtConfConicInfo(projPtr)
struct GeoProj *projPtr;
{
GeoProjInfo geoProjInfo;
GeoProjLambertConfConicParams *params = projPtr->params;
geoProjInfo.type = projPtr->type;
geoProjInfo.prm.refPt.lat = params->refLat;
geoProjInfo.prm.refPt.lon = params->refLon;
geoProjInfo.rotation = projPtr->rotation;
return geoProjInfo;
}
/*
*----------------------------------------------------------------------
*
* latLonToLambertConfConic --
*
* This is a value for a latLonToProjProc member of a GeoProj structure.
* It is usually called by the LatLonToProj procedure.
*
* Results:
* A map-point corresponding to a geopoint.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static MapPt
latLonToLambertConfConic(geoPt, projPtr)
GeoPt geoPt;
struct GeoProj *projPtr;
{
/*
* Formulas and notation from Snyder, p. 105
*/
MapPt mapPt;
double theta, rho;
GeoProjLambertConfConicParams *params = projPtr->params;
double lat, lon;
Angle d90 = AngleFmDeg(90.0);
GeoPtGetRad(DomainLonPt(geoPt, params->refLon), &lat, &lon);
if (AngleCmp(geoPt.lat, d90) == 0) {
/*
* Point at North Pole
*/
if (AngleCmp(params->refLat, 0) == South) {
return MapPtNowhere();
}
rho = 0.0;
} else if (AngleCmp(geoPt.lat, -d90) == 0) {
/*
* Point at South Pole
*/
if (AngleCmp(params->refLat, 0) == North) {
return MapPtNowhere();
}
rho = 0.0;
} else {
/*
* Point between poles
*/
rho = params->RF / pow(tan(M_PI_4 + 0.5 * lat), params->n);
}
theta = params->n * (lon - AngleToRad(params->refLon));
mapPt.abs = rho * sin(theta);
mapPt.ord = params->rho0 - rho * cos(theta);
return mapPt;
}
/*
*----------------------------------------------------------------------
*
* lambrtConfConicToLatLon --
*
* This is a value for a projToLatLonProc member of a GeoProj structure.
* It is usually called by the ProjToLatLon procedure.
*
* Results:
* A geopoint corresponding to a map-point.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
lambrtConfConicToLatLon(mapPt, projPtr)
MapPt mapPt;
struct GeoProj *projPtr;
{
/*
* Formulas and notation from Snyder, p. 105
*/
double lat, lon; /* Latitude and longitude in radians */
GeoProjLambertConfConicParams *params = projPtr->params;
double RF = params->RF;
double n = params->n;
double x = n > 0.0 ? mapPt.abs : -mapPt.abs;
double dy = n > 0.0 ? params->rho0 - mapPt.ord : mapPt.ord - params->rho0;
double rho = n > 0.0
? sqrt(mapPt.abs * mapPt.abs + dy * dy)
: -sqrt(mapPt.abs * mapPt.abs + dy * dy);
double theta = atan2(x, dy);
lat = 2.0 * atan(pow(RF/rho, 1.0/n)) - M_PI_2;
lon = AngleToRad(params->refLon) + theta / n;
return DomainLonPt(GeoPtFmRad(lat, lon), params->refLon);
}
/*
*----------------------------------------------------------------------
*
* SetLambertEqArea --
*
* This procedure sets a projection to Lambert Equal Area.
*
* Results:
* None.
*
* Side effects:
* The members of a GeoProj structure are modified.
*
*----------------------------------------------------------------------
*/
void
SetLambertEqArea(projPtr, refPt)
struct GeoProj *projPtr; /* Projection to modify */
GeoPt refPt; /* Reference point */
{
GeoProjRefPtParams *params;
double lat, lon;
Angle d90 = AngleFmDeg(90.0);
params = (GeoProjRefPtParams *)CKALLOC(sizeof(GeoProjRefPtParams));
projPtr->type = LambertEqArea;
params->refPt = GwchLonPt(refPt);
GeoPtGetRad(params->refPt, &lat, &lon);
params->cosRLat = cos(lat);
params->sinRLat = sin(lat);
if (projPtr->params) {
CKFREE(projPtr->params);
}
projPtr->params = params;
projPtr->projToLatLonProc = lambrtEqAreaToLatLon;
projPtr->latLonToProjProc = latLonToLambrtEqArea;
if (AngleCmp(refPt.lat, d90) == 0) {
sprintf(projPtr->descriptor, "LambertEqArea {90.0 0.0}");
} else if (AngleCmp(refPt.lat, -d90) == 0) {
sprintf(projPtr->descriptor, "LambertEqArea {-90.0 0.0}");
} else {
sprintf(projPtr->descriptor, "LambertEqArea {%9.3f %-9.3f}",
AngleToDeg(params->refPt.lat), AngleToDeg(params->refPt.lon));
}
projPtr->infoProc = refPtProjInfo;
}
/*
*----------------------------------------------------------------------
*
* latLonToLambrtEqArea --
*
* This is a value for a latLonToProjProc member of a GeoProj structure.
* It is usually called by the LatLonToProj procedure.
*
* Results:
* A map-point corresponding to a geopoint.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static MapPt
latLonToLambrtEqArea(geoPt, projPtr)
GeoPt geoPt;
struct GeoProj *projPtr;
{
MapPt mapPt;
GeoProjRefPtParams *params = projPtr->params;
double lat, lon;
double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
double dlon;
double k;
Angle dist;
double cosLat, sinLat, cosDLon;
Angle d90 = AngleFmDeg(90.0);
GeoPtGetRad(geoPt, &lat, &lon);
cosLat = cos(lat);
sinLat = sin(lat);
dlon = AngleToRad(geoPt.lon - params->refPt.lon);
cosDLon = cos(dlon);
dist = GeoDistance(params->refPt, geoPt);
if (AngleCmp(dist, d90) == 1) {
/*
* Point is further than 90 deg from reference point.
*/
return MapPtNowhere();
}
k = 1.0 + sinRLat * sinLat + cosRLat * cosLat * cosDLon;
k = sqrt(2.0 / k);
mapPt.abs = REarth() * k * cosLat * sin(dlon);
mapPt.ord = REarth() * k * (cosRLat * sinLat - sinRLat * cosLat * cosDLon);
return mapPt;
}
/*
*----------------------------------------------------------------------
*
* lambrtEqAreaToLatLon --
*
* This is a value for a projToLatLonProc member of a GeoProj structure.
* It is usually called by the ProjToLatLon procedure.
*
* Results:
* A geopoint corresponding to a map-point.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
lambrtEqAreaToLatLon(mapPt, projPtr)
MapPt mapPt;
struct GeoProj *projPtr;
{
GeoProjRefPtParams *params = projPtr->params;
double rho, c, cosC, sinC, ord;
double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
double lat, lon;
rho = sqrt(mapPt.abs * mapPt.abs + mapPt.ord * mapPt.ord);
if (rho == 0.0) {
return params->refPt;
}
if (rho > 2.0 * REarth()) {
return GeoPtNowhere();
}
c = 2.0 * asin(rho / (2.0 * REarth()));
cosC = cos(c);
sinC = sin(c);
ord = cosC * sinRLat + (mapPt.ord * sinC * cosRLat / rho);
if (ord > 1.0) {
return GeoPtNowhere();
}
lat = asin(ord);
lon = AngleToRad(params->refPt.lon)
+ atan2(mapPt.abs * sinC,
(rho * cosRLat * cosC - mapPt.ord * sinRLat * sinC));
return DomainLonPt(GeoPtFmRad(lat, lon), params->refPt.lon);
}
/*
*----------------------------------------------------------------------
*
* SetStereographic --
*
* This procedure sets a projection to Stereographic.
*
* Results:
* None.
*
* Side effects:
* The members of a GeoProj structure are modified.
*
*----------------------------------------------------------------------
*/
void
SetStereographic(projPtr, refPt)
struct GeoProj *projPtr; /* Projection to modify */
GeoPt refPt; /* Reference point */
{
GeoProjRefPtParams *params;
double lat, lon;
Angle d90 = AngleFmDeg(90.0);
params = (GeoProjRefPtParams *)CKALLOC(sizeof(GeoProjRefPtParams));
projPtr->type = Stereographic;
params->refPt = GwchLonPt(refPt);
GeoPtGetRad(params->refPt, &lat, &lon);
params->cosRLat = cos(lat);
params->sinRLat = sin(lat);
if (projPtr->params) {
CKFREE(projPtr->params);
}
projPtr->params = params;
projPtr->latLonToProjProc = latLonToStereographic;
projPtr->projToLatLonProc = stereographicToLatLon;
if (LatCmp(refPt.lat, d90) == Eq) {
sprintf(projPtr->descriptor, "Stereographic {90.0 0.0}");
} else if (LatCmp(refPt.lat, -d90) == Eq) {
sprintf(projPtr->descriptor, "Stereographic {-90.0 0.0}");
} else {
sprintf(projPtr->descriptor, "Stereographic {%9.3f %-9.3f}",
AngleToDeg(params->refPt.lat), AngleToDeg(params->refPt.lon));
}
projPtr->infoProc = refPtProjInfo;
}
/*
*----------------------------------------------------------------------
*
* latLonToStereographic --
*
* This is a value for a latLonToProjProc member of a GeoProj structure.
* It is usually called by the LatLonToProj procedure.
*
* Results:
* A map-point corresponding to a geopoint.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static MapPt
latLonToStereographic(geoPt, projPtr)
GeoPt geoPt;
struct GeoProj *projPtr;
{
MapPt mapPt;
GeoProjRefPtParams *params = projPtr->params;
GeoPt refPt = params->refPt;
double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
double lat, lon;
double dlon, cosDLon;
double k;
Angle dist;
double cosLat, sinLat;
Angle d90 = AngleFmDeg(90.0);
/*
* Go to 90 degrees from reference point (i.e. follow convention and
* treat as hemisphere projection)
*/
dist = GeoDistance(refPt, geoPt);
if (AngleCmp(dist, d90) == 1) {
return MapPtNowhere();
}
GeoPtGetRad(geoPt, &lat, &lon);
cosLat = cos(lat);
sinLat = sin(lat);
dlon = AngleToRad(geoPt.lon - params->refPt.lon);
cosDLon = cos(dlon);
k = 2.0 / (1.0 + sinRLat * sinLat + cosRLat * cosLat * cosDLon);
mapPt.abs = REarth() * k * cosLat * sin(dlon);
mapPt.ord = REarth() * k * (cosRLat * sinLat - sinRLat * cosLat * cosDLon);
return mapPt;
}
/*
*----------------------------------------------------------------------
*
* stereographicToLatLon --
*
* This is a value for a projToLatLonProc member of a GeoProj structure.
* It is usually called by the ProjToLatLon procedure.
*
* Results:
* A geopoint corresponding to a map-point.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
stereographicToLatLon(mapPt, projPtr)
MapPt mapPt;
struct GeoProj *projPtr;
{
GeoProjRefPtParams *params = projPtr->params;
double rho, c, cosC, sinC, ord;
double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
double lat, lon;
rho = sqrt(mapPt.abs * mapPt.abs + mapPt.ord * mapPt.ord);
if (rho == 0.0) {
return params->refPt;
}
c = 2.0 * atan2(rho, 2.0 * REarth());
cosC = cos(c);
sinC = sin(c);
ord = cosC * sinRLat + (mapPt.ord * sinC * cosRLat / rho);
if (ord > 1.0) {
return GeoPtNowhere();
}
lat = asin(ord);
lon = AngleToRad(params->refPt.lon) + atan2(mapPt.abs * sinC,
(rho * cosRLat * cosC - mapPt.ord * sinRLat * sinC));
return DomainLonPt(GeoPtFmRad(lat, lon), params->refPt.lon);
}
/*
*----------------------------------------------------------------------
*
* SetOrthographic --
*
* This procedure sets a projection to Orthographic.
*
* Results:
* None.
*
* Side effects:
* The members of a GeoProj structure are modified.
*
*----------------------------------------------------------------------
*/
void
SetOrthographic(projPtr, refPt)
struct GeoProj *projPtr; /* Projection to modify */
GeoPt refPt; /* Reference point */
{
double lat, lon;
GeoProjRefPtParams *params;
params = (GeoProjRefPtParams *)CKALLOC(sizeof(GeoProjRefPtParams));
projPtr->type = Orthographic;
params->refPt = GwchLonPt(refPt);
GeoPtGetRad(params->refPt, &lat, &lon);
params->cosRLat = cos(lat);
params->sinRLat = sin(lat);
if (projPtr->params) {
CKFREE(projPtr->params);
}
projPtr->params = params;
projPtr->latLonToProjProc = latLonToOrthographic;
projPtr->projToLatLonProc = orthographicToLatLon;
sprintf(projPtr->descriptor, "Orthographic {%9.3f %-9.3f}",
AngleToDeg(params->refPt.lat), AngleToDeg(params->refPt.lon));
projPtr->infoProc = refPtProjInfo;
}
/*
*----------------------------------------------------------------------
*
* latLonToOrthographic --
*
* This is a value for a latLonToProjProc member of a GeoProj structure.
* It is usually called by the LatLonToProj procedure.
*
* Results:
* A map-point corresponding to a geopoint.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static MapPt
latLonToOrthographic(geoPt, projPtr)
GeoPt geoPt;
struct GeoProj *projPtr;
{
MapPt mapPt;
GeoProjRefPtParams *params = projPtr->params;
GeoPt refPt = params->refPt;
Angle dist;
double lat, lon;
double coslat;
double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
double dlon;
Angle d90 = AngleFmDeg(90.0);
dist = GeoDistance(refPt, geoPt);
if (AngleCmp(dist, d90) == 1) {
/*
* Point further than 90 deg from refPt
*/
return MapPtNowhere();
}
GeoPtGetRad(geoPt, &lat, &lon);
coslat = cos(lat);
dlon = AngleToRad(geoPt.lon - params->refPt.lon);
mapPt.abs = REarth() * coslat * sin(dlon);
mapPt.ord = REarth() * (cosRLat * sin(lat) - sinRLat * coslat * cos(dlon));
return mapPt;
}
/*
*----------------------------------------------------------------------
*
* orthographicToLatLon --
*
* This is a value for a projToLatLonProc member of a GeoProj structure.
* It is usually called by the ProjToLatLon procedure.
*
* Results:
* A geopoint corresponding to a map-point.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
orthographicToLatLon(mapPt, projPtr)
MapPt mapPt;
struct GeoProj *projPtr;
{
GeoProjRefPtParams *params = projPtr->params;
double r; /* Polar radius */
double c; /* Angular distance from geoPt to ref point */
double cosC, sinC, ord; /* Computational constants */
double cosRLat = params->cosRLat, sinRLat = params->sinRLat;
double lat, lon;
r = sqrt(mapPt.abs * mapPt.abs + mapPt.ord * mapPt.ord);
if (r == 0.0) {
return params->refPt;
}
if (r / REarth() > 1.0) {
return GeoPtNowhere();
}
c = asin(r / REarth());
cosC = cos(c);
sinC = sin(c);
ord = cosC * sinRLat + (mapPt.ord * sinC * cosRLat / r);
if (ord > 1.0) {
return GeoPtNowhere();
}
lat = asin(ord);
lon = AngleToRad(params->refPt.lon) + atan2(mapPt.abs * sinC,
(r * cosRLat * cosC - mapPt.ord * sinRLat * sinC));
return DomainLonPt(GeoPtFmRad(lat, lon), params->refPt.lon);
}
/*
*----------------------------------------------------------------------
*
* refLonProjInfo --
*
* This is a value for an infoProc member in a GeoProj structure.
* It is usually called by the GeoProjGetInfo procedure
*
* Results:
* Projection information is copied into a GeoProjInfo structure.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoProjInfo
refLonProjInfo(projPtr)
struct GeoProj *projPtr;
{
GeoProjInfo geoProjInfo;
Angle *refLonPtr = projPtr->params;
geoProjInfo.type = projPtr->type;
geoProjInfo.prm.refLon = *refLonPtr;
geoProjInfo.rotation = projPtr->rotation;
return geoProjInfo;
}
/*
*----------------------------------------------------------------------
*
* refPtProjInfo --
*
* This is a value for an infoProc member in a GeoProj structure.
* It is usually called by the GeoProjGetInfo procedure
*
* Results:
* Projection information is copied into a GeoProjInfo structure.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoProjInfo
refPtProjInfo(projPtr)
struct GeoProj *projPtr;
{
GeoProjRefPtParams *params = projPtr->params;
GeoProjInfo geoProjInfo;
geoProjInfo.type = projPtr->type;
geoProjInfo.prm.refPt = params->refPt;
geoProjInfo.rotation = projPtr->rotation;
return geoProjInfo;
}
syntax highlighted by Code2HTML, v. 0.9.1