Cartography-Projection-GCTP
view release on metacpan or search on metacpan
gctpc/alconinv.c view on Meta::CPAN
fpxyi = yp * cr + xp * ci + ci;
den = fpxyr * fpxyr + fpxyi * fpxyi;
dxp = -(fxyr * fpxyr + fxyi * fpxyi) / den;
dyp = -(fxyi * fpxyr - fxyr * fpxyi) / den;
xp = xp + dxp;
yp = yp + dyp;
ds = fabs(dxp) + fabs(dyp);
nn++;
if (nn > 20)
{
p_error("Too many iterations in inverse","alcon-inv");
return(235);
}
}
while (ds > EPSLN);
/* convert Oblique Stereographic coordinates to LAT/LONG
------------------------------------------------------*/
rh = sqrt(xp * xp + yp * yp);
z = 2.0 * atan(rh / 2.0);
sincos(z,&sinz,&cosz);
gctpc/alconinv.c view on Meta::CPAN
phi = chi;
do
{
esphi = e * sin(phi);
dphi = 2.0 * atan(tan((HALF_PI + chi) / 2.0) *
pow(((1.0 + esphi) / (1.0 - esphi)),(e / 2.0))) - HALF_PI - phi;
phi += dphi;
nn++;
if (nn > 20)
{
p_error("Too many iterations in inverse","alcon-inv");
return(236);
}
}
while(fabs(dphi) > EPSLN);
*lat = phi;
*lon = adjust_lon (lon_center + atan2((xp * sinz), (rh * cos_p26 * cosz - yp *
sin_p26 * sinz)));
gctpc/cproj.c view on Meta::CPAN
{
dphi = (ml + e1 * sin(2.0 * phi) - e2 * sin(4.0 * phi) + e3 * sin(6.0 * phi))
/ e0 - phi;
phi += dphi;
if (fabs(dphi) <= .0000000001)
{
*flag = 0;
return(phi);
}
}
p_error("Latitude failed to converge after 15 iterations","PHI3Z-CONV");
*flag = 3;
return(3);
}
/* Function to compute, phi4, the latitude for the inverse of the
Polyconic projection.
------------------------------------------------------------*/
double phi4z (eccent,e0,e1,e2,e3,a,b,c,phi)
double eccent; /* Spheroid eccentricity squared */
gctpc/goodfor.c view on Meta::CPAN
PURPOSE: Transforms input longitude and latitude to Easting and
Northing for the Goode's Homolosine projection. The
longitude and latitude must be in radians. The Easting
and Northing values will be returned in meters.
PROGRAMMER DATE
---------- ----
D. Steinwand, EROS May, 1991; Updated Sept, 1992; Updated Feb 1993
S. Nelson, EDC Jun, 1993 Make changes in precision and number
of iterations.
ALGORITHM REFERENCES
1. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
U.S. Geological Survey Professional Paper 1453 , United State Government
Printing Office, Washington D.C., 1989.
2. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
State Government Printing Office, Washington D.C., 1987.
gctpc/molwfor.c view on Meta::CPAN
PURPOSE: Transforms input longitude and latitude to Easting and
Northing for the MOllweide projection. The
longitude and latitude must be in radians. The Easting
and Northing values will be returned in meters.
PROGRAMMER DATE
---------- ----
D. Steinwand, EROS May, 1991; Updated Sept, 1992; Updated Feb, 1993
S. Nelson, EDC Jun, 2993; Made corrections in precision and
number of iterations.
ALGORITHM REFERENCES
1. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
U.S. Geological Survey Professional Paper 1453 , United State Government
Printing Office, Washington D.C., 1989.
2. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
State Government Printing Office, Washington D.C., 1987.
gctpc/robinv.c view on Meta::CPAN
NAME ROBINSON
PURPOSE: Transforms input Easting and Northing to longitude and
latitude for the Robinson projection. The
Easting and Northing must be in meters. The longitude
and latitude values will be returned in radians.
PROGRAMMER DATE REASON
---------- ---- ------
T. Mittan March, 1993 Converted from FORTRAN to C.
S. Nelson Nov, 1993 Changed number of iterations from 20
to 75. This seemed to give a valid
Latitude with less fatal errors.
This function was adapted from the Robinson projection code (FORTRAN)
in the General Cartographic Transformation Package software which is
available from the U.S. Geological Survey National Mapping Division.
ALGORITHM REFERENCES
1. "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,
gctpc/robinv.c view on Meta::CPAN
* p2 * (pr[ip1 + 3] - 2.0 * pr[ip1 + 2] + pr[ip1 + 1])/2.0)
* PI / 2.0;
else
y1 = -R * (pr[ip1 +2] + p2 *(pr[ip1 + 3] - pr[ip1 +1]) / 2.0 + p2
* p2 * (pr[ip1 + 3] - 2.0 * pr[ip1 + 2] + pr[ip1 + 1])/2.0)
* PI / 2.0;
phid += -180.0 * (y1 - y) / PI / R;
i++;
if (i > 75)
{
p_error("Too many iterations in inverse","robinv-conv");
return(234);
}
}
while (fabs(y1 - y) > .00001);
break;
}
else
{
ip1 -= 1;
if (ip1 < 0)
{
p_error("Too many iterations in inverse","robinv-conv");
return(234);
}
}
}
*lat = phid * .01745329252;
/* calculate LONG. using final LAT. with transposed forward Stirling's
interpolation formula.
---------------------------------------------------------------------*/
*lon = lon_center + x / R / (xlr[ip1 + 2] + p2 * (xlr[ip1 + 3] - xlr[ip1 + 1])
gctpc/somfor.c view on Meta::CPAN
/* Adjust for confusion at beginning and end of landsat orbits
-----------------------------------------------------------*/
L250: rlm=PI*LANDSAT_RATIO;
rlm2=rlm+2.0*PI;
n++;
if(n>=3) goto L300;
if(tlam>rlm&&tlam<rlm2) goto L300;
if(tlam<rlm)tlamp=2.50*PI;
if(tlam>=rlm2) tlamp=HALF_PI;
goto L230;
L260: sprintf(errorbuf,"50 iterations without conv\n");
p_error(errorbuf,"som-forward");
return(214);
/* tlam computed - now compute tphi
--------------------------------*/
L300: ds=sin(tlam);
dd=ds*ds;
dp=sin(radlt);
tphi=asin(((1.0-es)*ca*dp-sa*cos(radlt)*sin(xlamt))/sqrt(1.0-es*dp*dp));
gctpc/sominv.c view on Meta::CPAN
sdsq=sd*sd;
s=p21*sa*cos(tlon)*sqrt((1.0+t*sdsq)/((1.0+w*sdsq)*(1.0+q*sdsq)));
blon=(x/a)+(y/a)*s/xj-a2*sin(2.0*tlon)-a4*sin(4.0*tlon)-(s/xj)*(c1*
sin(tlon)+c3*sin(3.0*tlon));
tlon=blon/b;
dif=tlon-sav;
if(fabs(dif)<conv)break;
}
if(inumb>=50)
{
p_error("50 iterations without convergence","som-inverse");
return(214);
}
/* Compute transformed lat.
------------------------*/
st=sin(tlon);
defac=exp(sqrt(1.0+s*s/xj/xj)*(y/a-c1*st-c3*sin(3.0*tlon)));
actan=atan(defac);
tlat=2.0*(actan-(PI/4.0));
gctpc/tminv.c view on Meta::CPAN
double y; /* (I) Y projection coordinate */
double *lon; /* (O) Longitude */
double *lat; /* (O) Latitude */
{
double con,phi; /* temporary angles */
double delta_phi; /* difference between longitudes */
long i; /* counter variable */
double sin_phi, cos_phi, tan_phi; /* sin cos and tangent values */
double c, cs, t, ts, n, r, d, ds; /* temporary variables */
double f, h, g, temp; /* temporary variables */
long max_iter = 6; /* maximun number of iterations */
/* fortran code for spherical form
--------------------------------*/
if (ind != 0)
{
f = exp(x/(r_major * scale_factor));
g = .5 * (f - 1/f);
temp = lat_origin + y/(r_major * scale_factor);
h = cos(temp);
con = sqrt((1.0 - h * h)/(1.0 + g * g));
gctpc/utminv.c view on Meta::CPAN
double y; /* (I) Y projection coordinate */
double *lon; /* (O) Longitude */
double *lat; /* (O) Latitude */
{
double con,phi; /* temporary angles */
double delta_phi; /* difference between longitudes */
long i; /* counter variable */
double sin_phi, cos_phi, tan_phi; /* sin cos and tangent values */
double c, cs, t, ts, n, r, d, ds; /* temporary variables */
double f, h, g, temp; /* temporary variables */
long max_iter = 6; /* maximun number of iterations */
/* fortran code for spherical form
--------------------------------*/
if (ind != 0)
{
f = exp(x/(r_major * scale_factor));
g = .5 * (f - 1/f);
temp = lat_origin + y/(r_major * scale_factor);
h = cos(temp);
con = sqrt((1.0 - h * h)/(1.0 + g * g));
( run in 1.190 second using v1.01-cache-2.11-cpan-96521ef73a4 )