|
Size: 3124
Comment:
|
Size: 3053
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 12: | Line 12: |
| * Download netcdf ice concentration data [[ftp://ftp.ifremer.fr/ifremer/cersat/products/gridded/psi-concentration/data|ftp server]] in [[http://nsidc.org/data/grids/ps_grid.html|polar stereographic projection]] * Read in the data using [[http://www.pyngl.ucar.edu/Nio.shtml|PyNIO]] |
* 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]] |
| Line 15: | Line 15: |
| * Loop over latitude 75N to 85N * 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]] |
* 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]] |
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
lat,lon=mapll(X,Y,sgn) converts from X,Y [km] to latitude and longitude
X,Y=mapxy(lat,lon,sgn) converts from latitude and longitude to X,Y
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]