Changeset 508 for trunk


Ignore:
Timestamp:
Mar 25, 2007, 6:17:41 PM (13 years ago)
Author:
Nick Burch
Message:

Support IE stuff

Location:
trunk/npemap.org.uk/scripts/freethepostcode.org-importer
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/npemap.org.uk/scripts/freethepostcode.org-importer/geo_helper.py

    r33 r508  
    55#  XYZ -> Lat+Long+Height
    66#  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)
    99#  OS easting/northing -> OS 6 figure ref
    1010#
     
    1313# GPL
    1414#
    15 # Nick Burch - v0.03 (25/06/2006)
     15# Nick Burch - v0.05 (25/03/2007)
    1616
    1717import math
     
    8080}
    8181
    82 ##############################################################
    83 #           OS Specific Helpers for Generic Methods          #
     82# Cassini Projection Origins
     83# List is lat (rad), long (rad), false easting, false northing
     84cassini_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
     91feet_per_meter = 1.0 / 0.3048007491   # 3.28083
     92
     93##############################################################
     94#         OS GB Specific Helpers for Generic Methods         #
    8495##############################################################
    8596
     
    94105        osgb_latlong = turn_xyz_into_llh(
    95106                                        osgb_xyz[0],osgb_xyz[1],osgb_xyz[2],'osgb')
    96 
    97107        return osgb_latlong
     108
    98109def turn_osgb36_into_wgs84(lat_dec,long_dec,height):
    99110        """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."""
     
    116127        """Turn OSGB36 easting and northing values into (decimal) lat/long values inOSGB36."""
    117128        return turn_eastingnorthing_into_latlong(easting,northing,'osgb')
     129
     130##############################################################
     131#         OS IE Specific Helpers for Generic Methods         #
     132##############################################################
     133
     134def 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
     146def 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
     159def 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
     163def 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')
    118166
    119167##############################################################
     
    323371
    324372        return (lat,long)
     373
     374##############################################################
     375#         Cassini Easting/Northing Transform Methods         #
     376##############################################################
     377
     378def turn_latlong_into_cassini_en(lat_dec,long_dec,scheme):
     379        """Latitude and Longitude, into Cassini-Soldner easting and northing co-ordinates, 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
     426def turn_cassini_en_into_latlong(easting,northing,scheme):
     427        """Cassini-Soldner 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 - ((1-e2) ** 0.5)) / (1 + ((1-e2) ** 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)
    325484
    326485##############################################################
  • trunk/npemap.org.uk/scripts/freethepostcode.org-importer/importer.py

    r179 r508  
    1010
    1111from pyPgSQL import PgSQL
    12 from geo_helper import turn_wgs84_into_osgb36, turn_osgb36_into_eastingnorthing
     12from geo_helper import turn_wgs84_into_osgb36, turn_osgb36_into_eastingnorthing, \
     13                                                turn_wgs84_into_osie36, turn_osie36_into_eastingnorthing
    1314import os
    1415
     
    1718dbname = "npemaps"
    1819dbhost = "localhost"
    19 dbuser = "npemap"
    20 dbpass = "npemap"
    21 
    22 
    23 raise Exception("Script needs to be updated as per ticket #39")
     20dbuser = "npemaps"
     21dbpass = ""
    2422
    2523
     
    8785                print "Invalid line '%s'" % line
    8886
    89         # Turn lat+long into easting+northing
    90         osll = turn_wgs84_into_osgb36(parts[0], parts[1], 0)
    91         en = turn_osgb36_into_eastingnorthing(osll[0], osll[1])
    92 
    9387        pc = {}
    94         pc["easting"] = en[0]
    95         pc["northing"] = en[1]
    9688        pc["outer"] = parts[2]
    9789        pc["inner"] = parts[3]
     
    9991        pc["raw_outer"] = parts[2]
    10092        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]
    101111
    102112        postcodes.append(pc)
     
    138148
    139149# 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)"
     150sql = "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)"
    141151sth = dbh.cursor()
    142152worked = 0
     153
     154# Add the postcodes
    143155for 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))
    145157        worked = worked + sth.rowcount
     158
    146159print "Added %d entries" % worked
    147160
Note: See TracChangeset for help on using the changeset viewer.