Subject: Re: Geoid Height Calculation
From: Tony Adams
Date: Thu, 05 Dec 1996 16:53:16 -0800
Jim Frame wrote:
>
> I'm interested in being able to approximate geoid heights for areas
> outside the U.S. and its possessions. Can anyone point me to a utility
> similar to GEOID96 (offered by NGS) with this capability?
>
> ---------------------------------------------------------------------
> Jim Frame jhframe@dcn.davis.ca.us (916) 756-8584 756-8201 (FAX)
> Frame Surveying & Mapping 609 A Street Davis, CA 95616
> -----------------------< Davis Community Network >-------------------
Below is some 'C' code for calculating Geoidial height for the WGS-84
datum. This uses a linear interpolation method fo exctrat heigh from a
surveyed set of readings. This works fine for all latitudes between 80'N
and 80'S and should be accurate to 0.5 metres
I hope this is of some use.
/*********************************************************************/
/* Definitions */
#define M_PI 3.1415927
/**********************************************************************/
/* Function prototypes */
float CalcGeodialHeight(float Lat,float Long);
float Weight(float Delta_Lon,float Delta_Lat);
/**********************************************************************/
/* Code table for look up parameters */
char code GeoidTable[19][36]=
{
-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,
-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,
-53,-54,-55,-52,-48,-42,-38,-38,-29,-26,-26,-24,-23,-21,-19,-16,-12, -8,
-4, -1, 1, 4, 4, 6, 5, 4, 2, -6,-15,-24,-33,-40,-48,-50,-53,-52,
-60,-60,-61,-55,-49,-44,-38,-31,-25,-16, -6, 1, 4, 5, 4, 2, 6, 12,
16, 16, 17, 21, 20, 26, 26, 22, 16, 10, -1,-16,-29,-36,-46,-55,-54,-59,
-45,-43,-37,-32,-30,-26,-23,-22,-16,-10, -2, 10, 20, 20, 21, 24, 22, 17,
16, 19, 25, 30, 35, 35, 33, 30, 27, 10, -2,-14,-23,-30,-33,-29,-35,-43,
-15,-18,-18,-16,-17,-15,-10,-10, -8, -2, 6, 14, 13, 3, 3, 10, 20, 27,
25, 26, 34, 39, 45, 45, 38, 39, 28, 13, -1,-15,-22,-22,-18,-15,-14,-10,
21, 6, 1, -7,-12,-12,-12,-10, -7, -1, 8, 23, 15, -2, -6, 6, 21, 24,
18, 26, 31, 33, 39, 41, 30, 24, 13, -2,-20,-32,-33,-27,-14, -2, 5, 20,
46, 22, 5, -2, -8,-13,-10, -7, -4, 1, 9, 32, 16, 4, -8, 4, 12, 15,
22, 27, 34, 29, 14, 15, 15, 7, -9,-25,-37,-39,-23,-14, 15, 33, 34, 45,
51, 27, 10, 0, -9,-11, -5, -2, -3, -1, 9, 35, 20, -5, -6, -5, 0, 13,
17, 23, 21, 8, -9,-10,-11,-20,-40,-47,-45,-25, 5, 23, 45, 58, 57, 63,
36, 22, 11, 6, -1, -8,-10, -8,-11, -9, 1, 32, 4,-18,-13, -9, 4, 14,
12, 13, -2,-14,-25,-32,-38,-60,-75,-63,-26, 0, 35, 52, 68, 76, 64, 52,
22, 16, 17, 13, 1,-12,-23,-20,-14, -3, 14, 10,-15,-27,-18, 3, 12, 20,
18, 12,-13, -9,-28,-49,-62,-89,-102,-63, -9, 33, 58, 73, 74, 63, 50,
32,
13, 12, 11, 2,-11,-28,-38,-29,-10, 3, 1,-11,-41,-42,-16, 3, 17, 33,
22, 23, 2, -3, -7,-36,-59,-90,-95,-63,-24, 12, 53, 60, 58, 46, 36, 26,
5, 10, 7, -7,-23,-39,-47,-34, -9,-10,-20,-45,-48,-32, -9, 17, 25, 31,
31, 26, 15, 6, 1,-29,-44,-61,-67,-59,-36,-11, 21, 39, 49, 39, 22, 10,
-7, -5, -8,-15,-28,-40,-42,-29,-22,-26,-32,-51,-40,-17, 17, 31, 34, 44,
36, 28, 29, 17, 12,-20,-15,-40,-33,-34,-34,-28, 7, 29, 43, 20, 4, -6,
-12,-10,-13,-20,-31,-34,-21,-16,-26,-34,-33,-35,-26, 2, 33, 59, 52, 51,
52, 48, 35, 40, 33, -9,-28,-39,-48,-59,-50,-28, 3, 23, 37, 18, -1,-11,
-8, 8, 8, 1,-11,-19,-16,-18,-22,-35,-40,-26,-12, 24, 45, 63, 62, 59,
47, 48, 42, 28, 12,-10,-19,-33,-43,-42,-43,-29, -2,-17, 23, 22, 6, 2,
2, 9, 17, 10, 13, 1,-14,-30,-39,-46,-42,-21, 6, 29, 49, 65, 60, 57,
47, 41, 21, 18, 14, 7, -3,-22,-29,-32,-32,-26,-15, -2, 13, 17, 19, 6,
2, 2, 1, -1, -3, -7,-14,-24,-27,-25,-19, 3, 24, 37, 47, 60, 61, 58,
51, 43, 29, 20, 12, 5, -2,-10,-14,-12,-10,-14,-12, -6, -2, 3, 6, 4,
3, 1, -2, -3, -3, -3, -1, 3, 1, 5, 9, 11, 19, 27, 31, 34, 33, 34,
33, 34, 28, 23, 17, 13, 9, 4, 4, 1, -2, -2, 0, 2, 3, 2, 1, 1,
13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13,
};
/**********************************************************************
FUNCTION NAME -
FLOAT WEIGHT(FLOAT DELTA_LON,FLOAT DELTA_LAT)
DESCRIPTION -
Computes the value of a weighted polynomial used in the calculation of
the altitude difference between the WGS-84 ellipsoid and the MSL geoid
INPUTS -
Latitude and longitude increments
OUTPUTS -
Weighted function value
**********************************************************************/
float Weight(float Delta_Lon,float Delta_Lat)
{
float Temp,res;
Temp=Delta_Lon*Delta_Lat;
res= (Temp * Temp) * ( 9- ( 6 * (Delta_Lon +Delta_Lat) )
+ (4 * Temp ) );
return(res);
}
/**********************************************************************
FUNCTION NAME -
FLOAT CALCGEODIALHEIGHT(FLOAT LAT,FLOAT LONG)
DESCRIPTION -
Computes the separation between the WGS-84 ellipsoid and the MSL geoid
INPUTS -
Latitude and longitude in radians
OUTPUTS -
Returns the altitude difference between WGS-84 and the MSL geoid
**********************************************************************/
float CalcGeodialHeight(float Lat,float Long)
{
float LatDeg,LongDeg,xx,yy,diff;
float term1,term2,term3,term4;
int Phi1,Phi2,Lambda1,Lambda2;
LatDeg =((18.0/M_PI)*Lat ) + 9.0;
LongDeg =((18.0/M_PI)*Long) + 18.0;
Phi1 =(int)LatDeg;
Lambda1 =(int)LongDeg;
if ( (Phi1 > -1 ) && (Phi1<19) && (Lambda1>-1) && (Lambda1<37) )
{
yy = LatDeg -Phi1;
xx = LongDeg -Lambda1;
if(LatDeg< 0 )yy=0;
Phi2=Phi1+1;
if(Phi1==18)
{
Phi2=18;
yy=0;
}
if(Lambda1==36)
Lambda1=0;
Lambda2=Lambda1+1;
if(LongDeg<0)
{
Lambda1=35;
xx=1+LongDeg;
}
if(Lambda1==35)
Lambda2=0;
term1=Weight(xx ,yy ) * (float)GeoidTable[Phi2][Lambda2];
term2=Weight(1-xx,yy ) * (float)GeoidTable[Phi2][Lambda1];
term3=Weight(1-xx,1-yy) * (float)GeoidTable[Phi1][Lambda1];
term4=Weight(xx ,1-yy) * (float)GeoidTable[Phi1][Lambda2];
diff=term1+term2+term3+term4;
}
return(diff);
}
/**********************************************************************/
