LehreWiki

Data formats

Netcdf

What is the northernmost ice free position on a given longitude? Calculate the latitude for the prime meridian on the 15th September 2008.


Polar stereographic projection

   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]

Solution

   1 # imports
   2 from numpy import *
   3 from scipy import array,arange,nan
   4 from PyNGL import Ngl
   5 from PyNGL import Nio
   6 
   7 # Polar stereographic (NSDIDC grid)
   8 def mapxy(x, y, sgn):
   9     cdr  = 57.29577951
  10     slat = 70.0
  11     re   = 6378.273
  12     ec2   = .006693883
  13     ec    =  sqrt(ec2)
  14     if sgn == 1:
  15         delta = 45.0
  16     else:
  17         delta = 0.0
  18     sl = slat * pi/180.0
  19     rho = sqrt(x**2 + y**2)
  20     cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2))
  21     t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0)
  22     if  (absolute(slat-90.0) < 1.e-5):
  23         t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re
  24     else:
  25         t = rho * t / (re * cm)
  26     chi = (pi / 2.0) - 2.0 * arctan(t)
  27     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 * 
  28 chi)+ (7.0 * ec2**3.0 / 120.0) * sin(6.0 * chi)
  29     alat = sgn * alat
  30     along = arctan2(sgn * x, -sgn * y)
  31     along = sgn * along
  32 
  33     along = along * 180. / pi
  34     alat  = alat * 180. / pi
  35     along = along - delta
  36     return [alat,along]
  37 
  38 def  mapll(lat, lon,sgn):
  39     cdr  = 57.29577951
  40     slat = 70.
  41     re   = 6378.273
  42     ec2   = .006693883
  43     ec    =  sqrt(ec2)
  44     if sgn == 1:
  45         delta = 45.0
  46     else:
  47         delta = 0.0
  48     latitude  = absolute(lat) * pi/180.
  49     longitude = (lon + delta) * pi/180.
  50     t =   tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2)
  51     if ((90-slat) < 1.e-5):
  52         rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec))
  53     else:
  54         sl  = slat*pi/180.
  55         tc  = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.)
  56         mc  = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2))
  57         rho = re*mc*t/tc
  58 
  59     y=-rho*sgn*cos(sgn*longitude)
  60     x= rho*sgn*sin(sgn*longitude)
  61     return [x,y]
  62 
  63 def  mapij(x,y,sgn):
  64     #i=(x+3850)/12.5
  65     #j=(abs(y)+5850)/12.5
  66     #ir=round((x+3850)/12.5)
  67     #jr=round((abs(y)+5850)/12.5)
  68     i=(x+5850)/12.5
  69     j=(abs(y)+3850)/12.5
  70     ir=round((x+5850)/12.5)
  71     jr=round((abs(y)+3850)/12.5)
  72     return [i,j,ir,jr]
  73 
  74 # Main
  75 #nc = Nio.open_file('20080915.nc')
  76 nc = Nio.open_file('20081122.nc')
  77 C=array(nc.variables['concentration'])[0,:,:].astype(float)
  78 C1=array(nc.variables['concentration'])[0,:,:].astype(float)
  79 C1.fill(128.0)
  80 latList = linspace(75.0, 85.0, 100)
  81 resultList=[]
  82 print nc
  83 for i in latList:
  84     stereoCoords=mapll(i,0.0,1)
  85     print stereoCoords
  86     gridCoords=mapij(stereoCoords[0],stereoCoords[1],1)
  87     print gridCoords
  88     if C[gridCoords[2],gridCoords[3]] != 0.0 and C[gridCoords[2],gridCoords[3]] != 128.0   :
  89          resultList.append([i, C[gridCoords[2],gridCoords[3]],gridCoords[2],gridCoords[3]])
  90          C1[gridCoords[2],gridCoords[3]]=20.0
  91 resultList.sort()
  92 print resultList[0]

LehreWiki: Python/Exercise4 (last edited 2008-12-10 12:54:38 by SimonSchoof)