source: trunk/npemap.org.uk/scripts/gmc-london/find-best-match.py @ 537

Last change on this file since 537 was 537, checked in by Nick Burch, 13 years ago

For working with gmc maps

  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 4.0 KB
Line 
1#!/usr/bin/python
2import math
3
4places = {
5        'Hither Green Sta':   ( (539016,174424), (22751,1390) ),
6        'Stratford Broadway': ( (539000,184488), (22898,13708) ),
7        'Mill Hill Road':     ( (522000,176056), (2105,3856) ),
8        'Park Drive':         ( (520896,174936), (646,2487) ),
9        'Lewisham High St':   ( (538152,175104), (21708,2254) ),
10        'Hither Green':       ( (538152,174928), (21713,2050) ),
11        'SW Corner Church':   ( (522584,173712), (2725,925) ),
12        'Charlton Sta':       ( (541116,178384), (25442,6112) ),
13        'Wanstead Park St':   ( (540560,185620), (24920,15035) ),
14        'Blackbird Hill Jnct':( (520328,186628), (349,16845) ),
15        'Pond':               ( (522000,185944), (2383,15898) ),
16        'QueenVicMemorial':   ( (529160,179760), (10866,8193) ),
17        'AlbertMemorial':     ( (526592,179744), (7729,8263) ),
18        'RegentsParkTube':    ( (528680,182112), (10362,5008) ),
19        'EustonRdCorner':     ( (529752,182592), (11697,5595) ),
20        'TowerOfLondon':      ( (533648,180544), (16268,7075) ),
21        'CanningTownSta':     ( (539304,181544), (23239,6190) ),
22}
23
24# Track the errors
25errors = {}
26
27def calculate_slant(slant_base,slant_adj_e,slant_adj_n,gmc_e,gmc_n):
28        # The further east we go, the smaller the slant is
29        slant = slant_base + \
30                                slant_adj_e*((26000.0-gmc_e)/26000.0) + \
31                                slant_adj_n*((17000.0-gmc_n)/17000.0)
32        slant_rad = slant / 360 * 2 * math.pi
33        return slant_rad
34
35# To do a convertion from gmc e+n to offset e+n
36def convert_to_offset(slant_rad,scale,gmc_e,gmc_n):
37        # The grid squares are supposed to be 0.5 miles across
38        # In actual fact, they are a tiny bit too small
39        adj_gmc_e = gmc_e * scale
40        adj_gmc_n = gmc_n * scale
41
42        # Turn them into meters from the start of the map
43        m_e = adj_gmc_e / 2 * 1.6
44        m_n = adj_gmc_n / 2 * 1.6
45
46        # Apply the slant factor
47        real_e = (math.sin(slant_rad) * m_n) + (math.cos(slant_rad) * m_e)
48        real_n = (math.cos(slant_rad) * m_n) + (math.sin(slant_rad) * m_e)
49
50        return (real_e,real_n)
51
52
53# We can vary the slant (base and easterly dependence)
54# We can also vary the scale factor
55for slant_base_int in range(35,141):
56        slant_base = slant_base_int / 100.0
57        print "Slant base is %0.2f, adjust is -0.50 to 1.20" % slant_base
58
59        for slant_adj_e_int in range(-50,121):
60                slant_adj_e = slant_adj_e_int / 100.0
61                print "  Trying slant easting dependence of %0.2f" % slant_adj_e
62
63                for slant_adj_n_int in range(-50,121):
64                        slant_adj_n = slant_adj_n_int / 100.0
65                        #print "    Trying slant northing dependence of %0.2f" % slant_adj_n
66
67                        # Get the slants
68                        slants = {}
69                        for place,values in places.items():
70                                slants[place] = calculate_slant(slant_base,slant_adj_e,slant_adj_n,values[1][0],values[1][1])
71                        slant_00 = calculate_slant(slant_base,slant_adj_e,slant_adj_n,646,2487)
72
73                        for scale_int in range(950,1251):
74                                scale = scale_int / 1000.0
75                                error = 0
76
77                                # Calculate the left corner e+n
78                                # (We know that 0,646 2,487 is 520,896 174,936)
79                                ngr_offset = convert_to_offset(slant_00,scale,646,2487)
80                                ngr_e_0 = 520896 - ngr_offset[0]
81                                ngr_n_0 = 174936 - ngr_offset[1]
82
83                                # Convert each one, and figure out the error
84                                for place,values in places.items():
85                                        ngr_values = values[0]
86                                        gmc_values = values[1]
87
88                                        base_pos = convert_to_offset(slants[place],scale,gmc_values[0],gmc_values[1])
89                                        found_ngr = (base_pos[0] + ngr_e_0, base_pos[1] + ngr_n_0)
90
91                                        #print "Should be %d %d - was %d %d" % (ngr_values[0],ngr_values[1], found_ngr[0], found_ngr[1])
92                                        this_error = math.sqrt( 
93                                                                                ((ngr_values[0]-found_ngr[0]) ** 2) +
94                                                                                ((ngr_values[1]-found_ngr[1]) ** 2)
95                                        )
96                                        error += this_error
97                                # Save the values and their error, if it's low enough
98                                #print "  Adj are %0.2f,%0.2f, scale is %0.3f, error is %d" % (slant_adj_e,slant_adj_n,scale,error)
99                                if error < 5000:
100                                        errors[(slant_base,slant_adj_e,slant_adj_n,scale)] = error
101
102# What had the lowest error?
103values = errors.keys()
104print "Ordering by error:"
105values.sort(lambda a,b: cmp(errors[a],errors[b]))
106
107print ""
108for this_val in values[0:15]:
109        print "%d (%d): %0.2f , %0.2f, %0.3f" % (errors[this_val], errors[this_val] / len(places.keys()), this_val[0], this_val[1], this_val[2])
Note: See TracBrowser for help on using the repository browser.