android经纬度转wev墨卡托,[转载]UTM与经纬度之间的转换

I get enough inquiries on this subject that I decided to create a

page for it.

Caution! Unlike latitude and longitude, there is no physical frame

of reference for UTM grids. Latitude is determined by the earth's

polar axis. Longitude is determined by the earth's rotation. If you

can see the stars and have a sextant and a good clock set to

Greenwich time, you can find your latitude and longitude. But there

is no way to determine your UTM coordinates except by

calculation.

UTM grids, on the other hand, are created by laying a square grid

on the earth. This means that different maps will have different

grids depending on the datum used (model of the shape of the

earth). I saw US military maps of Germany shift their UTM grids by

about 300 meters when a more modern datum was used for the maps.

Also, old World War II era maps of Europe apparently used a single

grid for all of Europe and grids in some areas are wildly tilted

with respect to latitude and longitude.

The two basic references for converting UTM and geographic

coordinates are U.S. Geological Survey Professional Paper 1395 and

U. S. Army Technical Manual TM 5-241-8 (complete citations below).

Each has advantages and disadvantages.

For converting latitude and longitude to UTM, the Army publication

is better. Its notation is more consistent and the formulas more

clearly laid out, making code easier to debug. In defense of the

USGS, their notation is constrained by space, and the need to be

consistent with cartographic literature and descriptions of several

dozen other map projections in the book.

For converting UTM to latitude and longitude, however, the Army

publication has formulas that involve latitude (the quantity to be

found) and which require reverse interpolation from tables. Here

the USGS publication is the only game in town.

Some extremely tiny terms that will not seriously affect

meter-scale accuracy have been omitted.

Converting Between Decimal Degrees, Degrees, Minutes and Seconds,

and Radians

(dd + mm/60 +ss/3600) to Decimal degrees (dd.ff)

dd = whole degrees, mm = minutes, ss = seconds

dd.ff = dd + mm/60 + ss/3600

Example: 30 degrees 15 minutes 22 seconds = 30 + 15/60 + 22/3600 =

30.2561

Decimal degrees (dd.ff) to (dd + mm/60 +ss/3600)

For the reverse conversion, we want to convert dd.ff to dd mm ss.

Here ff = the fractional part of a decimal degree.

mm = 60*ff

ss = 60*(fractional part of mm)

Use only the whole number part of mm in the final result.

30.2561 degrees = 30 degrees

.2561*60 = 15.366 minutes

.366 minutes = 22 seconds, so the final result is 30 degrees 15

minutes 22 seconds

Decimal degrees (dd.ff) to Radians

Radians = (dd.ff)*pi/180

Radians to Decimal degrees (dd.ff)

(dd.ff) = Radians*180/pi

Degrees, Minutes and Seconds to Distance

A degree of longitude at the equator is 111.2 kilometers. A minute

is 1853 meters. A second is 30.9 meters. For other latitudes

multiply by cos(lat). Distances for degrees, minutes and

secondsin latitude are very similar and

differ very slightly with latitude. (Before satellites, observing

those differences was a principal method for determining the exact

shape of the earth.)

Converting Latitude and Longitude to UTM

Okay, take a deep breath. This will

get very complicated, but

the math, although tedious, is only algebra and

trigonometry. It would sure be nice if someone

wrote aspreadsheet to

do this.

P = point under consideration

F = foot of perpendicular from P to the central meridian. The

latitude of F is called thefootprint latitude.

O = origin (on equator)

OZ = central meridian

LP = parallel of latitude of P

ZP = meridian of P

OL = k0S = meridional arc from equator

LF = ordinate of curvature

OF = N = grid northing

FP = E = grid distance from central meridian

GN = grid north

C = convergence of meridians = angle between true and grid

north

Another thing you need to know is the datum being used:

Datum

Equatorial Radius, meters (a)

Polar Radius, meters (b)

Flattening (a-b)/a

Use

NAD83/WGS84

6,378,137

6,356,752.3142

1/298.257223563

Global

GRS 80

6,378,137

6,356,752.3141

1/298.257222101

US

WGS72

6,378,135

6,356,750.5

1/298.26

NASA, DOD

Australian 1965

6,378,160

6,356,774.7

1/298.25

Australia

Krasovsky 1940

6,378,245

6,356,863.0

1/298.3

Soviet Union

International (1924) -Hayford (1909)

6,378,388

6,356,911.9

1/297

Global except as listed

Clake 1880

6,378,249.1

6,356,514.9

1/293.46

France, Africa

Clarke 1866

6,378,206.4

6,356,583.8

1/294.98

North America

Airy 1830

6,377,563.4

6,356,256.9

1/299.32

Great Britain

Bessel 1841

6,377,397.2

6,356,079.0

1/299.15

Central Europe, Chile, Indonesia

Everest 1830

6,377,276.3

6,356,075.4

1/300.80

South Asia

Don't interpret the chart to mean there is radical disagreement

about the shape of the earth. The earth isn't perfectly round, it

isn't even a perfect ellipsoid, and slightly different shapes work

better for some regions than for the earth as a whole. The top

three are based on worldwide data and are truly global. Also, you

are very unlikely to find UTM grids based on any of the earlier

projections.

The most modern datums (jars my Latinist sensibilities since the

plural of datum in Latin

is data, but that has a different meaning

to us) are NAD83 and WGS84. These are based in turn on GRS80.

Differences between the three systems derive mostly from

redetermination of station locations rather than differences in the

datum. Unless you are locating a first-order station to

sub-millimeter accuracy (in which case you are way beyond the scope

of this page) you can probably regard them as essentially

identical.

I have no information on the NAD83 and WGS84 datums, nor can my

spreadsheet calculate differences between those datums and

WGS84.

Formulas For Converting Latitude and Longitude to UTM

These formulas are slightly modified from Army (1973). They are

accurate to within less than a meter within a given grid zone. The

original formulas include a now obsolete term that can be handled

more simply - it merely converts radians to seconds of arc. That

term is omitted here but discussed below.

Symbols

lat = latitude of point

long = longitude of point

long0 = central meridian of zone

k0 = scale along

long0 = 0.9996. Even though it's a

constant, we retain it as a separate symbol to keep the numerical

coefficients simpler, also to allow for systems that might use a

different Mercator projection.

e = SQRT(1-b2/a2) = .08 approximately. This

is the eccentricity of the earth's elliptical cross-section.

e'2 =

(ea/b)2 =

e2/(1-e2) = .007 approximately. The quantity

e' only occurs in even powers so it need only be calculated as

e'2.

n = (a-b)/(a+b)

rho =

a(1-e2)/(1-e2sin2(lat))3/2.

This is the radius of curvature of the earth in the meridian

plane.

nu = a/(1-e2sin2(lat))1/2. This is

the radius of curvature of the earth perpendicular to the meridian

plane. It is also the distance from the point in question to the

polar axis, measured perpendicular to the earth's surface.

p = (long-long0) in

radians(This differs from the treatment in

the Army reference)

Calculate the Meridional Arc

S is the meridional arc through the point in question (the distance

along the earth's surface from the equator). All angles are in

radians.

S = A'lat - B'sin(2lat) + C'sin(4lat) - D'sin(6lat) + E'sin(8lat),

where lat is in radians and

A' = a[1 - n + (5/4)(n2 -

n3) + (81/64)(n4 -

n5) ...]

B' = (3 tan/2)[1 - n + (7/8)(n2 -

n3) + (55/64)(n4 -

n5) ...]

C' = (15 tan2/16)[1 - n +

(3/4)(n2 -

n3) ...]

D' = (35 tan3/48)[1 - n +

(11/16)(n2 -

n3) ...]

E' = (315 tan4/512)[1 - n ...]

The USGS gives this form, which may be more appealing to some.

(They use M where the Army uses S)

M = a[(1 - e2/4 - 3e4/64 - 5e6/256

....)lat

- (3e2/8 + 3e4/32 +

45e6/1024...)sin(2lat) + (15e4/256 + 45e6/1024 +

....)sin(4lat)

- (35e6/3072 + ....) sin(6lat) +

....)] where lat is in radians

This is the hard part. Calculating the arc length of an ellipse

involves functions called elliptic

