1 from numpy import *
   2 
   3 # Polar stereographic (NSDIDC grid)
   4 def mapxy(x, y, sgn):
   5     cdr  = 57.29577951
   6     slat = 70.0
   7     re   = 6378.273
   8     ec2   = .006693883
   9     ec    =  sqrt(ec2)
  10     if sgn == 1:
  11         delta = 45.0
  12     else:
  13         delta = 0.0
  14     sl = slat * pi/180.0
  15     rho = sqrt(x**2 + y**2)
  16     cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2))
  17     t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0)
  18     if  (absolute(slat-90.0) < 1.e-5):
  19         t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re
  20     else:
  21         t = rho * t / (re * cm)
  22     chi = (pi / 2.0) - 2.0 * arctan(t)
  23     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)
  24     alat = sgn * alat
  25     along = arctan2(sgn * x, -sgn * y)
  26     along = sgn * along
  27 
  28     along = along * 180. / pi
  29     alat  = alat * 180. / pi
  30     along = along - delta
  31     return [alat,along]
  32 
  33 def  mapll(lat, lon,sgn):
  34     cdr  = 57.29577951
  35     slat = 70.
  36     re   = 6378.273
  37     ec2   = .006693883
  38     ec    =  sqrt(ec2)
  39     if sgn == 1:
  40         delta = 45.0
  41     else:
  42         delta = 0.0
  43     latitude  = absolute(lat) * pi/180.
  44     longitude = (lon + delta) * pi/180.
  45     t =   tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2)
  46     if ((90-slat) < 1.e-5):
  47         rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec))
  48     else:
  49         sl  = slat*pi/180.
  50         tc  = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.)
  51         mc  = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2))
  52         rho = re*mc*t/tc
  53 
  54     y=-rho*sgn*cos(sgn*longitude)
  55     x= rho*sgn*sin(sgn*longitude)
  56     return [x,y]