Changeset d5dacc74a117f96a578399fc0a885e9b09731390


Ignore:
Timestamp:
09/21/07 18:30:38 (4 years ago)
Author:
flupke <flupke@…>
Children:
f1f525d6978dce4ec2af6d9a1b1202eaacc9b905
Parents:
fe16450ffbc9addda1be4cfcaad78b7c4ded9435
git-committer:
flupke <flupke@…> (09/21/07 18:30:38)
Message:

fixed projection code

git-svn-id:  http://kawa.selfip.org/svn/projects@444 532e76a9-45ae-c241-9c6d-8309ac6440cd

File:
1 edited

Legend:

Unmodified
Added
Removed

  • lambert/lambert.py

    r6832e21 rd5dacc7  
     1#-*- encoding: utf-8 -*- 
    12from math import * 
    23 
    … …  
    78PHI = 0       # Latitude 
    89LAMBDA = 1    # Longitude 
    9  
     10EPS10 = 1.e-10 
     11HALFPI = pi / 2.0 
     12QUARTERPI = pi / 4.0 
    1013 
    1114def gradesToRadians(grades): 
    … …  
    2124  """ Clamp negative numbers to 0.0 """ 
    2225  if number < 0.0: 
    23     print "Clamped", number 
    2426    return 0.0 
    2527  else: 
    … …  
    5557  """ Lambert conic conformal projection (2 standard parallels) parameters """ 
    5658   
    57   def __init__(self, ellipsoid, firstParallel, secondParallel, naturalOrigin, falseOrigin, shift): 
     59  def __init__(self, ellipsoid, firstParallel, secondParallel, origin, shift): 
    5860    self.ellipsoid = ellipsoid 
    5961    self.firstParallel = firstParallel 
    6062    self.secondParallel = secondParallel 
    61     self.naturalOrigin = naturalOrigin 
    62     self.falseOrigin = falseOrigin 
     63    self.origin = origin 
    6364    self.shift = shift 
     65   
     66  def __str__(self): 
     67    return "first parallel: %f\nsecond_parallel: %f" % \ 
     68        (degrees(self.firstParallel), degrees(self.secondParallel)) 
     69 
     70 
     71def msfn(sinphi, cosphi, es): 
     72  return (cosphi / sqrt(1.0 - es * sinphi * sinphi)) 
     73 
     74 
     75def tsfn(phi, sinphi, e): 
     76  sinphi *= e; 
     77  return (tan(0.5 * ((pi / 2.0) - phi)) / 
     78      pow((1. - sinphi) / (1. + sinphi), .5 * e)) 
    6479 
    6580 
    … …  
    7186  def __init__(self, parameters): 
    7287    self.setParameters(parameters) 
    73     self.project = self._projectOpti 
    74     self.unProject = self._unProjectOpti 
    7588 
    7689  def setParameters(self, parameters): 
    7790    self.parameters = parameters 
     91    projectionLatitude1 = parameters.firstParallel 
     92    projectionLatitude2 = parameters.secondParallel 
     93    projectionLatitude = parameters.origin[1] 
     94    e = self.parameters.ellipsoid.excentricity 
     95    es = self.parameters.ellipsoid.squaredExcentricity 
     96    self.n = sinphi = sin(projectionLatitude1) 
     97    cosphi = cos(projectionLatitude1) 
     98    secant = abs(projectionLatitude1 - projectionLatitude2) >= 1e-10 
     99    spherical = (es == 0.0) 
     100    if not spherical: 
     101      m1 = msfn(sinphi, cosphi, es) 
     102      ml1 = tsfn(projectionLatitude1, sinphi, e) 
     103      if secant: 
     104        sinphi = sin(projectionLatitude2) 
     105        self.n = log(m1 / 
     106            msfn(sinphi, cos(projectionLatitude2), es)) 
     107        self.n /= log(ml1 / tsfn(projectionLatitude2, sinphi, e)) 
     108      self.c = self.rho0 = m1 * pow(ml1, -self.n) / self.n 
     109      if abs(abs(projectionLatitude) - HALFPI) < 1e-10: 
     110        self.rho0 *= 0.0 
     111      else: 
     112        self.rho0 *= pow(tsfn(projectionLatitude, sin(projectionLatitude), e), self.n) 
     113    else: 
     114      if secant: 
     115        self.n = log(cosphi / cos(projectionLatitude2)) / \ 
     116            log(tan(QUARTERPI + .5 * projectionLatitude2) / \ 
     117            tan(QUARTERPI + .5 * projectionLatitude1)) 
     118      self.c = cosphi * pow(tan(QUARTERPI + .5 * projectionLatitude1), self.n) / self.n 
     119      if abs(abs(projectionLatitude) - HALFPI) < 1e-10: 
     120        self.rho0 = 0.0 
     121      else: 
     122        self.rho0 = self.c * pow(tan(QUARTERPI + .5 * projectionLatitude), -n) 
    78123 
    79     # Precalculate constant parameters 
    80     e = self.parameters.ellipsoid.excentricity 
    81     ee = self.parameters.ellipsoid.squaredExcentricity 
    82      
    83     f = lambda x: cos(x) / sqrt(1.0 - (ee * (sin(x)**2))) 
    84     m = {1: f(self.parameters.firstParallel), 2: f(self.parameters.secondParallel)} 
    85      
    86     f = lambda x: tan((pi / 4.0) - (x / 2.0)) / (((1.0 - (e * sin(x))) / (1.0 + (e * sin(x)))) ** (e / 2.0)) 
    87     t = {1: f(self.parameters.firstParallel), 
    88          2: f(self.parameters.secondParallel), 
    89          "F": clampPositive(f(self.parameters.falseOrigin[PHI]))} 
    90           
    91     self.n = (log(m[1]) - log(m[2])) / (log(t[1]) - log(t[2])) 
    92     self.F = m[1] / (self.n * (t[1]**self.n)) 
    93  
    94     f = lambda x: self.parameters.ellipsoid.semiMajor * self.F * (x**self.n) 
    95     self.r = {"F": f(t["F"])} 
    96  
    97   def _projectOpti(self, coordinates): 
    98     """ Projects geographical coordinates to (x, y) using Lambert projection. Optimized version. """ 
    99     e = self.parameters.ellipsoid.excentricity 
    100     t = clampPositive(tan((pi / 4.0) - (coordinates[PHI] / 2.0)) / (((1.0 - (e * sin(coordinates[PHI]))) / (1.0 + (e * sin(coordinates[PHI])))) ** (e / 2.0))) 
    101     r = self.parameters.ellipsoid.semiMajor * self.F * (t ** self.n) 
    102     theta = (self.n * (coordinates[LAMBDA] - self.parameters.naturalOrigin[LAMBDA])) - self._thetaAdjustment 
    103     return (self.parameters.shift[X] + (r * sin(theta)), self.parameters.shift[Y] + self.r["F"] - (r * cos(theta))) 
    104  
    105   def _projectNoOpti(self, coordinates): 
    106     """ Projects geographical coordinates to (x, y) using Lambert projection """ 
    107     e = self.parameters.ellipsoid.excentricity 
    108     ee = self.parameters.ellipsoid.squaredExcentricity 
    109     f = lambda x: cos(x) / sqrt(1.0 - (ee * (sin(x)**2))) 
    110     m = {1: f(self.parameters.firstParallel), 2: f(self.parameters.secondParallel)} 
    111     f = lambda x: tan((pi / 4.0) - (x / 2.0)) / (((1.0 - (e * sin(x))) / (1.0 + (e * sin(x)))) ** (e / 2.0)) 
    112     t = {1: f(self.parameters.firstParallel), 
    113          2: f(self.parameters.secondParallel), 
    114          "F": f(self.parameters.falseOrigin[PHI]), 
    115          "phi": f(coordinates[PHI])} 
    116     n = (log(m[1]) - log(m[2])) / (log(t[1]) - log(t[2])) 
    117     F = m[1] / (n * (t[1]**n)) 
    118     t["F"] = clampPositive(t["F"]) 
    119     t["phi"] = clampPositive(t["phi"]) 
    120     f = lambda x: self.parameters.ellipsoid.semiMajor * F * (x**n) 
    121     r = {"F": f(t["F"]), "phi": f(t["phi"])} 
    122     theta = (n * (coordinates[LAMBDA] - self.parameters.naturalOrigin[LAMBDA])) + self._thetaAdjustment 
    123     x = self.parameters.shift[X] + (r["phi"] * sin(theta)) 
    124     y = self.parameters.shift[Y] + r["F"] - (r["phi"] * cos(theta)) 
     124  def project(self, coordinates): 
     125    """ Projects geographical coordinates to (x, y) using Lambert projection. """ 
     126    coordinates = (coordinates[0] - self.parameters.origin[0], coordinates[1]) 
     127    if fabs(fabs(coordinates[1]) - (pi / 2.0)) < EPS10: 
     128      if coordinates[1] * self.n <= 0.0: 
     129        raise Exception 
     130      rho = 0.0 
     131    else: 
     132      if self.parameters.ellipsoid.squaredExcentricity == 0.0: 
     133        rho = self.c * pow(tan((pi / 4.0) + 0.5 * coordinates[1]), -n) 
     134      else: 
     135        rho = self.c * pow(tsfn(coordinates[1], sin(coordinates[1]), self.parameters.ellipsoid.excentricity), self.n) 
     136    a = coordinates[0] * self.n 
     137    x = self.parameters.ellipsoid.semiMajor * (rho * sin(a)) + self.parameters.shift[0] 
     138    y = self.parameters.ellipsoid.semiMajor * (self.rho0 - rho * cos(a)) + self.parameters.shift[1] 
    125139    return (x, y) 
    126140 
    127   def _unProjectOpti(self, coordinates, epsilon=0.0000000001): 
    128     e = self.parameters.ellipsoid.excentricity 
    129     r = sqrt(((coordinates[X] - self.parameters.shift[X]) ** 2) + ((self.parameters.shift[Y] - coordinates[Y] + self.r["F"]) ** 2)) 
    130     t = (r / (self.parameters.ellipsoid.semiMajor * self.F)) ** (1.0 / self.n) 
    131     phi = (pi / 2.0) - (2.0 * atan(t)) 
    132     while True: 
    133       prevPhi = phi 
    134       phi = (pi / 2.0) - (2.0 * atan(t * ((1.0 - e * sin(prevPhi)) / (1.0 + e * sin(prevPhi))) ** (e / 2.0))) 
    135       if phi - prevPhi <= epsilon: 
    136         break 
    137     theta = atan((coordinates[X] - self.parameters.shift[X]) / (self.parameters.shift[Y] - coordinates[Y] + self.r["F"])) + self._thetaAdjustment 
    138     lambd = theta / self.n + self.parameters.naturalOrigin[LAMBDA] 
    139     return (phi, lambd)     
     141 
     142  def inverse_project(self, coordinates): 
     143    """ Projects geographical coordinates to (x, y) using Lambert projection """ 
    140144 
    141145 
    … …  
    147151 
    148152STD_PROJECTIONS_PARAMETERS = { 
    149   "Lambert I": LambertProjectionParameters2SP( 
     153  "Lambert II Carto": LambertProjectionParameters2SP( 
    150154      ellipsoid=STD_ELLIPSOIDS["Clarke 1880 (IGN)"], 
    151       firstParallel=radians(sexagesimalToDecimal(degrees=48.0, minutes=35.0, seconds=54.682)),   # 48°35´54.682" 
    152       secondParallel=radians(sexagesimalToDecimal(degrees=50.0, minutes=23.0, seconds=45.282)),  # 50°23´45.282" 
    153       naturalOrigin=(gradesToRadians(55.0), 0.0), 
    154       falseOrigin=(gradesToRadians(55.0), 0.0), 
    155       shift=(600000.0, 200000.0)), 
    156   "Lambert II": LambertProjectionParameters2SP( 
    157       ellipsoid=STD_ELLIPSOIDS["Clarke 1880 (IGN)"], 
    158       firstParallel=radians(sexagesimalToDecimal(degrees=45.0, minutes=53.0, seconds=56.108)),   # 45°53´56.108" 
    159       secondParallel=radians(sexagesimalToDecimal(degrees=47.0, minutes=41.0, seconds=45.652)),  # 47°41´45.652" 
    160       naturalOrigin=(gradesToRadians(52.0), 0.0), 
    161       falseOrigin=(gradesToRadians(52.0), 0.0), 
    162       shift=(600000.0, 200000.0)), 
    163   "Lambert III": LambertProjectionParameters2SP( 
    164       ellipsoid=STD_ELLIPSOIDS["Clarke 1880 (IGN)"], 
    165       firstParallel=radians(sexagesimalToDecimal(degrees=43.0, minutes=11.0, seconds=57.449)),   # 43°11´57.449" 
    166       secondParallel=radians(sexagesimalToDecimal(degrees=44.0, minutes=59.0, seconds=45.938)),  # 44°59´45.938" 
    167       naturalOrigin=(gradesToRadians(49.0), 0.0), 
    168       falseOrigin=(gradesToRadians(49.0), 0.0), 
    169       shift=(600000.0, 200000.0)), 
    170   "Lambert IV": LambertProjectionParameters2SP( 
    171       ellipsoid=STD_ELLIPSOIDS["Clarke 1880 (IGN)"], 
    172       firstParallel=radians(sexagesimalToDecimal(degrees=41.0, minutes=33.0, seconds=37.396)),   # 41°33´37.396" 
    173       secondParallel=radians(sexagesimalToDecimal(degrees=42.0, minutes=46.0, seconds=3.588)),   # 42°46´03.588" 
    174       naturalOrigin=(gradesToRadians(46.85), 0.0), 
    175       falseOrigin=(gradesToRadians(46.85), 0.0), 
    176       shift=(234.358, 185861.369)), 
    177   "Lambert II extended": LambertProjectionParameters2SP( 
    178       ellipsoid=STD_ELLIPSOIDS["Clarke 1880 (IGN)"], 
    179       firstParallel=radians(sexagesimalToDecimal(degrees=45.0, minutes=53.0, seconds=56.108)),   # 44°53´56.108" 
    180       secondParallel=radians(sexagesimalToDecimal(degrees=47.0, minutes=41.0, seconds=45.652)),  # 47°41´45.652" 
    181       naturalOrigin=(gradesToRadians(52.0), 0.0), 
    182       falseOrigin=(gradesToRadians(52.0), 0.0), 
     155      firstParallel=radians(sexagesimalToDecimal(degrees=45.0, minutes=53.0, seconds=56.108)),   # 45°53'56.108" 
     156      secondParallel=radians(sexagesimalToDecimal(degrees=47.0, minutes=41.0, seconds=45.652)),  # 47°41'45.652" 
     157      origin=(radians(sexagesimalToDecimal(2.0, 20.0, 14.025)), radians(sexagesimalToDecimal(46.0, 48.0))), 
    183158      shift=(600000.0, 2200000.0)), 
    184   "Lambert 93": LambertProjectionParameters2SP( 
    185       ellipsoid=STD_ELLIPSOIDS["GRS 1980"], 
    186       firstParallel=radians(sexagesimalToDecimal(degrees=44.0)), 
    187       secondParallel=radians(sexagesimalToDecimal(degrees=49.0)), 
    188       naturalOrigin=(radians(sexagesimalToDecimal(degrees=46.0, minutes=30.0)), radians(sexagesimalToDecimal(degrees=3.0))), 
    189       falseOrigin=(radians(sexagesimalToDecimal(degrees=46.0, minutes=30.0)), radians(sexagesimalToDecimal(degrees=3.0))), 
    190       shift=(700000.0, 6600000.0)), 
    191159  } 
    192160 
    193161STD_PROJECTIONS = { 
    194   "Lambert I": LambertConicConformalProjection2SP(STD_PROJECTIONS_PARAMETERS["Lambert I"]), 
    195   "Lambert II": LambertConicConformalProjection2SP(STD_PROJECTIONS_PARAMETERS["Lambert II"]), 
    196   "Lambert III": LambertConicConformalProjection2SP(STD_PROJECTIONS_PARAMETERS["Lambert III"]), 
    197   "Lambert IV": LambertConicConformalProjection2SP(STD_PROJECTIONS_PARAMETERS["Lambert IV"]), 
    198   "Lambert II extended": LambertConicConformalProjection2SP(STD_PROJECTIONS_PARAMETERS["Lambert II extended"]), 
    199   "Lambert 93": LambertConicConformalProjection2SP(STD_PROJECTIONS_PARAMETERS["Lambert 93"]), 
     162  "Lambert II Carto": LambertConicConformalProjection2SP(STD_PROJECTIONS_PARAMETERS["Lambert II Carto"]), 
    200163  } 
    201164       
    … …  
    205168 
    206169if __name__ == "__main__": 
    207   parameters = LambertProjectionParameters2SP( 
    208     ellipsoid=STD_ELLIPSOIDS["International 1924"], 
    209     firstParallel=0.86975574, 
    210     secondParallel=0.89302680, 
    211     naturalOrigin=(1.57079633, 0.07604294), 
    212     falseOrigin=(1.57079633, 0.07604294), 
    213     shift=(150000.01, 5400088.44)) 
    214   projection = LambertConicConformalProjection2SPBelgium(parameters) 
    215   coords = (0.88452540, 0.10135773) 
     170  projection = STD_PROJECTIONS["Lambert II Carto"] 
     171  print projection.parameters 
     172  coords = (radians(5.82), radians(43.95)) 
    216173  result = projection.project(coords) 
    217   result2 = projection.unProject(result) 
    218   #result3 = projection.unProject((251763.20, 153034.13)) 
    219   print "coords: %s,\nprojected: %s,\nand back: %s" % (coords, result, result2) 
    220   assert str(result) == "(251763.20703772391, 153034.10206263792)" 
     174  #result = ((result[0] - 840005) / 16.666 , (result[1] - 1889995) / 16.666) 
     175  print "coords: %s,\nprojected: %s" % (coords, result) 
     176  #print "...", [degrees(x) for x in projection.unProject((840005, 1889995))] 
Note: See TracChangeset for help on using the changeset viewer.