integrals, which don't reduce to neat closed formulas. So they

have to be represented as series.

Converting Latitude and Longitude to UTM

All angles are in radians.

y = northing = K1 + K2p2 +

K3p4, where

K1 = Sk0,

K2 = k0 nu sin(lat)cos(lat)/2 =

k0 nu sin(2 lat)/4

K3 = [k0 nu

sin(lat)cos3(lat)/24][(5 - tan2(lat) +

9e'2cos2(lat) +

4e'4cos4(lat)]

x = easting = K4p + K5p3, where

K4 = k0 nu cos(lat)

K5 = (k0 nu cos3(lat)/6)[1 -

tan2(lat) + e'2cos2(lat)]

Easting x is relative to the central meridian. For conventional UTM

easting add 500,000 meters to x.

What the Formulas Mean

The hard part, allowing for the oblateness of the Earth, is taken

care of in calculating S (or M). So K1 is simply the arc length

along the central meridian of the zone corrected by the scale

factor. Remember, the scale is a hair less than 1 in the middle of

the zone, and a hair more on the outside.

All the higher K terms involve nu, the local radius of curvature

(roughly equal to the radius of the earth or roughly 6,400,000 m),

trig functions, and powers of e'2 ( =

.007 ). So basically they are never much larger than nu. Actually

the maximum value of K2 is about nu/4 (1,600,000), K3 is about

nu/24 (267,000) and K5 is about nu/6 (1,070,000). Expanding the

expressions will show that the tangent terms don't affect

anything.

If we were just to stop with the K2 term in the northing, we'd have

a quadratic in p. In other words, we'd approximate the parallel of

latitude as a parabola. The real curve is more complex. It will be

more like a hyperbola equatorward of about 45 degrees and an

ellipse poleward, at least within the narrow confines of a UTM

zone. (At any given latitude we're cutting the cone of latitude

vectors with an inclined plane, so the resulting intersection will

be a conic section. Since the projection cylinder has a curvature,

the exact curve is not a conic but the difference across a

six-degree UTM zone is pretty small.) Hence the

need for higher order terms. Now p will never be more than +/-3

degrees = .05 radians, so p2 is always

less than .0025 (1/400) and p4 is

always less than .00000625 (1/160000). Using a spreadsheet, it's

easy to see how the individual terms vary with latitude.

K2p2 never exceeds 4400 and

K3p4 is at most a bit over 3. That is,

the curvature of a parallel of latitude across a UTM zone is at

most a little less than 4.5 km and the maximum departure from a

parabola is at most a few meters.

K4 is what we'd calculate for easting in a simple-minded way, just

by calculating arc distance along the parallel of latutude. But, as

we get farther from the central meridian, the meridians curve

inward, so our actual easting will be less than K4. That's what K5

does. Since p is never more than +/-3 degrees = .05 radians,

p3 is always less than .000125

(1/8000). The maximum value of K5p3 is

about 150 meters.

That Weird Sin 1" Term in the Original Army Reference

The Army reference defines p in seconds of arc and includes a sin

1" term in the K formulas. The Sin 1" term is a holdover from the

days when this all had to be done on mechanical desk calculators

(pre-computer) and terms had to be kept in a range that

would retain sufficient precision at intermediate

steps. For that small an angle the difference between sin 1" and 1"

in radians is negligible. If p is in seconds of arc, then (psin 1")

merely converts it to radians.

The sin 1" term actually included an extra factor of 10,000, which

was then corrected by multiplying by large powers of ten

afterward.

The logic is a bit baffling. If I were doing this on a desk

calculator, I'd factor out as many terms as possible rather than

recalculate them for each term. But perhaps in practice the

algebraically obvious way created overflows or underflows, since

calculators could only handle limited ranges.

In any case, the sin1" term is not needed any more. Calculate p in

radians and omit the sin1" terms and the large power of ten

multipliers.

Converting UTM to Latitude and Longitude

y = northing, x = easting (relative to central meridian; subtract

500,000 from conventional UTM coordinate).

Calculate the Meridional Arc

This is easy: M = y/k0.

Calculate Footprint Latitude

mu = M/[a(1 - e2/4 - 3e4/64 -

5e6/256...)

e1 = [1 - (1 -

e2)1/2]/[1 + (1 -

e2)1/2]

footprint latitude fp = mu + J1sin(2mu) +

J2sin(4mu) + J3sin(6mu) + J4sin(8mu), where:

J1 = (3e1/2 - 27e13/32 ..)

J2 = (21e12/16 -

55e14/32 ..)

J3 = (151e13/96 ..)

J4 = (1097e14/512 ..)

Calculate Latitude and Longitude

e'2 =

(ea/b)2 =

e2/(1-e2)

C1 = e'2cos2(fp)

T1 = tan2(fp)

R1 =

a(1-e2)/(1-e2sin2(fp))3/2.

This is the same as rho in the forward conversion formulas above,

but calculated for fp instead of lat.

N1 = a/(1-e2sin2(fp))1/2. This is

the same as nu in the forward conversion formulas above, but

calculated for fp instead of lat.

D = x/(N1k0)

lat = fp - Q1(Q2 - Q3 + Q4), where:

Q1 = N1 tan(fp)/R1

Q2 = (D2/2)

Q3 = (5 + 3T1 + 10C1 -

4C12 -9e'2)D4/24

Q4 = (61 + 90T1 + 298C1 +45T12 -

3C12 -252e'2)D6/720

long = long0 + (Q5 - Q6 + Q7)/cos(fp), where:

Q5 = D

Q6 = (1 + 2T1 + C1)D3/6

Q7 = (5 - 2C1 + 28T1 - 3C12 +

8e'2 +

24T12)D5/120

What Do The Formulas Mean?

As the sketch above shows, because of the poleward curve of

parallels, the footprint latitude is always greater than the true

latitude. Q1 is just a scaling coefficient and is constant for any

given fp. The tangent term basically means the closer to the pole

you are, the faster the parallels curve. Q2 is a quadratic term in

x. Again, as with converting from geographic coordinates to UTM, we

approximate the parallel as a parabola and add higher order

corrections.

To determine longitude, we make a simple minded approximation that

longitude is proportional to easting, but then, since fp is too

large, the true longitude is smaller, since it lies on a parallel

closer to the the equator. The divisor cos(fp) corrects for the

varying length of degrees of longitude as latitude varies.

Before linking to the program, note:

It is an Excel spreadsheet. If you have Excel on your computer, it

will (or should) open when you click the link. Most major

spreadsheet programs can read spreadsheets in other formats.

A spreadsheet is not an applet or program. In particular, you can't

manually enter data into a cell and preserve any formulas that are

there. That's why some aspects of data entry are clunkier than they

might otherwise be with, say, a Visual Basic program.

There are three computation pages, one for single conversions, the

other two for batch conversions. Other pages contain information on

datums and the specific conversion formulas. To use the batch

conversions you need to be somewhat proficient in spreadsheets as

you will have to copy data and cell formulas.

For our mutual peace of mind, run a virus check.

You may copy the program for your own non-commercial use and for

non-commercial distribution to others, but not for commercial use.

Please give appropriate credit when citing your calculations. You

may also modify it as needed for your personal use. In Internet

Explorer, right click on the link and select Save Target As... to

save the spreadsheet to your own disk.

References

Snyder, J. P., 1987; Map Projections - A Working Manual. U.S.

Geological Survey Professional Paper 1395, 383

p. If you are at all serious about maps you

need this book.

Army, Department of, 1973; Universal Transverse Mercator Grid, U.

S. Army Technical Manual TM 5-241-8, 64 p.

NIMA Technical Report 8350.2, "Department of Defense World Geodetic

System 1984, Its Definition and Relationships with Local Geodetic

Systems," Second Edition, 1 September 1991 and its supplements. The

report is available from the NIMA Combat Support Center and its

stock number is DMATR83502WGS84. Non-DoD requesters may obtain the

report as a public sale item from the U.S. Geological Survey, Box

25286, Denver Federal Center, Denver, Colorado 80225 or by phone at

1-800-USA-MAPS.

Return to Professor

Dutch's Home Page

Created 12 September 2003, Last Update 14

December 2009

Not an official UW Green Bay site