/*
* geoLnArrToMap.c --
*
* This file defines functions that convert the geopoints in geolinearrays
* to map-points.
*
* 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: geoLnArrToMap.c,v 1.7 2007/06/29 15:59:53 tkgeomap Exp $
*
********************************************
*
*/
#include <stdio.h>
#include <math.h>
#include <tcl.h>
#include "geography.h"
#include "geoLines.h"
#include "mapLines.h"
#include "geoProj.h"
#include "geoLnArrToMap.h"
/*
* Declarations of local functions.
*/
static GeoPt noAdjust _ANSI_ARGS_((GeoPt geoPt));
static GeoPt offPole _ANSI_ARGS_((GeoPt geoPt));
static GeoPt nToEq _ANSI_ARGS_((GeoPt geoPt));
static GeoPt sToEq _ANSI_ARGS_((GeoPt geoPt));
static MapLnArr geoLnArrToCylMap _ANSI_ARGS_((GeoLnArr geoLnArr,
GeoProj proj));
static MapLnArr geoLnArrToHemisphere _ANSI_ARGS_((GeoLnArr geoLnArr,
GeoProj proj));
/*
*----------------------------------------------------------------------
*
* GeoLnArrToMap --
*
* This procedure converts the geopoints in a geolinearray to map-points.
*
* Results:
* Return value is a new MapLnArr or NULL if something goes wrong.
*
* Side effects:
* Storage for the returned array is dynamically allocated and should
* eventually be freed with a call to MapLnArrDestroy;
*
*----------------------------------------------------------------------
*/
MapLnArr
GeoLnArrToMap(geoLnArr, proj)
GeoLnArr geoLnArr;
GeoProj proj;
{
GeoProjInfo projInfo; /* GeoProj information */
MapLnArr mapLnArr; /* Return value */
if (!proj) {
return NULL;
}
projInfo = GeoProjGetInfo(proj);
switch (projInfo.type) {
case CylEqDist:
case LambertConfConic:
case CylEqArea:
case Mercator:
mapLnArr = geoLnArrToCylMap(geoLnArr, proj);
break;
case Orthographic:
case Stereographic:
case LambertEqArea:
mapLnArr = geoLnArrToHemisphere(geoLnArr, proj);
break;
default:
return NULL;
}
if ( !mapLnArr ) {
return NULL;
}
MapLnArrSet(mapLnArr, geoLnArr, proj);
return mapLnArr;
}
/*
*----------------------------------------------------------------------
*
* noAdjust --
*
* This procedure just returns its input. It is a place holder.
*
* Results:
* The input value.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
noAdjust(geoPt)
GeoPt geoPt;
{
return geoPt;
}
/*
*----------------------------------------------------------------------
*
* offPole --
*
* This procedure moves a point off the North or South pole. This is
* necessary for cylindrical projections which are undefined at the poles.
*
* Results:
* If geoPt is not on a pole, return value is geoPt. If geoPt is on a
* pole, return value is a point near geoPt.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
offPole(geoPt)
GeoPt geoPt;
{
Angle nearAng = AngleFmDeg(89.0);
Angle d90 = AngleFmDeg(90.0);
geoPt.lat = (LatCmp(geoPt.lat, d90) == Eq) ? nearAng
: (LatCmp(geoPt.lat, -d90) == Eq) ? -nearAng : geoPt.lat;
return geoPt;
}
/*
*----------------------------------------------------------------------
*
* sToEq --
*
* This procedure moves a point in the Southern Hemisphere to the equator.
*
* Results:
* A point with the same longitude as the given geopoint, with a latitude
* not in the Southern Hemisphere.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
sToEq(geoPt)
GeoPt geoPt;
{
geoPt.lat = (LatCmp(geoPt.lat, 0) == South) ? 0 : geoPt.lat;
return geoPt;
}
/*
*----------------------------------------------------------------------
*
* nToEq --
*
* This procedure moves a point in the Northern Hemisphere to the equator.
*
* Results:
* A point with the same longitude as the given geopoint, with a latitude
* not in the Northern Hemisphere.
*
* Side effects:
* None.
*
*----------------------------------------------------------------------
*/
static GeoPt
nToEq(geoPt)
GeoPt geoPt;
{
geoPt.lat = (LatCmp(geoPt.lat, 0) == North) ? 0 : geoPt.lat;
return geoPt;
}
/*
*----------------------------------------------------------------------
*
* geoLnArrToCylMap --
*
* This procedure applies a cylindrical projection to the geopoints in a
* geolinearray. Then it rotates the projected points.
*
* In a cylindrical projection, the lat-lon's are projected onto a cylinder
* that wraps around the Earth with a longitude discontinuity 180 degrees
* from the reference longitude. Segments that cross the discontinuity are
* split.
*
* Results:
* Return value is a new mapLnArr or NULL if something goes wrong.
*
* Side effects:
* Storage for the returned array is dynamically allocated and should
* eventually be freed with a call to MapLnArrDestroy;
*
*----------------------------------------------------------------------
*/
static MapLnArr
geoLnArrToCylMap(geoLnArr, proj)
GeoLnArr geoLnArr;
GeoProj proj;
{
MapLnArr mapLnArr = NULL; /* Return value */
unsigned nl, nls, np; /* Loop parameters */
GeoProjInfo projInfo; /* Projection information */
GeoPt (*adjust)(GeoPt) = noAdjust; /* Function to bring a point into
* projection domain, if necessary */
Angle refLon = 0.0; /* Reference longitude */
Angle iWest, iEast; /* West and east end of domain */
Angle lon0, lon1; /* Longitudes of last and current pt */
GeoPt geoPt0, geoPt1; /* Points in the current segment */
GeoPt cross; /* Point on domain boundary */
enum LonSgn side0, side1; /* Which side of rLon geoPt0 and geoPt1
* are on */
enum LonSgn side; /* If line has touched rLon, this stores
* which side it approached from */
GeoLn geoLn = NULL; /* Current line in latlon coordinates */
MapLn mapLn[2] = {NULL, NULL}; /* Current line in map coordinates */
MapLn currLn; /* Line we are adding points to */
int i; /* Counter to toggle mapLn index */
Angle d180 = AngleFmDeg(180.0);
nls = geoLnArr->nLines;
if ( !(mapLnArr = MapLnArrCreate(nls))
|| !(mapLn[0] = MapLnCreate(1000))
|| !(mapLn[1] = MapLnCreate(1000)) ) {
MapLnDestroy(mapLn[0]);
MapLnDestroy(mapLn[1]);
MapLnArrDestroy(mapLnArr);
return NULL;
}
projInfo = GeoProjGetInfo(proj);
switch (projInfo.type) {
case CylEqDist:
refLon = projInfo.prm.refPt.lon;
adjust = noAdjust;
break;
case LambertConfConic:
refLon = projInfo.prm.refPt.lon;
if (LatCmp(projInfo.prm.refPt.lat, 0) == North) {
adjust = sToEq;
} else {
adjust = nToEq;
}
break;
case CylEqArea:
refLon = projInfo.prm.refLon;
adjust = noAdjust;
break;
case Mercator:
refLon = projInfo.prm.refLon;
adjust = offPole;
break;
case Orthographic:
case Stereographic:
case LambertEqArea:
return NULL;
}
iWest = refLon - d180;
iEast = refLon + d180;
if (refLon > 0) {
/*
* If reference longitude is in the eastern hemisphere, split lines
* that cross the western boundary.
*/
for (nl = 0; nl < nls; nl++) {
geoLn = GeoLnArrGetLine(geoLnArr, nl);
currLn = mapLn[(i = 0)];
side = PrMd;
/*
* Process the first segment.
*/
geoPt0 = adjust(GeoLnGetPt(geoLn, 0));
geoPt1 = adjust(GeoLnGetPt(geoLn, 1));
lon0 = DomainLon(geoPt0.lon, iWest);
lon1 = DomainLon(geoPt1.lon, iWest);
side0 = LonCmp(geoPt0.lon, iWest);
side1 = LonCmp(geoPt1.lon, iWest);
if (side0 == PrMd && side1 == West) {
/*
* Line is departing from western boundary to the west.
* Draw the segment at the east end.
*/
geoPt0.lon = iEast;
}
if (side0 != PrMd && side1 == PrMd) {
/*
* If line has just intersected west, make a note
* of which side it is approaching from.
*/
side = side0;
if (side0 == West) {
/*
* If line has arrived at west side from beyond, draw
* second point at east end.
*/
geoPt1.lon = iEast;
}
}
if (LonBtwn(iWest, lon0, lon1)) {
/*
* Segment is crossing western boundary.
*/
cross.lat = geoPt0.lat + (double)(geoPt1.lat - geoPt0.lat)
/ (lon1 - lon0) * (iWest - lon0);
if (LonCmp(lon0, iWest) == West) {
/*
* Crossing west to east.
*/
cross.lon = iEast;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
currLn = mapLn[++i % 2];
cross.lon = iWest;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
} else {
/*
* Crossing east to west.
*/
cross.lon = iWest;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
currLn = mapLn[++i % 2];
cross.lon = iEast;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
}
} else {
/*
* No crossing => no special treatment.
*/
MapLnAddPt(LatLonToProj(geoPt0, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
}
geoPt0 = geoPt1;
for (np = 2; np < geoLn->nPts; np++) {
/*
* Process the rest of the segments.
*/
geoPt1 = adjust(GeoLnGetPt(geoLn, np));
lon0 = DomainLon(geoPt0.lon, iWest);
lon1 = DomainLon(geoPt1.lon, iWest);
side0 = LonCmp(geoPt0.lon, iWest);
side1 = LonCmp(geoPt1.lon, iWest);
if (side0 != PrMd && side1 == PrMd) {
/*
* If line has just arrived at west, make a note
* of which side it is approaching from.
*/
side = side0;
}
if (side != PrMd && side0 == PrMd
&& side1 != PrMd && side1 != side) {
/*
* Line has intersected west and
* then left on the opposite side.
*/
currLn = mapLn[++i % 2];
geoPt0.lon = (side == East) ? iEast : iWest;
MapLnAddPt(LatLonToProj(geoPt0, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
side = PrMd;
} else if (side != PrMd && side0 == PrMd
&& side1 != PrMd && side1 == side) {
/*
* Line has traveled along rLon and
* then left on the same side.
*/
if (side == West) {
geoPt1.lon = DomainLon(geoPt1.lon, iEast);
}
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
side = PrMd;
} else if (LonBtwn(iWest, lon0, lon1)) {
/*
* Segment is crossing western boundary.
*/
cross.lat = geoPt0.lat
+ (double)(geoPt1.lat - geoPt0.lat) / (lon1 - lon0)
* (iWest - lon0);
if (LonCmp(lon0, iWest) == West) {
/*
* Crossing west to east.
*/
cross.lon = iEast;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
currLn = mapLn[++i % 2];
cross.lon = iWest;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
} else {
/*
* Crossing east to west.
*/
cross.lon = iWest;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
currLn = mapLn[++i % 2];
cross.lon = iEast;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
}
} else {
/*
* Segment is not crossing a boundary. Just add the point.
*/
if (side == West) {
geoPt1.lon = DomainLon(geoPt1.lon, iEast);
}
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
}
geoPt0 = geoPt1;
}
/*
* If projected rotated line has any points, add it to the
* output array.
*/
if (mapLn[0]->nPts > 1) {
MapLnArrAddLine(mapLn[0], mapLnArr);
}
if (mapLn[1]->nPts > 1) {
MapLnArrAddLine(mapLn[1], mapLnArr);
}
MapLnClear(mapLn[0]);
MapLnClear(mapLn[1]);
}
} else {
/*
* Reference longitude <= 0, i.e. reference longitude is in the
* western hemisphere. Split lines that cross the eastern boundary.
*/
for (nl = 0; nl < nls; nl++) {
geoLn = GeoLnArrGetLine(geoLnArr, nl);
currLn = mapLn[(i = 0)];
side = PrMd;
/*
* Process the first segment.
*/
geoPt0 = adjust(GeoLnGetPt(geoLn, 0));
geoPt1 = adjust(GeoLnGetPt(geoLn, 1));
lon0 = DomainLon(geoPt0.lon, iEast);
lon1 = DomainLon(geoPt1.lon, iEast);
side0 = LonCmp(geoPt0.lon, iEast);
side1 = LonCmp(geoPt1.lon, iEast);
if (side0 == PrMd && side1 == East) {
/*
* Line is departing from eastern boundary to the east.
* Draw the segment at the west end.
*/
geoPt0.lon = iWest;
}
if (side0 != PrMd && side1 == PrMd) {
/*
* If line has just intersected east, make a note
* of which side it is approaching from.
*/
side = side0;
if (side0 == East) {
/*
* If line has arrived at east side from beyond, draw
* second point at west end.
*/
geoPt1.lon = iWest;
}
}
if (LonBtwn(iEast, lon0, lon1)) {
/*
* Segment is crossing eastern boundary.
*/
cross.lat = geoPt0.lat + (double)(geoPt1.lat - geoPt0.lat)
/ (lon1 - lon0) * (iEast - lon0);
if (LonCmp(lon0, iEast) == West) {
/*
* Crossing west to east
*/
cross.lon = iEast;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
cross.lon = iWest;
currLn = mapLn[++i % 2];
MapLnAddPt(LatLonToProj(cross, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
} else {
/*
* Crossing east to west
*/
cross.lon = iWest;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
cross.lon = iEast;
currLn = mapLn[++i % 2];
MapLnAddPt(LatLonToProj(cross, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
}
} else {
/*
* No crossing => no special treatment.
*/
MapLnAddPt(LatLonToProj(geoPt0, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
}
geoPt0 = geoPt1;
for (np = 2; np < geoLn->nPts; np++) {
/*
* Process the rest of the segments.
*/
geoPt1 = adjust(GeoLnGetPt(geoLn, np));
lon0 = DomainLon(geoPt0.lon, iEast);
lon1 = DomainLon(geoPt1.lon, iEast);
side0 = LonCmp(geoPt0.lon, iEast);
side1 = LonCmp(geoPt1.lon, iEast);
if (side0 != PrMd && side1 == PrMd) {
/*
* If line has just intersected east, make a note
* of which side it is approaching from.
*/
side = side0;
}
if (side != PrMd && side0 == PrMd
&& side1 != PrMd && side1 != side) {
/*
* Line has intersected east and
* then left on the opposite side.
*/
currLn = mapLn[++i % 2];
geoPt0.lon = (side == West) ? iWest : iEast;
MapLnAddPt(LatLonToProj(geoPt0, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
side = PrMd;
} else if (side != PrMd && side0 == PrMd
&& side1 != PrMd && side1 == side) {
/*
* Line has intersected east and
* then left on the same side.
*/
if (side == East) {
geoPt1.lon = DomainLon(geoPt1.lon, iWest);
}
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
side = PrMd;
} else if (LonBtwn(iEast, lon0, lon1)) {
/*
* Segment is crossing eastern boundary.
*/
cross.lat = geoPt0.lat
+ (double)(geoPt1.lat - geoPt0.lat)
/ (lon1 - lon0) * (iEast - lon0);
if (LonCmp(lon0, iEast) == West) {
/*
* Crossing west to east
*/
cross.lon = iEast;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
cross.lon = iWest;
currLn = mapLn[++i % 2];
MapLnAddPt(LatLonToProj(cross, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
} else {
/*
* Crossing east to west
*/
cross.lon = iWest;
MapLnAddPt(LatLonToProj(cross, proj), currLn);
cross.lon = iEast;
currLn = mapLn[++i % 2];
MapLnAddPt(LatLonToProj(cross, proj), currLn);
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
}
} else {
/*
* Segment is not crossing a boundary. Just add the point.
*/
if (side == East) {
geoPt1.lon = DomainLon(geoPt1.lon, iWest);
}
MapLnAddPt(LatLonToProj(geoPt1, proj), currLn);
}
geoPt0 = geoPt1;
}
/*
* If projected rotated line has any points, add it to the
* output array.
*/
if (mapLn[0]->nPts > 1) {
MapLnArrAddLine(mapLn[0], mapLnArr);
}
if (mapLn[1]->nPts > 1) {
MapLnArrAddLine(mapLn[1], mapLnArr);
}
MapLnClear(mapLn[0]);
MapLnClear(mapLn[1]);
}
}
MapLnDestroy(mapLn[0]);
MapLnDestroy(mapLn[1]);
MapLnArrSetAlloc(mapLnArr, mapLnArr->nLines);
return mapLnArr;
}
/*
*----------------------------------------------------------------------
*
* geoLnArrToHemisphere --
*
* This procedure applies a hemispheric projection to the geopoints in
* a geolinearray. Then it rotates the projected points.
*
* In a hemispheric projection, points are defined for the hemisphere
* centered on the projection's reference point. Segments that cross the
* great circle 90 degrees from the reference point are truncated.
*
* Results:
* Return value is a new mapLnArr or NULL if something goes wrong.
*
* Side effects:
* Storage for the returned array is dynamically allocated and should
* eventually be freed with a call to MapLnArrDestroy;
*
*----------------------------------------------------------------------
*/
static MapLnArr
geoLnArrToHemisphere(geoLnArr, proj)
GeoLnArr geoLnArr;
GeoProj proj;
{
MapLnArr mapLnArr = NULL; /* Return value */
unsigned nl, nls, np, nps; /* Loop parameters */
GeoProjInfo projInfo; /* GeoProj info */
GeoPt refPt; /* Reference point */
GeoPt ccl0, ccl1; /* Two points on domain boundary */
GeoPt geoPt0, geoPt1; /* Points in the current segment */
GeoPt bndPt, exitPt, entr0; /* Points on domain boundary */
int in0, in1 = 0; /* True if point in domain */
int in00; /* True if 1st point in line in domain */
Angle az0, az1, az; /* Arc angle and increment */
Angle dAz = AngleFmDeg(2.0);/* Angle increment in loop */
MapLn mapLn = NULL; /* Line being made right now */
MapPt mapPt; /* Point to add to mapLn */
Angle d90 = AngleFmDeg(90.0);
/*
* Get reference point from projection
*/
projInfo = GeoProjGetInfo(proj);
refPt = projInfo.prm.refPt;
ccl0 = GeoStep(refPt, 0, d90);
ccl1 = GeoStep(refPt, d90, d90);
nls = geoLnArr->nLines;
if (!(mapLnArr = MapLnArrCreate(nls)) || !(mapLn = MapLnCreate(1000))) {
MapLnDestroy(mapLn);
MapLnArrDestroy(mapLnArr);
return NULL;
}
for (nl = 0; nl < nls; nl++) {
/*
* Convert a line from the array to a rotated projection coordinates
*/
GeoLn geoPts = GeoLnArrGetLine(geoLnArr, nl);
nps = geoPts->nPts;
entr0 = exitPt = GeoPtNowhere();
/*
* Handle first point in line
*/
geoPt0 = GeoLnGetPt(geoPts, 0);
in0 = in00 = (AngleCmp(GeoDistance(refPt, geoPt0), d90) == -1);
if (in0) {
mapPt = LatLonToProj(geoPt0, proj);
MapLnAddPt(mapPt, mapLn);
}
/*
* Loop through the rest of the points
*/
for (np = 1; np < nps; np++) {
geoPt1 = GeoLnGetPt(geoPts, np);
in1 = (AngleCmp(GeoDistance(refPt, geoPt1), d90) == -1);
if (in0 && in1) {
/*
* Segment is in projection domain
*/
mapPt = LatLonToProj(geoPt1, proj);
MapLnAddPt(mapPt, mapLn);
} else if (in0 && !in1) {
/*
* Segment is leaving projection domain.
*/
bndPt = GCircleX(ccl0, ccl1, geoPt0, geoPt1);
az = Azimuth(bndPt, refPt);
bndPt = GeoStep(bndPt, az, 100);
exitPt = bndPt;
mapPt = LatLonToProj(bndPt, proj);
MapLnAddPt(mapPt, mapLn);
} else if (!in0 && in1) {
/*
* Segment is entering projection domain.
*/
bndPt = GCircleX(ccl0, ccl1, geoPt0, geoPt1);
az = Azimuth(bndPt, refPt);
bndPt = GeoStep(bndPt, az, 100);
/*
* If entering for first time, note location so we can close
* the polygon when done
*/
if ( !in00 && GeoPtIsNowhere(entr0) ) {
entr0 = bndPt;
}
/*
* If line has left domain and is coming back in,
* draw arc from exit point to entry point (bndPt)
*/
if (GeoPtIsSomewhere(exitPt)) {
az0 = Azimuth(refPt, exitPt);
az1 = DomainLon(Azimuth(refPt, bndPt), az0);
if (az1 > az0) {
for (az = az0 + dAz; az < az1; az += dAz) {
GeoPt pt = GeoStep(refPt, az, d90 - 100);
mapPt = LatLonToProj(pt, proj);
MapLnAddPt(mapPt, mapLn);
}
} else {
for (az = az0 - dAz; az > az1; az -= dAz) {
GeoPt pt = GeoStep(refPt, az, d90 - 100);
mapPt = LatLonToProj(pt, proj);
MapLnAddPt(mapPt, mapLn);
}
}
}
/*
* Add entry point and current point
*/
mapPt = LatLonToProj(bndPt, proj);
MapLnAddPt(mapPt, mapLn);
mapPt = LatLonToProj(geoPt1, proj);
MapLnAddPt(mapPt, mapLn);
}
geoPt0 = geoPt1;
in0 = in1;
}
/*
* If line started and finished outside domain
* draw arc from last exit to first entry to close the polygon
*/
if ( !in00 && !in1 && GeoPtIsSomewhere(entr0) ) {
az0 = Azimuth(refPt, exitPt);
az1 = DomainLon(Azimuth(refPt, entr0), az0);
if (az1 > az0) {
for (az = az0 + dAz; az < az1; az += dAz) {
GeoPt pt = GeoStep(refPt, az, d90 - 100);
mapPt = LatLonToProj(pt, proj);
MapLnAddPt(mapPt, mapLn);
}
} else {
for (az = az0 - dAz; az > az1; az -= dAz) {
GeoPt pt = GeoStep(refPt, az, d90 - 100);
mapPt = LatLonToProj(pt, proj);
MapLnAddPt(mapPt, mapLn);
}
}
}
if (mapLn->nPts > 1) {
MapLnArrAddLine(mapLn, mapLnArr);
}
MapLnClear(mapLn);
}
MapLnDestroy(mapLn);
return mapLnArr;
}
syntax highlighted by Code2HTML, v. 0.9.1