| UTM conversion [message #22342] |
Tue, 07 November 2000 00:00  |
Kelly Dean
Messages: 92
Registered: March 1997
|
Member |
|
|
Is there an IDL routine available to convert UTM ( easting and northing
) into earth (lat/long) coordinates?
Kelly Dean
|
|
|
|
| Re: UTM conversion [message #22424 is a reply to message #22342] |
Wed, 08 November 2000 00:00  |
Pavel A. Romashkin
Messages: 531
Registered: November 2000
|
Senior Member |
|
|
I think I have it. Is it it, below?
Cheers,
Pavel
> Luis Alonso wrote:
>
>> Hello!
>>
>> I need to change coordinates from Latitude/Longitude to UTM.
>> I know that with MAP one can proyect on either one, but i don't know how to
>> transform a vector from one into the other.
>>
>> Thanks
>> Luis Alonso
>
> Hello,
>
> The book by Snyder (referenced below) is very much worth the $40US. I have
> used the following code
> to convert Lon/Lat to UTM in meters.
>
> Ben
>
> --------SNIP-------
> ;+
> ; NAME:
> ; ll_to_utm
> ;
> ; PURPOSE:
> ; This function converts latitude/longitude measurements to UTM coordinates.
> ; Two map data are are available: NAD27 and WGS84.
> ;
> ; The formulae come from: J.P. Snyder, "Map projections - A working
> manual',1987,
> ; U.S.G.S. Professional Paper 1395, Supt. of Docs No: I 19.16:1395,
> ; U.S. Govt Prinitng Office,Washington, DC 20402.
> ;
> ; Constants used in formulae come from this source also unless otherwise
> noted.
> ;
> ;
> ; CATEGORY:
> ; Mapping.
> ;
> ; CALLING SEQUENCE:
> ; Result = ll_to_utm(lon,lat)
> ;;
> ; INPUTS:
> ; LAT/LON Identically sized vectors containning the decimal degrees of
> latitude and longitude.
> ; Degrees west and south are negative.
> ;
> ; OPTIONAL INPUTS:
> ; None.
> ;
> ; KEYWORD PARAMETERS:
> ; DATUM Set this keyword to 1 for NAD27 (default) and to 2 for WGS84.
> ; NOFALSE Set this keyword to prevent the return of False Easting. Default is
> to
> ; return a False Easting (NoFalse = 0).
> ;
> ; OUTPUTS:
> ; A 2,n elements double precision array where [0,*] = Easting
> ; and [1,*] = False Northing.
> ;
> ; OPTIONAL OUTPUTS:
> ; None.
> ;
> ; COMMON BLOCKS:
> ; None.
> ;
> ; SIDE EFFECTS:
> ; None known.
> ;
> ; RESTRICTIONS:
> ; None known.
> ;
> ; EXAMPLE:
> ; Lat =double( [44.0551666000000000, 44.0551999,44.0552166])
> ; Lon =double( [-69.5231832999999990, -69.5231999,-69.5232166])
> ; UTM = ll_to_utm(Lon, Lat, Datum = 2)
> ; For i = 0,2 Do Print, Lon[i],UTM[0,i], Lat[i],UTM[1,i]
> ; -69.523186 458092.92 44.055168 4878133.2
> ; -69.523201 458091.72 44.055199 4878136.6
> ; -69.523216 458090.51 44.055218 4878138.8
>
> ;
> ; MODIFICATION HISTORY:
> ; Written by: Ben Tupper Spring 1998
> ; Pemaquid River Company
> ; email: pemaquidriver@tidewater.net
> ; tel: (207) 563 - 1048
> ; 248 Lower Round Pond Road
> ; POB 106
> ; Bristol, ME 04539-0106
> ;
> ; SEP 1999 Cleaned up documentation and added NoFalse keyword.
> ; DEC 2, 1999 Reversed order of input and output arguments to fit familiar
> x,y ordered pair sense.
> ; 4 FEB 2000 Documentation change on 2 DEC 1999.
> ;-
>
> FUNCTION ll_to_utm, lonV,latV,datum=datum, NoFalse = NoFalse
>
> On_Error,2
>
> Num = n_elements(LatV) & Utm = DblArr(2,Num)
> Lon=Double(LonV) & Lat=Double(LatV)
>
> If n_elements(NoFalse) Eq 0 Then NoFalse = 0
> If n_elements(Datum) EQ 0 Then Datum = 1
> Data = ['Blank','NAD27','WGS84']
>
> ;A emimajor axis of the of ellipsoid earth
> ;eccsq eccentricity of ellipsoid (squared)
> ;k0 relative scale factor along parallel of latitude
> ;lon longitude in +/-ddd.dddddd
> ;lat latitude in +/-dd.ddddddd
> ;X the rectangualar coordinate output...measured in meters from central
> meridian
> ; this is converted to a false easting by adding half the distance between
> ; the bounding meridians
> ;Y the rectangular coordinate output.... measured in meters from the equator
>
> Case Data[datum] of
>
> 'NAD27': BEGIN
> ; based on Clarke 1866 (aka NAD27) or
> ; assign constant values for NAD 27 to UTM
> A = 6378206.4d
> eccsq = 0.00676866d
> k0 = 0.9996d
>
> END
>
> 'WGS84':BEGIN
> ; WGS84 (Aka NAD83)
> ; from Chuck Gantz via http://gpsy.com/gpsinfo/geotoutm.htm
> ; assign constant values for WGS 84 datum for lat/lon to UTM
> A = 6378137.d
> eccsq = 0.00669439d
> k0 = 0.9996d
>
> END
>
> EndCase ;which datum
>
> ;determine central meridian (lambda0 in Snyder)
> central_meridian = double(long(lon/6)*6+3*long(lon/6.)/abs(long(lon/6.)))
>
> ;determine lon origin in radians
> lon_orig = central_meridian*!dtor
>
> ; convert to radians
> lon = lon*!dtor & lat = lat*!dtor
>
> ;
> eccsq_prime =eccsq/(1.-eccsq)
>
> N = a/(1.-eccsq*sin(lat)*sin(lat))^.5
>
> T = tan(lat)*tan(lat)
>
> C = eccsq_prime * cos(lat)*cos(lat)
>
> AA = cos(lat)*(lon-lon_orig)
>
> M = a*((1-eccsq/4.-3.*eccsq^2/64.-5.*eccsq^3/256.)*lat - $
> (3.*eccsq/8. + 3.*eccsq^2/32.+45.*eccsq^3/1024)*sin(2.*lat)+$
> (15.*eccsq^2/256.+45.*eccsq^3/1024.)*sin(4.*lat)- $
> (35.*eccsq^3/3072.)*sin(6.*lat))
>
> ; k is unused for UTM... default to k0
> ;k = k0*(1.d+ (1.d+C)*A^2/2.d + $
> ; (5.d - 4.d*T + 42.d*C + 13.d*C^2 -28.d*eccsq_prime)*A^4/24.d + $
> ; (61.d - 148.d*T + 16.d*T^2)*A^6/720.d)
>
> ; Easting or X in Snyder
> UTM[0,*]= k0*N*(AA+(1.-T-C)*AA^3/6.d $
> + (5.d - 18.d*T + T^2 + 72.d*C - 58.d*eccsq_prime)*AA^5/120.d)
>
> ; make False Easting if desired (default)
> If NoFalse EQ 0 Then UTM[0,*] = UTM[0,*] + 500000.0d
>
> ;Northing or Y in Snyder
> UTM[1,*] = k0*(M+N*tan(lat)*( AA^2 / 2.d + ( 5.d - T + 9.d *C + 4.d * C^2 ) *
> $
> $ AA^4/24.d + (61d - 58d*T + T^2 + 600d*C - 330d*eccsq_prime)*
> AA^6 / 720.d ))
>
> Return,UTM
>
> end
>
> ------SNIP------
>
> --
> Ben Tupper
>
> Bigelow Laboratory for Ocean Science
> tupper@seadas.bigelow.org
>
> pemaquidriver@tidewater.net
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
