/* * 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 #include #include #include #include #include #include #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; }