1:- module(geohashing, [hash_point/3,
2 code_graticule/2,
3 code_graticule/3,
4 minesweeper/3,
5 globalhash/2]).
21:- use_module(weblog(info/stock/crox/croxstock)). 22:- use_module(library(semweb/rdf_db)). 23:- use_module(library(dcg/basics)).
@Code see above @Grat opaque graticule structure */
37code_graticule(Code, graticule(Lat, Long)) :- 38 var(Code), 39 float(Lat), 40 float(Long), 41 Lat >= -90.0, 42 Lat =< 90.0, 43 Long >= -180.0, 44 Long =< 180.0, 45 codes_gnum(CLat, Lat), 46 codes_gnum(CLong, Long), 47 format(codes(Code), '~s, ~s', [CLat, CLong]). 48code_graticule(Code, graticule(Lat, Long)) :- 49 is_list(Code), 50 phrase(grat_syntax(RLat, RLong), Code), 51 wrap_long(RLong, Long), 52 bound_lat(RLat, Lat), 53 Lat >= -90.0, 54 Lat =< 90.0, 55 Long >= -180.0, 56 Long =< 180.0 .
@Code see above @Grat opaque graticule structure */
71code_graticule(CodeLat, CodeLong, graticule(RLat, RLong)) :- 72 var(CodeLat), 73 bound_lat(RLat, Lat), 74 wrap_long(RLong, Long), 75 codes_gnum(CodeLat, Lat), 76 codes_gnum(CodeLong, Long). 77code_graticule(CodeLat, CodeLong, graticule(Lat, Long)) :- 78 is_list(CodeLat), 79 is_list(CodeLong), 80 phrase(half_grat(RLat), CodeLat), 81 phrase(half_grat(RLong), CodeLong), 82 bound_lat(RLat, Lat), 83 wrap_long(RLong, Long). 84 85wrap_long(RLong, Long) :- 86 RLong >= -180.0, 87 wrap_long_upper(RLong, Long). 88wrap_long(RLong, Long) :- 89 RLong < -180.0, 90 NLong is RLong + 360.0, 91 wrap_long(NLong, Long). 92wrap_long_upper(RLong, RLong) :- 93 RLong < 180.0. 94wrap_long_upper(RLong, Long) :- 95 RLong >= 180.0, 96 NLong is RLong - 360.0, 97 wrap_long_upper(NLong, Long). 98bound_lat(RLat, -90.0) :- 99 RLat < -90.0,! . 100bound_lat(RLat, 90.0) :- 101 RLat > 90.0,! . 102bound_lat(Lat, Lat).
108codes_gnum("-0", -0.0) :- !. 109codes_gnum(CG, G) :- 110 I is truncate(G), 111 format(codes(CG), '~d', [I]). 112 113grat_syntax(Lat, Long) --> 114 blanks, grat(Lat), frac_part, ",", 115 blanks, grat(Long), frac_part, blanks. 116 117half_grat(Coord) --> 118 blanks, grat(Coord), blanks. 119 120grat(-0.0) --> 121 "-0",!. 122grat(G) --> 123 integer(I), 124 { 125 G is 1.0 * I 126 }. 127 128frac_part --> []. 129frac_part --> ".", digits(_). 130 131 132:- dynamic a_hash_point/3.
point(Lat, Long)
on 2008-05-27 the rules for computation of points E of 30W changed
*/
140hash_point(Date, Grat, Point) :- 141 a_hash_point(Date, Grat, Point),!. 142 143hash_point(Date, Grat, Pt) :- 144 hash_atom( Date, Grat, HashAtom), 145 point_from_atom(HashAtom, Grat, Pt), 146 asserta(a_hash_point(Date, Grat, Pt)). 147 148hash_atom( YY - MM - DD, graticule(_, Long), ToHash) :- 149 LLong is float_integer_part(Long), 150 ( YY < 1000 -> YYYY is YY + 2000 ; YYYY is YY), 151 ( east_rule(YYYY - MM - DD, LLong) -> DDDD is DD - 1; DDDD is DD), 152 normalize_geohash_date(YYYY - MM - DD, ZY - ZM - ZD), 153 ( ZM < 10 -> MFill = '0' ; MFill = ''), 154 ( ZD < 10 -> DFill = '0' ; DFill = ''), 155 normalize_geohash_date(YYYY - MM - DDDD, StockY - StockM - StockD), 156 crstock_stats(djia, opening, StockY - StockM - StockD, Stock), 157 format(atom(ToHash), '~d-~w~d-~w~d-~2f', 158 [ZY, MFill, ZM, DFill, ZD, Stock]). 159 160point_from_atom(ToHash, graticule(Lat, Long), point(PLat, PLong)) :- 161 LLat is float_integer_part(Lat), 162 LLong is float_integer_part(Long), 163 % The only available MD5 *DIGEST* implementation is from semweb 164 rdf_atom_md5(ToHash, 1, HashedAtom), 165 atom_codes(HashedAtom, HashedCodes), 166 append(FirstHalf, LastHalf, HashedCodes), 167 same_length(FirstHalf, LastHalf),!, 168 hex_codes_frac_float(FirstHalf, 1.0, LatFrac), 169 hex_codes_frac_float(LastHalf, 1.0, LongFrac), 170 grat_geo(LLat , LatFrac, PLat), 171 grat_geo(LLong, LongFrac, PLong).
179 globalhash(Date, point(GLat, GLong)) :- 180 hash_point(Date, graticule(10, -35), point(PLat, PLong)), 181 GLat is abs(float_fractional_part(PLat)) * 180.0 - 90.0, 182 GLong is abs(float_fractional_part(PLong)) * 360.0 - 180.0.
191grat_geo(-0.0, Frac, Geo) :- 192 !, 193 Geo is -1.0 * Frac. 194grat_geo(Num, Frac, Geo) :- 195 Num >= 0.0, 196 Geo is Num + Frac. 197grat_geo(Num, Frac, Geo) :- 198 Num < 0.0, 199 Geo is Num - Frac.
point(Lat, Long)
structures for the 9 nearest hash points
note - at the poles returns a shorter list
*/
212minesweeper(Date, Grat, Points) :- 213 setof(Pt, a_ms_pt(Date, Grat, Pt), Points). 214 215a_ms_pt(Date, Grat, Pt) :- 216 member(OffLat, [-1, 0, 1]), 217 member(OffLong, [-1, 0, 1]), 218 hash_offset(Grat, OffLat, OffLong, NGrat), 219 hash_point(Date, NGrat, Pt).
225hash_offset(graticule(Lat, Long), OffLat, OffLong, graticule(NLat, NNLong)) :- 226 coord_offset(Lat, OffLat, NLat), 227 NLat >= -90.0, 228 NLat =< 90.0, 229 coord_offset(Long, OffLong, NLong), 230 ( NLong =< -180.0 -> 231 NNLong is NLong + 360.0 ; 232 NNLong is NLong 233 ), 234 ( NLong >= 180.0 -> 235 NNLong is NLong - 360.0 ; 236 NNLong is NLong 237 ). 238 239coord_offset(G, 0, G). 240coord_offset(-0.0, 1, 0.0). 241coord_offset(0.0, -1, -0.0). 242coord_offset(G, N, NG) :- NG is N + G. 243 244normalize_geohash_date(Y - M - D, YO - MO - DO) :- 245 date_time_stamp(date(Y, M, D, 0,0,0,0,-,-), Stamp), 246 stamp_date_time(Stamp, Dt, 0), 247 date_time_value(date, Dt, date(YO, MO, DO)).
254east_rule(YYYY - _ - _, IntLong) :- 255 IntLong > -30, 256 YYYY > 2008. 257east_rule(YYYY - MM - _, IntLong) :- 258 IntLong > -30, 259 YYYY = 2008, 260 MM > 5. 261east_rule(YYYY - MM - DD, IntLong) :- 262 IntLong > -30, 263 YYYY = 2008, 264 MM = 5, 265 DD >= 27.
Codes is expected to be a list of ascii codes of hex digits, Place is the fractional place to the extreme left, Value is the floating pt equivilent
*/
277hex_codes_frac_float([], _, 0.0). 278hex_codes_frac_float([H|T], Place, Value) :- 279 char_type(H, xdigit(D)), 280 NewPlace is Place / 16.0, 281 hex_codes_frac_float(T, NewPlace, RemVal), 282 Value is NewPlace * D + RemVal. 283 284 285/* Test suite below here */ 286 287% main entry point for testing 288test_geohashing :- 289 retractall(croxstock:stock_stat(_,_,_,_)), 290 retractall(a_hash_point(_,_,_)),!, 291 test_hash_end_end,!, 292 test_east_rule,!, 293 test_grat_geo,!, 294 test_hashatom,!, 295 test_hex_frac,!, 296 format('finished~n'). 297 298% test conversion from hex string fraction to float fraction 299test_hex_frac :- 300 test_case_hex_frac(Hex, Frac), 301 call(hex_codes_frac_float(Hex, 1.0, GotFrac)), 302 ( abs(GotFrac - Frac) < 0.0001 -> true 303 ; format('hex_codes ~s expected ~w got ~w~n', [Hex, Frac, GotFrac]) 304 ), 305 fail. 306test_hex_frac. 307 308test_case_hex_frac("db9318c2259923d0", 0.857713). 309test_case_hex_frac("8b672cb305440f97", 0.544543). 310test_case_hex_frac("ffffffffffffffff", 0.999999). 311test_case_hex_frac("0000000000000000", 0.0). 312test_case_hex_frac("8000000000000000", 0.5). 313 314% test converting graticule and fractional location to geo coord 315test_grat_geo :- 316 test_case_gg(GratNum, Frac, Geo), 317 call(grat_geo(GratNum, Frac, GotGeo)), 318 ( GotGeo =:= Geo -> true 319 ; format('grat_geo case ~w should be ~w~n', [grat_geo(GratNum, Frac, GotGeo), Geo]) 320 ), 321 fail. 322test_grat_geo. 323 324test_case_gg(40, 0.56, 40.56). 325test_case_gg(40.0, 0.87, 40.87). 326test_case_gg(0.0 , 0.66, 0.66). 327test_case_gg(-0.0, 0.66, -0.66). 328test_case_gg(-2.0, 0.45, -2.45). 329 330 331% test the east rule rule 332test_east_rule :- 333 test_east_data(Dt, Long, true), 334 east_rule(Dt, Long), 335 ( east_rule(Dt, Long) -> 336 true 337 ; 338 format('eastrule should apply to ~w ~w but doesnt~n', [Dt, Long]) 339 ), 340 fail. 341test_east_rule :- 342 test_east_data(Dt, Long, false), 343 ( east_rule(Dt, Long) -> 344 format('eastrule should not apply to ~w ~w but does~n', [Dt, Long]) 345 ), 346 fail. 347test_east_rule. 348 349test_east_data( 2008 - 5 - 20, 0, false). 350test_east_data( 2013 - 5 - 28, -31, false). 351test_east_data( 2013 - 5 - 28, 0, true). 352 353 354% Test the entire hashpoint generation 355test_hash_end_end :- 356 test_item(Y-M-D, GLat, GLong, Lat, Long), 357 hash_point(Y-M-D, graticule(GLat, GLong), point(PLat, PLong)), 358 Dist is sqrt((Lat - PLat) * (Lat - PLat) + 359 (Long - PLong) * (Long - PLong)), 360 ( Dist > 0.0001 -> % 36 ft at equator 361 format('~w-~w-~w_~w_~w should be ~w,~w but is ~w,~w~n', 362 [Y,M,D,GLat,GLong,Lat,Long,PLat,PLong]) 363 ; true), 364 fail. 365test_hash_end_end. 366 367test_item(2013 - 5 - 31, 50, 8, 50.47615093, 8.75973636). 368test_item(2013 - 5 - 31, 52, 1, 52.47615093, 1.75973636). 369test_item(2013 - 5 - 31, 52, -0.0, 52.47615093, -0.75973636). 370test_item(2013 - 5 - 31, 52, 0.0, 52.47615093, 0.75973636). 371test_item(2013 - 5 - 31, 52, -1.0, 52.47615093, -1.75973636). 372test_item(2013 - 6 - 6, 37, -120, 37.64757933, -120.17135344). 373test_item(2013 - 6 - 6, 52, -2, 52.13259239, -2.18860213). 374test_item(2013 - 6 - 6, 52, 2, 52.13259239, 2.18860213). 375 376% The 30W compliance test 377test_item(2008 - 5 - 20, 68, -30, 68.63099, -30.61895 ). 378test_item(2008 - 5 - 21, 68, -30, 68.17947, -30.86154 ). 379test_item(2008 - 5 - 22, 68, -30, 68.97287, -30.2387 ). 380test_item(2008 - 5 - 23, 68, -30, 68.40025, -30.72277 ). 381test_item(2008 - 5 - 24, 68, -30, 68.12665, -30.54753 ). 382test_item(2008 - 5 - 25, 68, -30, 68.94177, -30.18287 ). 383test_item(2008 - 5 - 26, 68, -30, 68.67313, -30.60731 ). 384test_item(2008 - 5 - 27, 68, -30, 68.20968, -30.10144 ). 385test_item(2008 - 5 - 28, 68, -30, 68.68745, -30.21221 ). 386test_item(2008 - 5 - 29, 68, -30, 68.4647, -30.03412 ). 387test_item(2008 - 5 - 30, 68, -30, 68.8531, -30.2446 ). 388test_item(2008 - 5 - 20, 68, -29, 68.63099, -29.61895 ). 389test_item(2008 - 5 - 21, 68, -29, 68.17947, -29.86154 ). 390test_item(2008 - 5 - 22, 68, -29, 68.97287, -29.2387 ). 391test_item(2008 - 5 - 23, 68, -29, 68.40025, -29.72277 ). 392test_item(2008 - 5 - 24, 68, -29, 68.12665, -29.54753 ). 393test_item(2008 - 5 - 25, 68, -29, 68.94177, -29.18287 ). 394test_item(2008 - 5 - 26, 68, -29, 68.67313, -29.60731 ). 395test_item(2008 - 5 - 27, 68, -29, 68.12537, -29.57711 ). 396test_item(2008 - 5 - 28, 68, -29, 68.71044, -29.11273 ). 397test_item(2008 - 5 - 29, 68, -29, 68.27833, -29.74114 ). 398test_item(2008 - 5 - 30, 68, -29, 68.32272, -29.70458 ). 399 400 401test_hashatom :- 402 hashatom_test_case(Dt, Grat, HashAtom), 403 call(hash_atom(Dt, Grat, GotAtom)), 404 ( GotAtom = HashAtom -> true ; 405 format('hash_atom ~w ~w expeected ~q got ~q~n', [Dt, Grat, HashAtom, GotAtom]) 406 ), 407 fail. 408test_hashatom. 409 410% East rule still uses current date for string! 411hashatom_test_case(2013 - 5 - 31, graticule(50, 8), '2013-05-31-15306.02'). % east rule 412hashatom_test_case(2013 - 6 - 6 , graticule(37, -120), '2013-06-06-14955.45'). 413 414% peeron map has small 'toggle debug info' in LL corner. Very useful. 415% 416% some stock prices 417% 2013-05-31 15322.22 418% 2013-05-30 15306.02 419% 2013-06-06 14955.45
Tools for geohashing
Tools for the sport of geohashing.
http://wiki.xkcd.com/geohashing/Main_Page
Because there are two zeros in graticules (if the fractional part of the hash longitude is 0.35, then the hash for +0 graticule is 0.35w and 0.35e for the -0 graticule. Thus So there's a bit of type here, with a graticule being a structure that handles all this.
*/