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