1 #module general/polar_projection.py
   2 # LK 14/08/2008,4/2/2009
   3 
   4 from numpy import *
   5 
   6 # Polar stereographic (NSDIDC grid)
   7 def mapxy(x, y, sgn):
   8     cdr  = 57.29577951
   9     slat = 70.0
  10     re   = 6378.273
  11     ec2   = .006693883
  12     ec    =  sqrt(ec2)
  13     if sgn == 1:
  14         delta = 45.0
  15     else:
  16         delta = 0.0
  17     sl = slat * pi/180.0
  18     rho = sqrt(x**2 + y**2)
  19     cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2))
  20     t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0)
  21     if  (absolute(slat-90.0) < 1.e-5):
  22         t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re
  23     else:
  24         t = rho * t / (re * cm)
  25     chi = (pi / 2.0) - 2.0 * arctan(t)
  26     alat = chi + ((ec2 / 2.0) + (5.0 * ec2**2.0 / 24.0) + (ec2**3.0 / 12.0)) * sin(2 * chi) + ((7.0 * ec2**2.0 / 48.0) + (29.0 * ec2**3 / 240.0)) *sin(4.0 * chi)+ (7.0 * ec2**3.0 / 120.0) * sin(6.0 * chi)
  27     alat = sgn * alat
  28     along = arctan2(sgn * x, -sgn * y)
  29     along = sgn * along
  30 
  31     along = along * 180. / pi
  32     alat  = alat * 180. / pi
  33     along = along - delta
  34     return [alat,along]
  35 
  36 def  mapll(lat, lon,sgn):
  37     cdr  = 57.29577951
  38     slat = 70.
  39     re   = 6378.273
  40     ec2   = .006693883
  41     ec    =  sqrt(ec2)
  42     if sgn == 1:
  43         delta = 45.0
  44     else:
  45         delta = 0.0
  46     latitude  = absolute(lat) * pi/180.
  47     longitude = (lon + delta) * pi/180.
  48     t =   tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2)
  49     if ((90-slat) < 1.e-5):
  50         rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec))
  51     else:
  52         sl  = slat*pi/180.
  53         tc  = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.)
  54         mc  = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2))
  55         rho = re*mc*t/tc
  56 
  57     y=-rho*sgn*cos(sgn*longitude)
  58     x= rho*sgn*sin(sgn*longitude)
  59     return [x,y]
  60 
  61 def  mapll2(lat, lon):
  62     if lat>0:
  63         sgn=1
  64     else:
  65         sgn=-1
  66     return mapll(lat,lon,sgn)
  67 
  68 def dist(lat1,lon1,lat2,lon2):
  69     x1,y1=mapll2(lat1,lon1)
  70     x2,y2=mapll2(lat2,lon2)
  71     #print x1,y1,x2,y2
  72     d=sqrt((x2-x1)**2+(y2-y1)**2)
  73     if sign(lat1)!=sign(lat2):
  74         d=10000.0
  75     return d