|
⇤ ← Revision 1 as of 2008-11-23 13:31:50
Size: 982
Comment:
|
Size: 3124
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 18: | Line 18: |
---- = 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 {{{#!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] }}} |
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 netcdf ice concentration data ftp server in polar stereographic projection
Read in the data using PyNIO
- To calculate the meridional ice concentration do the following:
- 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 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]