#acl AdminGroup:read,write,delete,revert  EditorGroup:read,write StudentGroup:read,write All:read
#format wiki
#language de
#pragma section-numbers off

= 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://ftp.ifremer.fr/ifremer/cersat/products/gridded/psi-concentration/data|ftp server]] 
 * Read in the netcdf data using [[http://www.pyngl.ucar.edu/Nio.shtml|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 [[http://nsidc.org/data/grids/ps_grid.html|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

{{{#!python
from numpy import *

# Polar stereographic (NSDIDC grid)
def mapxy(x, y, sgn):
    cdr  = 57.29577951
    slat = 70.0
    re   = 6378.273
    ec2   = .006693883
    ec    =  sqrt(ec2)
    if sgn == 1:
        delta = 45.0
    else:
        delta = 0.0
    sl = slat * pi/180.0
    rho = sqrt(x**2 + y**2)
    cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2))
    t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0)
    if  (absolute(slat-90.0) < 1.e-5):
        t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re
    else:
        t = rho * t / (re * cm)
    chi = (pi / 2.0) - 2.0 * arctan(t)
    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)
    alat = sgn * alat
    along = arctan2(sgn * x, -sgn * y)
    along = sgn * along

    along = along * 180. / pi
    alat  = alat * 180. / pi
    along = along - delta
    return [alat,along]

def  mapll(lat, lon,sgn):
    cdr  = 57.29577951
    slat = 70.
    re   = 6378.273
    ec2   = .006693883
    ec    =  sqrt(ec2)
    if sgn == 1:
        delta = 45.0
    else:
        delta = 0.0
    latitude  = absolute(lat) * pi/180.
    longitude = (lon + delta) * pi/180.
    t =   tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2)
    if ((90-slat) < 1.e-5):
        rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec))
    else:
        sl  = slat*pi/180.
        tc  = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.)
        mc  = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2))
        rho = re*mc*t/tc

    y=-rho*sgn*cos(sgn*longitude)
    x= rho*sgn*sin(sgn*longitude)
    return [x,y]

}}}

=== Solution ===
{{{#!python
# imports
from numpy import *
from scipy import array,arange,nan
from PyNGL import Ngl
from PyNGL import Nio

# Polar stereographic (NSDIDC grid)
def mapxy(x, y, sgn):
    cdr  = 57.29577951
    slat = 70.0
    re   = 6378.273
    ec2   = .006693883
    ec    =  sqrt(ec2)
    if sgn == 1:
        delta = 45.0
    else:
        delta = 0.0
    sl = slat * pi/180.0
    rho = sqrt(x**2 + y**2)
    cm = cos(sl) / sqrt(1.0 - ec2 * (sin(sl)**2))
    t = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ec * sin(sl)) / (1.0 + ec * sin(sl)))**(ec / 2.0)
    if  (absolute(slat-90.0) < 1.e-5):
        t = rho * sqrt((1. + ec)**(1. + ec) * (1. - ec)**(1. - ec)) / 2. / re
    else:
        t = rho * t / (re * cm)
    chi = (pi / 2.0) - 2.0 * arctan(t)
    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)
    alat = sgn * alat
    along = arctan2(sgn * x, -sgn * y)
    along = sgn * along

    along = along * 180. / pi
    alat  = alat * 180. / pi
    along = along - delta
    return [alat,along]

def  mapll(lat, lon,sgn):
    cdr  = 57.29577951
    slat = 70.
    re   = 6378.273
    ec2   = .006693883
    ec    =  sqrt(ec2)
    if sgn == 1:
        delta = 45.0
    else:
        delta = 0.0
    latitude  = absolute(lat) * pi/180.
    longitude = (lon + delta) * pi/180.
    t =   tan(pi/4-latitude/2)/((1-ec*sin(latitude))/(1+ec*sin(latitude)))**(ec/2)
    if ((90-slat) < 1.e-5):
        rho = 2.*re*t/sqrt((1.+ec)**(1.+ec)*(1.-ec)**(1.-ec))
    else:
        sl  = slat*pi/180.
        tc  = tan(pi/4.-sl/2.)/((1.-ec*sin(sl))/(1.+ec*sin(sl)))**(ec/2.)
        mc  = cos(sl)/sqrt(1.0-ec2*(sin(sl)**2))
        rho = re*mc*t/tc

    y=-rho*sgn*cos(sgn*longitude)
    x= rho*sgn*sin(sgn*longitude)
    return [x,y]

def  mapij(x,y,sgn):
    #i=(x+3850)/12.5
    #j=(abs(y)+5850)/12.5
    #ir=round((x+3850)/12.5)
    #jr=round((abs(y)+5850)/12.5)
    i=(x+5850)/12.5
    j=(abs(y)+3850)/12.5
    ir=round((x+5850)/12.5)
    jr=round((abs(y)+3850)/12.5)
    return [i,j,ir,jr]

# Main
#nc = Nio.open_file('20080915.nc')
nc = Nio.open_file('20081122.nc')
C=array(nc.variables['concentration'])[0,:,:].astype(float)
C1=array(nc.variables['concentration'])[0,:,:].astype(float)
C1.fill(128.0)
latList = linspace(75.0, 85.0, 100)
resultList=[]
print nc
for i in latList:
    stereoCoords=mapll(i,0.0,1)
    print stereoCoords
    gridCoords=mapij(stereoCoords[0],stereoCoords[1],1)
    print gridCoords
    if C[gridCoords[2],gridCoords[3]] != 0.0 and C[gridCoords[2],gridCoords[3]] != 128.0   :
         resultList.append([i, C[gridCoords[2],gridCoords[3]],gridCoords[2],gridCoords[3]])
	 C1[gridCoords[2],gridCoords[3]]=20.0
resultList.sort()
print resultList[0]

}}}