Navigation Menu

Skip to content

Instantly share code, notes, and snippets.

@greggirwin

greggirwin/r2-mapping-lib.r

Created July 18, 2019 18:09
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save greggirwin/2380315ed4912a2801b581ca3e6ada0e to your computer and use it in GitHub Desktop.
Save greggirwin/2380315ed4912a2801b581ca3e6ada0e to your computer and use it in GitHub Desktop.
Download ZIP
R2 Mapping Lib
Raw
r2-mapping-lib.r
REBOL []
;!! This is very much a work in progress, with many things untested. It needs
; to be clear what is functional and break out R&D bits into a working script
; for that purpose.
;-------------------------------------------------------------------------------
Arcsin: func [X [decimal!]] [
;'Arcsin(X) = Atn(X / Sqr(-X * X + 1)) [from MS Help - VB4, 1995]
arctangent (X / square-root (-X * X + 1))
]
as-radians: func [
"Convert a value given in degrees to radians."
degrees [number!]
] [
degrees * (pi / 180)
]
atan2: func [
"Arctangent2"
x
y
][
either zero? y [
either negative? X [pi] [0]
][
2.0 * arctangent ( Y / ( X + square-root ((X * X) + (Y * Y)) ) )
]
]
incr: func [
{Increment a value by 1.}
word [word!]
/by {Change by this amount} ; /skip ?
value
][
set word add get word any [value 1]
]
limit: func [
"Make sure val falls between lower and upper bounds, inclusive"
val lower upper
][
max lower min val upper
]
; http://math2.org/math/algebra/functions/sincos/properties.htm
; 'sin^2(X) = (1 - cos(2X))/2 [from Greer and Hancox, 1991]
sine2: func [x [decimal!]] [
power sine x 2 ; This is the fastest of these alternatives
;(1 - cosine (2 * x)) / 2
;(sine x) * (sine x)
;multiply sine x sine x
]
;-------------------------------------------------------------------------------
; For testing and debugging.
show-ddmm.mmmm-to-sexa: func [lat lon /local lat-d lon-d] [
print [tab lat lon]
print [tab "=" lat-d: ddmm.mmmm-to-dd.dddd lat lon-d: ddmm.mmmm-to-dd.dddd lon]
print [tab "=" decimal-to-sexagesimal lat-d decimal-to-sexagesimal lon-d]
]
; http://woodshole.er.usgs.gov/operations/sfmapping/gps.htm
; http://www.circuitcellar.com/library/print/1000/Stefan123/4.htm
; NMEA = National Marine Electronics Association
; see http://www.gpsinformation.org/dale/nmea.htm
; http://www.earthpoint.us/ExcelToKml.aspx
; http://www.tma.dk/gps/
; http://transition.fcc.gov/mb/audio/bickel/DDDMMSS-decimal.html
; Google maps API uses a pure number of degrees.
; Split your GPS values into (d)dd and mm.mmmm.
; Divide the mm.mmmm by 60 and add it to the ddd.
; Latitude and longitude data received from a GPS receiver in NMEA-0183 format is
; in units ddmm.mmmm, where dd equals degrees, mm equals minutes, and .mmmm is
; decimal minutes.
; 1. Separate and save dd from the incoming latitude and longitude
; 2. Divide mm.mmmm by 60, yielding 0.dddd
; 3. Add the saved dd to 0.dddd, yielding dd.dddd
dd.dddd-to-ddmm.mmmm: func [value [decimal!]] [
deg: round/down value
min-sec: abs (value - deg * 60) ;* remainder value 1
min: round/down min-sec
;sec: abs round 60 * remainder min-sec 1
sec: min-sec - min * 60
rejoin [deg pad-num min 2 "." pad-num sec 4]
]
dd.dddd-to-radians: func [dd.dddd] [
dd.dddd / 57.2957795
]
; ddmm.mmmm = NMEA = WGS84
ddmm.mmmm-to-dd.dddd: func [
"Convert lat/lon in NMEA/WGS864/ddmm.mmmm format to decimal degrees (dd.dddd)."
ddmm.mmmm [decimal!]
] [
deg: round/down ddmm.mmmm / 100
deg + ((ddmm.mmmm - (deg * 100)) / 60)
]
dms-compass-dir: func [
"Return a compass direction (NSEW) given a lat/lon degree value."
degrees [number!] "Negative degrees are South and West; positive values are North and East."
direction [word!] "'lat or 'lon"
][
pick select [lat "SN" lon "WE"] direction negative? degrees
]
decimal-to-sexagesimal: func [
value [decimal!]
/compass "Return the compass direction, rather than a signed value for degrees."
direction [word!] "'lat or 'lon"
/local deg min-sec min sec
] [
deg: round/down value
min-sec: abs (60 * remainder value 1)
min: round/down min-sec
sec: round 60 * remainder min-sec 1
;!! The Masculine Ordinal Indicator "º" was used in the first case here.
; I changed it to the degree sign after some research. I don't know
; why it was there originally, unless it was a character code typo.
rejoin either compass [
[abs deg "°" min {'} sec {"} dms-compass-dir deg direction]
][
[deg "°" min {'} sec {"}]
]
]
decimal-to-dms: :decimal-to-sexagesimal
dd.dddd-to-dms: :decimal-to-sexagesimal
ddmm.mmmm-to-dms: func [
"Convert lat/lon in NMEA/WGS864/ddmm.mmmm format to deg/min/sec."
ddmm.mmmm [decimal!]
] [
dd.dddd-to-dms ddmm.mmmm-to-dd.dddd ddmm.mmmm
]
ddmm.mmmm-to-dms: func [
"Convert lat/lon in NMEA/WGS864/ddmm.mmmm format to deg/min/sec."
ddmm.mmmm [decimal!]
/compass "Return the compass direction, rather than a signed value for degrees."
direction [word!] "'lat or 'lon"
] [
either compass [
dd.dddd-to-dms/compass ddmm.mmmm-to-dd.dddd ddmm.mmmm direction
][
dd.dddd-to-dms ddmm.mmmm-to-dd.dddd ddmm.mmmm
]
]
nautical-miles-to-statute-miles: func [naut-miles] [naut-miles * 1.150779]
radians-to-dd.dddd: func [radians] [radians * 57.2957795]
radians-to-natuical-miles: func [radians] [radians * 3437.7387]
sexagesimal-to-decimal: func [value [string!] /local deg min-sec min sec sign] [
dig=: charset "0123456789"
non-dig=: complement dig=
sign: 1
parse/all value [
opt [#"-" (sign: -1)]
copy deg some dig=
some non-dig=
copy min some dig=
some non-dig=
copy sec some dig=
any non-dig=
]
either all [deg min sec not empty? deg not empty? min not empty? sec] [
deg: to integer! deg
min: to integer! min
sec: to integer! sec
deg + (min / 60) + (sec / 3600) * sign
][
;print "Did you give me degrees, minutes, and seconds? I don't think so."
none
]
]
statue-miles-to-feet: func [stat-miles] [stat-miles * 5'280]
statute-miles-to-nautical-miles: func [stat-miles] [stat-miles / 1.150779]
;-------------------------------------------------------------------------------
; Steeve
; Project: func [
; {Orthogonal projection of a point P on a line AB, return coordinates [x y]}
; ax ay bx by px py
; /local sx sy ux uy ratio
; ][
; sx: bx - ax
; sy: by - ay
; ux: px - ax
; uy: py - ay
; ratio: sx * ux + (sy * uy) / (sx * sx + (sy * sy))
; reuse [ratio * sx + ax ratio * sy + ay]
; ]
;-------------------------------------------------------------------------------
WGS84-to-Google-Bing: func [
lon [decimal!]
lat [decimal!]
/local x y
] [
x: lon * 20037508.34 / 180
y: log-e (tangent ((90 + lat) * PI / 360)) / (PI / 180)
y: y * 20037508.34 / 180
reduce [x y]
]
Google-Bing-to-WGS84-Mercator: func [
x [decimal!]
y [decimal!]
/local lat lon
] [
lon: (x / 20037508.34) * 180
lat: (y / 20037508.34) * 180
lat: 180 / PI * (2 * arctangent (exp (lat * PI / 180)) - PI / 2)
reduce [lon lat]
]
; Bing Tile Maps info: http://msdn.microsoft.com/en-us/library/bb259689.aspx
; This logic should be largely compatible across mapping services, as they all
; use compatible or similar systems. Need to verify and call out where there
; are differences, then generalize this and make naming generic where possible.
; Bing uses a spherical mercator projection
bing-tile-maps: context [
; Equatorial map table
;"Level of Detail" "Map Width and Height (pixels)" "Ground Resolution (meters / pixel)" "Map Scale (at 96 dpi)"
equatorial-map-table: [
1 512 78'271.5170 "1:295,829,355.45"
2 1'024 39'135.7585 "1:147,914,677.73"
3 2'048 19'567.8792 "1:73,957,338.86"
4 4'096 9'783.9396 "1:36,978,669.43"
5 8'192 4'891.9698 "1:18,489,334.72"
6 16'384 2'445.9849 "1:9,244,667.36"
7 32'768 1'222.9925 "1:4,622,333.68"
8 65'536 611.4962 "1:2,311,166.84"
9 131'072 305.7481 "1:1,155,583.42"
10 262'144 152.8741 "1:577,791.71"
11 524'288 76.4370 "1:288,895.85"
12 1'048'576 38.2185 "1:144,447.93"
13 2'097'152 19.1093 "1:72,223.96"
14 4'194'304 9.5546 "1:36,111.98"
15 8'388'608 4.7773 "1:18,055.99"
16 16'777'216 2.3887 "1:9,028.00"
17 33'554'432 1.1943 "1:4,514.00"
18 67'108'864 0.5972 "1:2,257.00"
19 134'217'728 0.2986 "1:1,128.50"
20 268'435'456 0.1493 "1:564.25"
21 536'870'912 0.0746 "1:282.12"
22 1'073'741'824 0.0373 "1:141.06"
23 2'147'483'648 0.0187 "1:70.53"
]
EARTH-RADIUS: 6'378'137 ; In meters
MIN_LATITUDE: -85.05112878 ; Limits are an area where systems may differ.
MAX_LATITUDE: 85.05112878
MIN_LONGITUDE: -180
MAX_LONGITUDE: 180
LEVEL_OF_DETAIL_MIN: 0
LEVEL_OF_DETAIL_MAX: 23
clip: :limit
clip-latitude: func [
latitude [decimal!] "Latitude (in degrees)"
][
clip latitude MIN_LATITUDE MAX_LATITUDE
]
clip-longitude: func [
longitude [decimal!] "Longitude (in degrees)"
][
clip longitude MIN_LONGITUDE MAX_LONGITUDE
]
map-size: func [
"Returns the map width and height in pixels."
level [integer!] "Level of detail, from 1 (lowest detail) to 23 (highest detail)"
] [
2 ** level * 256
]
; Manually inspect results
; repeat i 23 [print [i map-size i]]
;Some specific tests:
; 512 = map-size 1
; 2'097'152 = map-size 13
; 134'217'728 = map-size 19
; 2'147'483'648 = map-size 23
ground-resolution: func [
"Returns the ground resolution (in meters per pixel) at a specified latitude and level of detail."
latitude [decimal!] "Latitude (in degrees) at which to measure the ground resolution."
level [integer!] "Level of detail, from 1 (lowest detail) to 23 (highest detail)"
] [
latitude: clip-latitude latitude
(cosine (latitude * PI / 180)) * 2 * PI * EARTH-RADIUS / map-size level
]
map-scale: func [
"Determines the map scale at a specified latitude, level of detail, and screen resolution."
; Returns the map scale, expressed as the denominator N of the ratio 1 : N.
latitude [decimal!] "Latitude (in degrees) at which to measure the ground resolution."
level [integer!] "Level of detail, from 1 (lowest detail) to 23 (highest detail)"
dpi [integer!] "Resolution of the screen, in dots per inch"
][
(ground-resolution latitude level) * dpi / 0.0254
]
lat-long-to-pixel: func [
"Converts a point from latitude/longitude WGS-84 coordinates (in degrees) into pixel XY coordinates at a specified level of detail."
latitude [decimal!] "Latitude of the point, in degrees."
longitude [decimal!] "Longitude of the point, in degrees."
level [integer!] "Level of detail, from 1 (lowest detail) to 23 (highest detail)"
/local x y sin-latitude map-sz
][
latitude: clip-latitude latitude
longitude: clip-longitude longitude
x: (longitude + 180) / 360
sin-latitude: sine/radians (latitude * PI / 180)
y: 0.5 - ((log-e ((1 + sin-latitude) / (1 - sin-latitude))) / (4 * PI))
map-sz: map-size level
x: round clip x * map-sz + 0.5 0 map-sz - 1
y: round clip y * map-sz + 0.5 0 map-sz - 1
as-pair x y
]
lat-long-to-pixel-XY: func [
"Converts a point from latitude/longitude WGS-84 coordinates (in degrees) into pixel XY coordinates at a specified level of detail."
latitude [decimal!] "Latitude of the point, in degrees."
longitude [decimal!] "Longitude of the point, in degrees."
level [integer!] "Level of detail, from 1 (lowest detail) to 23 (highest detail)"
/local pair
][
pair: lat-long-to-pixel latitude longitude level
reduce [pair/x pair/y]
]
pixel-XY-to-lat-long: func [
"Converts a pixel from pixel XY coordinates at a specified level of detail into latitude/longitude WGS-84 coordinates (in degrees)."
x [integer!] "X coordinate of the point, in pixels."
y [integer!] "Y coordinate of the point, in pixels."
level [integer!] "Level of detail, from 1 (lowest detail) to 23 (highest detail)"
/local map-sz
][
map-sz: map-size level
x: ((clip x 0 map-sz - 1) / map-sz) - 0.5
y: 0.5 - ((clip y 0 map-sz - 1) / map-sz)
reduce [
'latitude 90 - (360 * (arctangent/radians (exp ((negate y) * 2 * PI))) / PI)
'longitude 360 * x
]
]
pixel-to-tile: func [
"Converts pixel XY coordinates into tile XY coordinates of the tile containing the specified pixel."
coord [pair!] "Pixel XY coordinate."
][
as-pair coord/x / 256 coord/y / 256
]
pixel-XY-to-tile-XY: func [
"Converts pixel XY coordinates into tile XY coordinates of the tile containing the specified pixel."
x [integer!] "Pixel X coordinate."
y [integer!] "Pixel Y coordinate."
][
reduce [x / 256 y / 256]
]
tile-XY-to-pixel: func [
"Converts tile XY coordinates into pixel XY coordinates of the upper-left pixel of the specified tile."
x [integer!] "Tile X coordinate."
y [integer!] "Tile Y coordinate."
][
as-pair x * 256 y * 256
]
tile-XY-to-pixel-XY: func [
"Converts tile XY coordinates into pixel XY coordinates of the upper-left pixel of the specified tile."
x [integer!] "Tile X coordinate."
y [integer!] "Tile Y coordinate."
][
reduce [x * 256 y * 256]
]
tile-to-quad-key: func [
"Converts tile XY coordinates into a QuadKey (string) at a specified level of detail."
coord [pair!] "Tile XY coordinate."
level [integer!] "Level of detail, from 1 (lowest detail) to 23 (highest detail)"
/local quad-key digit mask
][
tile-XY-to-quad-key coord/x coord/y level
]
; Manually tested to level 3. Produces correct results.
tile-XY-to-quad-key: func [
"Converts tile XY coordinates into a QuadKey (string) at a specified level of detail."
x [integer!] "Tile X coordinate."
y [integer!] "Tile Y coordinate."
level [integer!] "Level of detail, from 1 (lowest detail) to 23 (highest detail)"
/local quad-key digit mask
][
quad-key: copy ""
for i level 1 -1 [ ; FOR is slow. May want to optimize this.
digit: #"0"
mask: shift/left 1 (i - 1) ; This is REBOL's native SHIFT, in newer versions.
if not zero? (X and mask) [
incr 'digit
]
if not zero? (Y and mask) [
incr/by 'digit 2
]
append quad-key digit
]
quad-key
]
; Manually tested to level 3. Produces correct results.
quad-key-to-tile-XY: func [
"Converts a QuadKey into tile XY coordinates."
quad-key [any-string!] "QuadKey of the tile."
/local x y level mask
][
x: y: 0
level: length? quad-key
for i level 1 -1 [ ; FOR is slow. May want to optimize this.
mask: shift/left 1 (i - 1) ; This is REBOL's native SHIFT, in newer versions.
switch/default digit: pick quad-key level - i + 1 [
#"0" []
#"1" [x: x OR mask]
#"2" [y: y OR mask]
#"3" [
x: x OR mask
y: y OR mask
]
][
print mold digit
make error! "Invalid QuadKey digit sequence."
]
]
reduce [
'tile reduce [x y] ;as-pair x y
'level level
]
]
]
;-------------------------------------------------------------------------------
; UTM/USNG
; By convention, the WGS 84 geoid describe Earth as an ellipsoid along
; North-South axis with an equatorial radius of a = 6378.137km and an eccentricity
; of e = 0.0818192.
; Latitude/longitude-to-UTM test values
; [ 0.0000 0.0000 ] "31 N 166021 0"
; [ 0.1300 -0.2324 ] "30 N 808084 14385"
; [-45.6456 23.3545 ] "34 G 683473 4942631"
; [-12.7650 -33.8765 ] "25 L 404859 8588690"
; [-80.5434 -170.6540] "02 C 506346 1057742"
; [ 90.0000 177.0000] "60 Z 500000 9997964"
; [-90.0000 -177.0000] "01 A 500000 2035"
; [ 90.0000 3.0000 ] "31 Z 500000 9997964"
; [ 23.4578 -135.4545] "08 Q 453580 2594272"
; [ 77.3450 156.9876] "57 X 450793 8586116"
; [-89.3454 -48.9306 ] "22 A 502639 75072"
; http://www.ngs.noaa.gov/cgi-bin/usng_getus.prl
; Latitude Longitude Datum USNG
; ---------- ----------- ------ -----------------------
; N385222.08 W0770206.86 NAD 83 -> 18SUJ2344204629(NAD 83) 1M
; N385222.1 W0770206.9 NAD 83 -> 18SUJ23440463(NAD 83) 10M
; N385222. W0770206. NAD 83 -> 18SUJ234046(NAD 83) 100M
; The combination of a zone and a latitude band defines a grid zone. The zone is
; always written first, followed by the latitude band.
; A position on the Earth is referenced in the UTM system by the UTM zone, and the
; easting and northing coordinate pair. The easting is the projected distance of
; the position eastward from the central meridian, while the northing is the
; projected distance of the point north from the equator (in the northern
; hemisphere). Eastings and northings are measured in meters. The point of origin
; of each UTM zone is the intersection of the equator and the zone's central
; meridian. In order to avoid dealing with negative numbers, the central meridian
; of each zone is given a "false easting" value of 500,000 meters.
; A USNG designation is comprised of 3 components; a grid zone designation (GZD), a 100 kilometer
; (100K) square designation, and a varying number of even digits representing an X and Y coordinate
; (XY) within the GZD and 100K square. The 100,000 and 1,000,000 place digits are replaced by the
; GZD and 100K characters. This means that the XY portion begins with the 10,000 place of the UTM
; coordinate. More digits represent increasing precision and are always paired and grouped as an X
; component and a Y component (XXXYYY, not XYXYXY).
UTM-zone-prefix=: [2 digit] ; 00-59
UTM-zone=: [1 char!]
UTM-specification: [UTM-zone-prefix= UTM-zone= 2 14 digit]
UTM-1-meter: [UTM-zone-prefix= UTM-zone= 14 digit]
UTM-10-meter: [UTM-zone-prefix= UTM-zone= 12 digit]
UTM-100-meter: [UTM-zone-prefix= UTM-zone= 10 digit]
UTM-1-km: [UTM-zone-prefix= UTM-zone= 8 digit]
UTM-10-km: [UTM-zone-prefix= UTM-zone= 6 digit]
; http://www.igorexchange.com/node/927
;-------------------------------------------------------------------------------
;-- UTM Zone
UTM-zone-letter: UTM-lat-zone: func [
{Return the UTM letter designator for the given latitude, "Z" if latitude is outside the UTM limits of 84N to 80S.}
lat [number!]
][
;?? Do we want to throw an error for a bad value, or return a known invalid zone letter?
; If we throw an error, it catches the mistake sooner, but makes them TRY in the caller.
; If we don't throw an error, the bad value propagates until someone chokes or handles it.
; if any [lat < -80 lat > 84] [
; make error! join [script invalid-arg] lat
; ]
either all [lat >= -80 lat <= 84 ] [
;!! Note the extra X at the end of the string, which accounts for the
; X zone being wider (72-84). Without that, the forumla runs off the
; end of the string and returns none for 80-84º values.
form pick "CDEFGHJKLMNPQRSTUVWXX" add 1 (round/floor (lat + 80)) / 8
][
; The latitude is outside the UTM limits.
"Z"
]
]
; Tests
; for i -80 -8 8 [
; print [
; i - .01 UTM-zone-letter i - .01
; tab
; i UTM-zone-letter i
; tab
; i + .01 UTM-zone-letter i + .01
; ]
; ]
; for i 0 88 8 [
; print [
; i - .01 UTM-zone-letter i - .01
; tab
; i UTM-zone-letter i
; tab
; i + .01 UTM-zone-letter i + .01
; ]
; ]
; Long -= int((Long+180)/360)*360; //ensure longitude within -180.00..179.9
UTM-zone-number: UTM-lon-zone: func [
{Return the UTM number designator for the given longitude.}
lat [number!] "Latitude. Needed due to special Svalbard zones."
lon [number!] "Longitude. Valid range is -180 to < 180"
/local zone
][
;?? Do we want to throw an error for a bad value, or return a known invalid zone letter?
; If we throw an error, it catches the mistake sooner, but makes them TRY in the caller.
; If we don't throw an error, the bad value propagates until someone chokes or handles it.
; if any [lon < -180 lon >= 180] [
; make error! join [script invalid-arg] lon
; ]
;?? Do we want to normalize the longitude? If so, we should handle all values.
; This works for -180.. values, but < -180 aren't changed.
; Normalize longitude to -180.00..179.999999999999 range.
;lon: lon - ((round/down ((lon + 180) / 360)) * 360)
either all [lon >= -180 lon < 180] [
zone: round/down ((lon + 180) / 6) + 1
; Svalbard Zone checks
if all [lat >= 56 lat < 64 lon >= 3 lon < 12] [
zone: 32
]
if all [lat >= 72.0 lat < 84.0 lon >= 0] [
zone: case [
all [lon >= 0 lon < 9] 31
all [lon >= 9 lon < 21] 33
all [lon >= 21 lon < 33] 35
all [lon >= 33 lon < 42] 37
; this is our 'else case
true [round/down ((lon + 180) / 6) + 1]
]
]
zone
][
; The longitude is outside the UTM limits.
0
]
]
; Tests
; for i -180 -6 6 [
; print [
; i - .01 attempt [UTM-zone-number 0.0 i - .01]
; tab
; i UTM-zone-number 0.0 i
; tab
; i + .01 UTM-zone-number 0.0 i + .01
; ]
; ]
; for i 0 180 6 [
; print [
; i - .01 UTM-zone-number 0.0 i - .01
; tab
; i attempt [UTM-zone-number 0.0 i]
; tab
; i + .01 attempt [UTM-zone-number 0.0 i + .01]
; ]
; ]
; ; Svalbard checks
; for i 0 12 3 [
; print [
; i - .01 UTM-zone-number 60.0 i - .01
; tab
; i attempt [UTM-zone-number 60.0 i]
; tab
; i + .01 attempt [UTM-zone-number 60.0 i + .01]
; ]
; ]
; for i 0 42 3 [
; print [
; i - .01 UTM-zone-number 80.0 i - .01
; tab
; i attempt [UTM-zone-number 80.0 i]
; tab
; i + .01 attempt [UTM-zone-number 80.0 i + .01]
; ]
; ]
UTM-zone: func [
{Return the UTM zone designator for the given coordinate.}
lat [number!]
lon [number!]
][
join pad-num UTM-lon-zone lat lon 2 UTM-lat-zone lat
]
; buff: copy ""
; for lat -80 88 8 [
; for lon -180 180 6 [
; repend buff reform [lat lon attempt [UTM-zone lat lon] newline]
; print [
; lat lon attempt [UTM-zone lat lon]
; ]
; ]
; ]
; cc buff
; buff: copy ""
; line: copy ""
; for lon -180 180 6 [
; append line rejoin [lon ", "]
; ]
; append buff join line newline
; for lat -80 88 8 [
; clear line
; insert line rejoin [lat ", "]
; for lon -180 180 6 [
; ;append repend line reform [lat lon attempt [UTM-zone lat lon]] ", "
; append repend line reform [attempt [UTM-zone lat lon]] ", "
; ]
; ; insert lines to the output matches how a map shows things
; insert buff join line newline
; print line
; ]
; cc buff
;-------------------------------------------------------------------------------
; http://ereimer.net/programs/LatLong-UTM.cpp
; http://robotics.ai.uiuc.edu/~hyoon24/LatLongUTMconversion.py
; void LLtoUTM(int eId, double Lat, double Long, double& Northing, double& Easting, int& Zone){
; // converts LatLong to UTM coords; 3/22/95: by ChuckGantz chuck.gantz@globalstar.com, from USGS Bulletin 1532.
; // Lat and Long are in degrees; North latitudes and East Longitudes are positive.
; double a = ellip[eId].EquatorialRadius;
; double ee= ellip[eId].eccSquared;
; Long -= int((Long+180)/360)*360; //ensure longitude within -180.00..179.9
; double N, T, C, A, M;
; double LatRad = Lat*deg2rad;
; double LongRad = Long*deg2rad;
;
; Zone = int((Long + 186)/6);
; if( Lat >= 56.0 && Lat < 64.0 && Long >= 3.0 && Long < 12.0 ) Zone = 32;
; if( Lat >= 72.0 && Lat < 84.0 ){ //Special zones for Svalbard
; if( Long >= 0.0 && Long < 9.0 ) Zone = 31;
; else if( Long >= 9.0 && Long < 21.0 ) Zone = 33;
; else if( Long >= 21.0 && Long < 33.0 ) Zone = 35;
; else if( Long >= 33.0 && Long < 42.0 ) Zone = 37;
; }
; double LongOrigin = Zone*6 - 183; //origin in middle of zone
; double LongOriginRad = LongOrigin * deg2rad;
;
; double EE = ee/(1-ee);
;
; N = a/sqrt(1-ee*sin(LatRad)*sin(LatRad));
; T = tan(LatRad)*tan(LatRad);
; C = EE*cos(LatRad)*cos(LatRad);
; A = cos(LatRad)*(LongRad-LongOriginRad);
;
; M= a*((1 - ee/4 - 3*ee*ee/64 - 5*ee*ee*ee/256 ) *LatRad
; - (3*ee/8 + 3*ee*ee/32 + 45*ee*ee*ee/1024) *sin(2*LatRad)
; + (15*ee*ee/256 + 45*ee*ee*ee/1024 ) *sin(4*LatRad)
; - (35*ee*ee*ee/3072 ) *sin(6*LatRad));
;
; Easting = k0*N*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*EE)*A*A*A*A*A/120) + 500000.0;
;
; Northing = k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
; + (61-58*T+T*T+600*C-330*EE)*A*A*A*A*A*A/720));
; }
gps-pos-ctx: context [
UTMNorthing: none
UTMEasting: none
UTMzone: copy ""
]
; LLtoUTM: func [
; ReferenceEllipsoid
; Lat [decimal!]
; Long [decimal!]
; /local gps-pos UTMEasting UTMNorthing UTMZone
;
; ] [
; ; //converts lat/long to UTM coords. Equations from USGS Bulletin 1532
; ; //East Longitudes are positive, West longitudes are negative.
; ; //North latitudes are positive, South latitudes are negative
; ; //Lat and Long are in decimal degrees
; ; //Written by Chuck Gantz- chuck.gantz@globalstar.com
; ;variable ReferenceEllipsoid, Lat, Long
; GPSpos: make gps-pos-ctx []
;
; ;variable UTMEasting, UTMNorthing
; ;string UTMZone
;
; make/o/d/n=(23) EquatorialRadius, SquareOfEccentricity
; make/o/n=23/T EllipsoidName
; EquatorialRadius={6377563, 6378160, 6377397, 6377484, 6378206, 6378249, 6377276, 6378166, 6378150, 6378160, 6378137, 6378200, 6378270, 6378388, 6378245, 6377340, 6377304, 6378155, 6378160, 6378165, 6378145, 6378135, 6378137}
; SquareOfEccentricity={.00667054,.006694542,.006674372,.006674372,.006768658,.006803511,.006637847,.006693422,.006693422,.006694605,.00669438,.006693422,.00672267,.00672267,.006693422,.00667054,.006637847,.006693422,.006694542,.006693422,.006694542,.006694318,.00669438}
; EllipsoidName={"Airy","Australian National","Bessel 1841","Bessel 1841 (Nambia) ","Clarke 1866","Clarke 1880","Everest","Fischer 1960 (Mercury) ","Fischer 1968","GRS 1967", "GRS 1980", "Helmert 1906","Hough", "International", "Krassovsky", "Modified Airy", "Modified Everest", "Modified Fischer 1960", "South American 1969","WGS 60","WGS 66","WGS-72", "WGS-84"}
;
;
; variable a = EquatorialRadius[ReferenceEllipsoid]
; variable eccSquared = SquareOfEccentricity[ReferenceEllipsoid];
; variable k0 = 0.9996;
;
; variable LongOrigin;
; variable eccPrimeSquared;
; variable N, T, C, AA, M;
;
; ;//Make sure the longitude is between -180.00 .. 179.9
; variable LongTemp = (Long+180)-floor((Long+180)/360)*360-180; // -180.00 .. 179.9;
;
; variable LatRad = Lat*DEG2RAD;
; variable LongRad = LongTemp*DEG2RAD;
; variable LongOriginRad;
; variable ZoneNumber;
;
; ZoneNumber = floor((LongTemp + 180)/6) + 1;
;
; if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
; ZoneNumber = 32;
; endif
; ; // Special zones for Svalbard
; if( Lat >= 72.0 && Lat < 84.0 )
; if ( LongTemp >= 0.0 && LongTemp < 9.0 )
; ZoneNumber = 31;
; elseif( LongTemp >= 9.0 && LongTemp < 21.0 )
; ZoneNumber = 33;
; elseif(LongTemp >= 21.0 && LongTemp < 33.0 )
; ZoneNumber = 35;
; elseif(LongTemp >= 33.0 && LongTemp < 42.0 )
; ZoneNumber = 37;
; endif
; endif
; LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
; LongOriginRad = LongOrigin * DEG2RAD;
;
; ; //compute the UTM Zone from the latitude and longitude
;
; eccPrimeSquared = (eccSquared)/(1-eccSquared);
;
; N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
; T = tan(LatRad)*tan(LatRad);
; C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
; AA = cos(LatRad)*(LongRad-LongOriginRad);
;
; M = (1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
; M -= (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
; M += (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
; M -= (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad)
; M *= a
;
; UTMEasting = (k0*N*(AA+(1-T+C)*AA*AA*AA/6+ (5-18*T+T*T+72*C-58*eccPrimeSquared)*AA*AA*AA*AA*AA/120)+ 500000.0);
;
; UTMNorthing = (k0*(M+N*tan(LatRad)*(AA*AA/2+(5-T+9*C+4*C*C)*AA*AA*AA*AA/24+ (61-58*T+T*T+600*C-330*eccPrimeSquared)*AA*AA*AA*AA*AA*AA/720)));
; if(Lat < 0)
; UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
; endif
;
; gps.UTMEasting = UTMEasting
; gps.UTMNorthing = UTMNorthing
; gps.UTMZone = num2str(ZoneNumber)+ UTMLetterDesignator(Lat)
; print gps
; ]
; USNG
;
; Grid Zone Designation. The U.S. geographic area is divided into 6-degree
; longitudinal zones designated by a number and 8-degree latitudinal bands
; designated by a letter. Each area is given a unique alpha-numeric Grid Zone
; Designator (GZD) (i.e. 18S).
;
; 100,000-meter Square Identification. Each GZD 6x8 degree area is covered by a
; specific scheme of 100,000-meter squares where each square is identified by two
; unique letters. (i.e. 18SUJ - Identifies a specific 100,000-meter square in the
; specified GZD).
;
; Grid Coordinates A point position within the 100,000-meter square shall be given
; by the UTM grid coordinates in terms of its Easting (E) and Northing (N). An
; equal number of digits shall be used for E and N where the number of digits
; depends on the precision desired in position referencing. In this convention,
; the reading shall be from left with Easting first and then Northing.
;-------------------------------------------------------------------------------
; http://www.russwill.com/library/computers/formulas/Lat.htm
; http://en.wikipedia.org/wiki/Great-circle_distance
; Solve for the great circle distance between two points where the coordinates are
; specified in decimal degrees of latitude and longitude, assuming no change in
; elevation. Longitudes in the western hemisphere must be specified as negative
; values, in the eastern hemisphere, they are positive. Use negative latitudes in
; the southern hemisphere, use positive latitudes the north. Coordinates must be
; entered in decimal degrees of longitude and latitude. Longitude for the western
; hemisphere and latitude for the southern hemisphere are expressed as negative
; values.
Lat-Lon-Distance: func [
lat-1 [decimal!]
lon-1 [decimal!]
lat-2 [decimal!]
lon-2 [decimal!]
units [string! word!] "MI, FT, YD, KM, or M (M = Meters)"
/local
earth-radius
delta-lat delta-lon
central-angle great-circle-dist
][
; Get the radius of the earth in the specified units
earth-radius: switch/default uppercase form units [
"MI" [3'956] ; Miles
"FT" [20'887'680] ; Feet
"YD" [6'962'560] ; Yards
"KM" [6'367] ; Kilometres
"M" [6'367'000] ; Metres 6'378'137
][ ; non-supported units
make error! join [script invalid-arg] units
]
delta-lon: (as-radians lon-2) - (as-radians lon-1)
delta-lat: (as-radians lat-2) - (as-radians lat-1)
; Haversine formula
central-angle: (sine2 (delta-lat / 2)) + ((cosine lat-1) * (cosine lat-2) * (sine2 (delta-lon / 2)))
great-circle-dist: 2 * arcsine min 1 square-root central-angle
earth-radius * great-circle-dist
]
;-------------------------------------------------------------------------------
; http://www.consultsarath.com/contents/articles/KB000012-distance-between-two-points-on-globe--calculation-using-cSharp.aspx
; http://www.consultsarath.com/contents/articles/KB000010-calculate-distance-between-two-points-on-globe-from-latitude-and-longitude-coordinates.aspx
; Haversine formula
GetDistanceBetweenPoints: func [
lat1 [decimal!]
long1 [decimal!]
lat2 [decimal!]
long2 [decimal!]
/local
distance dLat dLong a c radiusE radiusP nr dr radius
] [
distance: 0;
dLat: (lat2 - lat1) / 180 * PI
dLong: (long2 - long1) / 180 * PI
a: sine (dLat / 2) * sine (dLat / 2) + cosine (lat2) * sine (dLong/2) * sine (dLong/2)
c: 2 * atan2 square-root a square-root (1 - a)
;Calculate radius of earth
; For this you can assume any of the two points.
radiusE: 6378135 ; Equatorial radius, in metres
radiusP: 6356750 ; Polar Radius
;Numerator part of function
nr: power (radiusE * radiusP * Cosine (lat1 / 180 * PI)) 2
;Denominator part of the function
dr: power ((radiusE * Cosine (lat1 / 180 * PI)) 2) + power (radiusP * Sine (lat1 / 180 * PI)) 2
radius: square-root (nr / dr)
;Calaculate distance in metres.
distance = radius * c;
return distance;
]
;-------------------------------------------------------------------------------
comment {
The radius of curvature of an ellipsoidal Earth in the plane of the meridian is given by
R' = a * (1 - e^2) / (1 - e^2 * (sin(lat))^2)^(3/2)
where a is the equatorial radius, b is the polar radius, and e is the eccentricity of the ellipsoid = sqrt(1 - b^2/a^2)
and the radius of curvature in a plane perpendicular to the meridian and perpendicular to a plane tangent to the surface is given by
N = a/sqrt(1-e^2*(sin(lat)^2))
A Swedish book suggests use of the geometric mean of these two radii of curvature for all azimuths, as it produces errors of order of magnitude 0.1% for distances within 500 km (300 mi) at 60 degrees latitude. The formula for that average is no more complicated than either of its components. That is, r = sqrt(R' * N) becomes
r = a * sqrt(1 - e^2) / (1 - e^2 * (sin(lat))^2)
Using these formulas with
a = 6378 km (3963 mi) Equatorial radius (surface to center distance)
b = 6357 km (3950 mi) Polar radius (surface to center distance)
e = 0.081082 Eccentricity
gives the following table of values for the Radii of Curvature:
latitude...........r...................R'..................N..........
00 degrees . 6357 km (3950 mi) . 6336 km (3937 mi) . 6378 km (3963 mi)
15 degrees . 6360 km (3952 mi) . 6340 km (3940 mi) . 6379 km (3964 mi)
30 degrees . 6367 km (3957 mi) . 6352 km (3947 mi) . 6383 km (3966 mi)
45 degrees . 6378 km (3963 mi) . 6367 km (3957 mi) . 6389 km (3970 mi)
60 degrees . 6388 km (3970 mi) . 6383 km (3966 mi) . 6394 km (3973 mi)
75 degrees . 6396 km (3974 mi) . 6395 km (3974 mi) . 6398 km (3975 mi)
90 degrees . 6399 km (3976 mi) . 6399 km (3976 mi) . 6399 km (3976 mi)
}
;-------------------------------------------------------------------------------
; http://www.movable-type.co.uk/scripts/latlong-vincenty.html
;
; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
; /* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2010 */
; /* */
; /* from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the */
; /* Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975 */
; /* http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf */
; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
;
; /**
; * Calculates geodetic distance between two points specified by latitude/longitude using
; * Vincenty inverse formula for ellipsoids
; *
; * @param {Number} lat1, lon1: first point in decimal degrees
; * @param {Number} lat2, lon2: second point in decimal degrees
; * @returns (Number} distance in metres between points
; */
; distVincenty: func [
; lat1 [decimal!]
; lon1 [decimal!]
; lat2 [decimal!]
; lon2 [decimal!]
; /local
; a b f L U1 U2 sinU1 sinU2 cosU1 cosU2 lambda lambdaP iterLimit
; sinSigma sinLambda cosLambda cosSigma sigma sinAlpha cosSqAlpha
; cos2SigmaM C uSq deltaSigma s fwdAz revAz
; ][
; a: 6378137 ; WGS-84 ellipsoid params
; b: 6356752.314245
; f: 1 / 298.257223563
;
; L: as-radians (lon2 - lon1)
; U1: arctangent ((1 - f) * tangent (as-radians lat1))
; U2: arctangent ((1 - f) * tangent (as-radians lat2))
; sinU1: sine U1
; cosU1: cosine U1
; sinU2: sine U2
; cosU2: cosine U2
;
; lambda: L
; lambdaP: iterLimit: 100
;
; do [
; sinLambda: sine lambda
; cosLambda: cosine lambda
;
; sinSigma: square-root (
; ((cosU2 * sinLambda) * (cosU2 * sinLambda)) +
; (((cosU1 * sinU2) - (sinU1 * cosU2 * cosLambda)) * ((cosU1 * sinU2) - (sinU1 * cosU2 * cosLambda)))
; )
;
; if zero? sinSigma [return 0] ; co-incident points
;
; cosSigma: (sinU1 * sinU2) + (cosU1 * cosU2 * cosLambda)
; sigma: arctangent2 sinSigma cosSigma
; sinAlpha: cosU1 * cosU2 * sinLambda / sinSigma
; cosSqAlpha: 1 - sinAlpha * sinAlpha
; cos2SigmaM: cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
;
; if not number? cos2SigmaM [cos2SigmaM = 0] ; equatorial line: cosSqAlpha=0 (§6)
; C: f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
;
; lambdaP: lambda
; lambda: L + (1 - C) * f * sinAlpha *
; (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))
;
; ] while ( abs (lambda - lambdaP) > 1e-12 and decr iterLimit > 0);
;
; if zero? iterLimit [return none] ; formula failed to converge
;
; uSq: cosSqAlpha * ((a * a) - (b * b)) / (b * b);
; A: 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
; B: uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
; deltaSigma: B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)-
; B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma)*(-3 + 4 * cos2SigmaM * cos2SigmaM)));
; s: b * A * (sigma - deltaSigma);
;
; s: round/to s .001 ; round to 1mm precision
; return s
;
; ; note: to return initial/final bearings in addition to distance, use something like:
; fwdAz: arctangent2 cosU2 * sinLambda cosU1 * sinU2 - sinU1 * cosU2 * cosLambda
; revAz: arctangent2 cosU1 * sinLambda -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda
; return { distance: s, initialBearing: fwdAz.toDeg(), finalBearing: revAz.toDeg() }
; ]
;-------------------------------------------------------------------------------
comment {
Usage:
<p>Lat 1: <input name="lat1" value="53 09 02N"> Long 1: <input name="long1" value="001 50 40W"></p>
<p>Lat 2: <input name="lat2" value="52 12 19N"> Long 2: <input name="long2" value="000 08 33W"></p>
<input type="button" value="calculate distance"
onClick="result.value = distVincenty(Geo.parseDMS(lat1.value), Geo.parseDMS(long1.value),
Geo.parseDMS(lat2.value), Geo.parseDMS(long2.value)) + ' m'">
<input name="result" value="">
See the Lat/Long page for the deg/min/sec parsing method Geo.parseDMS()
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment