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.

Download ice concentration data ftp server
Read in the netcdf data using PyNIO
To calculate the meridional ice concentration do the following:
- Loop over latitude 75°N to 85°N
- Convert to polar stereographic coordinates [km].
- Convert to pixel indices (y,x). The resolution of the grid is 12.5 km. The offset to the pole is given in the definition polar stereographic projection

Polar stereographic projection

X,Y=mapll(lat,lon,sgn) converts from latitude and longitude to X,Y
lat,lon=mapxy(X,Y,sgn) converts from X,Y [km] to latitude and longitude
Set sgn=1 for north and sgn=-1 for southern hemisphere

   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]