مشکلم حل شد
عنوان دقیق مسئله تبدیل مختصات latLong به UTM هستش که از فرمول زیر میشه اینکارو انجام داد:
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.
* 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.
از این توضیحات سایت IBM
هم برای برنامه نویسیش میتونید استفاده کنید.