1:- module(geohashing, [hash_point/3,
    2		       code_graticule/2,
    3		       code_graticule/3,
    4		       minesweeper/3,
    5		      globalhash/2]).

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.

*/

   21:- use_module(weblog(info/stock/crox/croxstock)).   22:- use_module(library(semweb/rdf_db)).   23:- use_module(library(dcg/basics)).
 code_graticule(?Code:codes, ?Grat:graticule) is semidet
Attempts to convert codes to a graticule structure. Codes should look something like "-14,3" (the 14 deg N latitude, 3 deg E longitude graticule) "52, 0" (Cambridge, UK graticule) "52, -0" (Northampton, UK graticule)

@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_graticule(?CodeLat:codes, ?CodeLong:codes, ?Grat:graticule) is semidet
Attempts to convert a pair of codes to a graticule structure. Codes should look something like "-14" (the 14 deg coord) "0" (the 0 coord) " -0" (the -0 coord)

@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).
 codes_gnum(-C:codes, +G:float) is det
convert a graticule float value to codes */
  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.
 hash_point(+Date:term, +Lat:number, +Long:number, -Point:term) is det
given a date and lat/long, returns 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).
 globalhash(+Date:date, -Point:point) is semidet
find location of globalhash for this date */
  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.
 grat_geo(+GratNum, +Frac, -Geo) is det
get a real geocoord from a graticule number and fraction */
  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.
 minesweeper(+Date:term, +Grat:graticule, -Points:list) is det
given a date and graticule, returns a list of nine 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).
 hash_offset(+G:graticule, +OffLat:integer, +OffLong:integer, -NG:graticule) is semidet
offset a graticule. fails if new graticule isn't on earth. */
  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)).
 east_rule(+Date:geohashdate, +IntLong:int) is semidet
unifies if this location and date invokes the 'east of 30W' rule */
  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.
 hex_codes_frac_float(+Codes:codes, +Place:float, -Value:float) is semidet
converts a hex string to a fractional float value

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