Changeset 508
 Timestamp:
 Mar 25, 2007, 6:17:41 PM (15 years ago)
 Location:
 trunk/npemap.org.uk/scripts/freethepostcode.orgimporter
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/npemap.org.uk/scripts/freethepostcode.orgimporter/geo_helper.py
r33 r508 5 5 # XYZ > Lat+Long+Height 6 6 # Lat+Long > other Lat+Long (Helmert Transform) 7 # Lat+Long > easting/northing (OS Only)8 # easting/northing > Lat+Long (OS Only)7 # Lat+Long > easting/northing (OS GB+IE Only) 8 # easting/northing > Lat+Long (OS GB+IE Only) 9 9 # OS easting/northing > OS 6 figure ref 10 10 # … … 13 13 # GPL 14 14 # 15 # Nick Burch  v0.0 3 (25/06/2006)15 # Nick Burch  v0.05 (25/03/2007) 16 16 17 17 import math … … 80 80 } 81 81 82 ############################################################## 83 # OS Specific Helpers for Generic Methods # 82 # Cassini Projection Origins 83 # List is lat (rad), long (rad), false easting, false northing 84 cassini_values = { 85 'osgb' : [ (53.0 + (13.0 / 60.0) + (17.274 / 3600.0)) /360.0 *2.0*math.pi, 86 (2.0 + (41.0 / 60.0) + (3.562 / 3600.0)) /360.0 *2.0*math.pi, 87 0, 0 ] 88 } 89 90 # How many feet to the meter 91 feet_per_meter = 1.0 / 0.3048007491 # 3.28083 92 93 ############################################################## 94 # OS GB Specific Helpers for Generic Methods # 84 95 ############################################################## 85 96 … … 94 105 osgb_latlong = turn_xyz_into_llh( 95 106 osgb_xyz[0],osgb_xyz[1],osgb_xyz[2],'osgb') 96 97 107 return osgb_latlong 108 98 109 def turn_osgb36_into_wgs84(lat_dec,long_dec,height): 99 110 """See http://www.gps.gov.uk/guide6.asp#6.2 and http://www.gps.gov.uk/guide6.asp#6.6 for the calculations, and http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34h.html for some background.""" … … 116 127 """Turn OSGB36 easting and northing values into (decimal) lat/long values inOSGB36.""" 117 128 return turn_eastingnorthing_into_latlong(easting,northing,'osgb') 129 130 ############################################################## 131 # OS IE Specific Helpers for Generic Methods # 132 ############################################################## 133 134 def turn_wgs84_into_osie36(lat_dec,long_dec,height): 135 """As per turn_wgs84_into_osgb36, but for Irish grid""" 136 137 wgs84_xyz = turn_llh_into_xyz(lat_dec,long_dec,height,'wgs84') 138 139 osie_xyz = turn_xyz_into_other_xyz( 140 wgs84_xyz[0],wgs84_xyz[1],wgs84_xyz[2],'wgs84','osie') 141 142 osie_latlong = turn_xyz_into_llh( 143 osie_xyz[0],osie_xyz[1],osie_xyz[2],'osie') 144 return osie_latlong 145 146 def turn_osie36_into_wgs84(lat_dec,long_dec,height): 147 """As per turn_osgb36_into_wgs84, but for Irish grid""" 148 149 osie_xyz = turn_llh_into_xyz(lat_dec,long_dec,height,'osie') 150 151 wgs84_xyz = turn_xyz_into_other_xyz( 152 osie_xyz[0],osie_xyz[1],osie_xyz[2],'osie','wgs84') 153 154 wgs84_latlong = turn_xyz_into_llh( 155 wgs84_xyz[0],wgs84_xyz[1],wgs84_xyz[2],'wgs84') 156 157 return wgs84_latlong 158 159 def turn_osie36_into_eastingnorthing(lat_dec,long_dec): 160 """Turn OSIE36 (decimal) lat/long values into OS IE easting and northing values.""" 161 return turn_latlong_into_eastingnorthing(lat_dec,long_dec,'osie') 162 163 def turn_eastingnorthing_into_osie36(easting,northing): 164 """Turn OSIE36 easting and northing values into (decimal) lat/long values inOSIE36.""" 165 return turn_eastingnorthing_into_latlong(easting,northing,'osie') 118 166 119 167 ############################################################## … … 323 371 324 372 return (lat,long) 373 374 ############################################################## 375 # Cassini Easting/Northing Transform Methods # 376 ############################################################## 377 378 def turn_latlong_into_cassini_en(lat_dec,long_dec,scheme): 379 """Latitude and Longitude, into CassiniSoldner easting and northing coordinates, in the given scheme. See http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34g.html for details of the calculation used""" 380 381 a = abe_values[scheme][0] 382 b = abe_values[scheme][1] 383 e2 = abe_values[scheme][2] 384 385 e4 = e2 * e2 386 e6 = e2 * e2 * e2 387 388 theta = float(lat_dec) /360.0 *2.0*math.pi 389 landa = float(long_dec) /360.0 *2.0*math.pi 390 391 theta0 = cassini_values[scheme][0] 392 landa0 = cassini_values[scheme][1] 393 false_easting = cassini_values[scheme][2] 394 false_northing = cassini_values[scheme][3] 395 396 # Compute intermediate values 397 A = (landa  landa0) * math.cos(theta) 398 T = math.tan(theta) * math.tan(theta) 399 C = e2 / (1.0  e2) * math.cos(theta) * math.cos(theta) 400 v = a / math.sqrt( 1  (e2 * math.sin(theta) * math.sin(theta)) ) 401 402 A2 = A ** 2 403 A3 = A ** 3 404 A4 = A ** 4 405 A5 = A ** 5 406 407 # And M, which is how far along the meridian our latitude is from the origin 408 def makeM(picked_theta): 409 return a * ( 410 (1.0  e2/4.0  3.0*e4/64.0  5.0*e6/256.0) * picked_theta 411  (3.0*e2/8.0 + 3.0*e4/32.0 + 45.0*e6/1024.0) * math.sin(2.0*picked_theta) 412 + (15.0*e4/256.0 + 45.0*e6/1024.0) * math.sin(4.0*picked_theta) 413  (35.0*e6/3072.0) * math.sin(6.0*picked_theta) 414 ) 415 M = makeM(theta) 416 M0 = makeM(theta0) 417 418 # Now calculate 419 easting = false_easting + v * ( 420 A  T * A3 / 6.0  (8.0  T + 8.0*C) * T * A5 / 120.0 ) 421 northing = false_northing + M  M0 + v * math.tan(theta) * ( 422 A2 / 2.0 + (5.0  T + 6.0*C) * A4 / 24.0 ) 423 424 return (easting,northing) 425 426 def turn_cassini_en_into_latlong(easting,northing,scheme): 427 """CassiniSoldner easting and northing, into Latitude and Longitude, in the given scheme. See http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34g.html for details of the calculation used""" 428 429 a = abe_values[scheme][0] 430 b = abe_values[scheme][1] 431 e2 = abe_values[scheme][2] 432 433 e4 = e2 * e2 434 e6 = e2 * e2 * e2 435 436 theta0 = cassini_values[scheme][0] 437 landa0 = cassini_values[scheme][1] 438 false_easting = cassini_values[scheme][2] 439 false_northing = cassini_values[scheme][3] 440 441 def makeM(picked_theta): 442 return a * ( 443 (1.0  e2/4.0  3.0*e4/64.0  5.0*e6/256.0) * picked_theta 444  (3.0*e2/8.0 + 3.0*e4/32.0 + 45.0*e6/1024.0) * math.sin(2.0*picked_theta) 445 + (15.0*e4/256.0 + 45.0*e6/1024.0) * math.sin(4.0*picked_theta) 446  (35.0*e6/3072.0) * math.sin(6.0*picked_theta) 447 ) 448 449 # Compute first batch of intermediate values 450 M1 = makeM(theta0) + (northing  false_northing) 451 mu1 = M1 / (a * (1.0  e2/4.0  3.0*e4/64.0  5.0*e6/256.0) ) 452 e1 = (1  ((1e2) ** 0.5)) / (1 + ((1e2) ** 0.5)) 453 454 e1_2 = e1 ** 2 455 e1_3 = e1 ** 3 456 e1_4 = e1 ** 4 457 458 # Now compute theta1 at T1 459 theta1 = mu1 + ( 460 + (3.0*e1 / 2.0  27.0*e1_3 / 32.0) * math.sin(2.0*mu1) 461 + (21.0*e1_2 / 16.0  55.0*e1_4 / 32.0) * math.sin(4.0*mu1) 462 + (151.0*e1_3 / 96.0) * math.sin(6.0*mu1) 463 + (1097.0*e1_4 / 512.0) * math.sin(8.0*mu1) 464 ) 465 T1 = (math.tan(theta1)) ** 2 466 467 # Now we can find v1, ro1 and D 468 v1 = a / math.sqrt( 1.0  (e2 * math.sin(theta1) * math.sin(theta1)) ) 469 ro1 = a * (1  e2) / ((1  e2 * math.sin(theta1) * math.sin(theta1)) ** 1.5) 470 D = (easting  false_easting) / v1 471 472 # And finally the lat and long 473 lat = theta1  (v1 * math.tan(theta1)) / ro1 * ( 474 D*D/2.0  (1.0 + 3.0 * T1) * ( (D**4) / 24.0 ) ) 475 long = landa0 + ( 476 D  T1 * (D**3) / 3.0 + (1 + 3.0 * T1) * T1 * (D**5) / 15.0 477 ) / math.cos(theta1) 478 479 # Now make decimal versions 480 lat_dec = lat * 360.0 / 2.0 / math.pi 481 long_dec = long * 360.0 / 2.0 / math.pi 482 483 return (lat_dec,long_dec) 325 484 326 485 ############################################################## 
trunk/npemap.org.uk/scripts/freethepostcode.orgimporter/importer.py
r179 r508 10 10 11 11 from pyPgSQL import PgSQL 12 from geo_helper import turn_wgs84_into_osgb36, turn_osgb36_into_eastingnorthing 12 from geo_helper import turn_wgs84_into_osgb36, turn_osgb36_into_eastingnorthing, \ 13 turn_wgs84_into_osie36, turn_osie36_into_eastingnorthing 13 14 import os 14 15 … … 17 18 dbname = "npemaps" 18 19 dbhost = "localhost" 19 dbuser = "npemap" 20 dbpass = "npemap" 21 22 23 raise Exception("Script needs to be updated as per ticket #39") 20 dbuser = "npemaps" 21 dbpass = "" 24 22 25 23 … … 87 85 print "Invalid line '%s'" % line 88 86 89 # Turn lat+long into easting+northing90 osll = turn_wgs84_into_osgb36(parts[0], parts[1], 0)91 en = turn_osgb36_into_eastingnorthing(osll[0], osll[1])92 93 87 pc = {} 94 pc["easting"] = en[0]95 pc["northing"] = en[1]96 88 pc["outer"] = parts[2] 97 89 pc["inner"] = parts[3] … … 99 91 pc["raw_outer"] = parts[2] 100 92 pc["raw_inner"] = parts[3] 93 94 pc["easting"] = None 95 pc["northing"] = None 96 pc["ie_easting"] = None 97 pc["ie_northing"] = None 98 99 # Turn lat+long into easting+northing 100 # All NI postcodes are BT 101 if pc["outer"][0:2] == 'BT': 102 osll = turn_wgs84_into_osie36(parts[0], parts[1], 0) 103 en = turn_osie36_into_eastingnorthing(osll[0], osll[1]) 104 pc["ie_easting"] = en[0] 105 pc["ie_northing"] = en[1] 106 else: 107 osll = turn_wgs84_into_osgb36(parts[0], parts[1], 0) 108 en = turn_osgb36_into_eastingnorthing(osll[0], osll[1]) 109 pc["easting"] = en[0] 110 pc["northing"] = en[1] 101 111 102 112 postcodes.append(pc) … … 138 148 139 149 # Add the latest list to the database 140 sql = "INSERT INTO postcodes (outward, inward, raw_postcode_outward, raw_postcode_inward, easting, northing, source) VALUES (%s, %s, %s, %s, %s, %s, %s)"150 sql = "INSERT INTO postcodes (outward, inward, raw_postcode_outward, raw_postcode_inward, easting, northing, ie_easting, ie_northing, source) VALUES (%s, %s, %s, %s, %s, %s, %s, %s, %s)" 141 151 sth = dbh.cursor() 142 152 worked = 0 153 154 # Add the postcodes 143 155 for postcode in postcodes: 144 sth.execute(sql, (postcode["outer"], postcode["inner"], postcode["raw_outer"], postcode["raw_inner"], postcode["easting"], postcode["northing"], ftp_source))156 sth.execute(sql, (postcode["outer"], postcode["inner"], postcode["raw_outer"], postcode["raw_inner"], postcode["easting"], postcode["northing"], postcode["ie_easting"], postcode["ie_northing"], ftp_source)) 145 157 worked = worked + sth.rowcount 158 146 159 print "Added %d entries" % worked 147 160
Note: See TracChangeset
for help on using the changeset viewer.