Skip to content

Commit 9c9664e

Browse files
mayakbeckermayakbecker
committed
09_DataCorrection_Becker_Adusumilli added
1 parent c2c29af commit 9c9664e

3 files changed

Lines changed: 1205 additions & 0 deletions

File tree

‎09_DataCorrection_Becker_Adusumilli/README.md‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
Correcting ICESat-2 data and related applications
2+
-------------------------------------------------
3+
4+
Instructors: Maya Becker (m3becker@ucsd.edu) and Susheel Adusumilli (suadusum@ucsd.edu)
5+
6+
## Tutorials:
7+
8+
Part 1: Application of geophysical corrections over ice shelves
9+
Part 2: Derivation of ice-shelf height changes and conversion to mass changes
10+
11+
## Learning objectives:
12+
13+
- Establish an understanding of which geophysical corrections are suitable in which scenarios
14+
- Recognize limitations of geophysical corrections attached to ICESat-2 data products as well as external model outputs over Antarctic ice shelves
15+
- Mask grounded ice so as to focus on floating ice shelves, with the ability to extend this technique to other surfaces/terrains)
16+
- Load atmospheric products in NetCDF format to generate firn corrections
17+
- Convert height observations to ice-shelf thickness after applying the relevant corrections to ATL06 data over an ice-shelf area,
18+
- Estimate mass change from time series of ice-shelf height changes

‎09_DataCorrection_Becker_Adusumilli/data_correction_rendered.ipynb‎

Lines changed: 904 additions & 0 deletions
Large diffs are not rendered by default.

‎09_DataCorrection_Becker_Adusumilli/is2_utils.py‎

Lines changed: 283 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,283 @@
1+
import numpy as np
2+
from datetime import datetime as dt
3+
import time
4+
5+
CDR = 180. / np.pi
6+
E2 = 6.694379852*1e-3
7+
E = np.sqrt(E2)
8+
EARTH_RADIUS_KM = RE = 6378.1370
9+
EARTH_RADIUS_M = 6378137.0
10+
11+
def date2fracyr(dates):
12+
fracyr = np.empty_like(dates)
13+
i = 0
14+
for date in dates:
15+
def sinceEpoch(date): # returns seconds since epoch
16+
return time.mktime(date.timetuple())
17+
s = sinceEpoch
18+
19+
year = date.year
20+
startOfThisYear = dt(year=year, month=1, day=1)
21+
startOfNextYear = dt(year=year+1, month=1, day=1)
22+
23+
yearElapsed = s(date) - s(startOfThisYear)
24+
yearDuration = s(startOfNextYear) - s(startOfThisYear)
25+
fraction = yearElapsed/yearDuration
26+
fracyr[i] = date.year + fraction
27+
i+=1
28+
return fracyr
29+
30+
def ll2ps(lon, lat, slat=71, slon=0, hemi='s', units='km'):
31+
""" Spherical lon/lat -> Polar Steregraphic x/y.
32+
33+
This function converts from geodetic latitude and longitude to
34+
polar stereographic 'x/y' coordinates for the polar regions. The
35+
equations are from Snyder, J.P., 1982, Map Projections Used by
36+
the U.S. Geological Survey, Geological Survey Bulletin 1532, U.S.
37+
Government Printing Office. See JPL Technical Memorandum
38+
3349-85-101 for further details.
39+
40+
Parameters
41+
----------
42+
lon, lat : array-like (1d or 2d) or float
43+
Geodetic longitude and latitude (degrees, -/+180 or 0/360 and -/+90).
44+
slat : float
45+
Standard latitude (e.g., 70 S), see Notes.
46+
slon : float
47+
Standard longitude (e.g., 0), see Notes.
48+
hemi : string
49+
Hemisphere: 'n' or 's' (not case-sensitive).
50+
units : string
51+
Polar Stereographic x/y units: 'm' or 'km' (not case-sensitive).
52+
53+
Returns
54+
-------
55+
x, y : ndarray (1d or 2d) or float
56+
Polar stereographic x and y coordinates (in 'm' or 'km').
57+
Notes
58+
-----
59+
SLAT is is the "true" latitude in the plane projection
60+
(the map), so there is no deformation over this latitude;
61+
e.g., using the same SLON but changing SLAT from 70 to 71
62+
degrees, will move things in polar stereo. The goal is to
63+
locally preserve area and angles. Most users use 71S but
64+
the sea ice people use 70S.
65+
66+
SLON provides a "vertical" coordinate for plotting and for
67+
rectangle orientation. E.g., for Arctic sea ice, NSIDC use
68+
SLON=45 in order to make a grid that is optimized for where
69+
sea ice occurs. Ex: CATS2008a has SLON=-70 (AP roughly up),
70+
so that the grid can be long enough to include South Georgia.
71+
Other examples are:
72+
MOA Image Map (the GeoTIFF): SLAT=-71, SLON=0
73+
MOA mask grid (from Laurie): SLAT=-71, SLON=-70
74+
Scripps mask grid (from GeoTIFF): SLAT=-71, SLON=0
75+
History
76+
-------
77+
Written in Fortran by C.S. Morris - Apr 29, 1985
78+
Revised by C.S. Morris - Dec 11, 1985
79+
Revised by V.J. Troisi - Jan 1990
80+
SGN - provides hemisphere dependency (+/- 1)
81+
Revised by Xiaoming Li - Oct 1996
82+
Corrected equation for RHO
83+
Converted to Matlab by L. Padman - Oct 25, 2006
84+
Updated for slon by L. Padman - Nov 21, 2006
85+
Converted to Python by F.S. Paolo - Mar 23, 2010
86+
87+
Example
88+
-------
89+
>>> lon = [-150.3, 66.2, 5.3]
90+
>>> lat = [70.2, 75.5, 80.3]
91+
>>> x, y = ll2ps(lon, lat, slat=71, slon=-70, hemi='s', units='m')
92+
Original (Matlab) documentation
93+
-------------------------------
94+
ARGUMENTS:
95+
96+
Variable I/O Description
97+
98+
lat I Geodetic Latitude (degrees, +90 to -90)
99+
lon I Geodetic Longitude (degrees, 0 to 360)
100+
SLAT I Standard latitude (typ. 71, or 70)
101+
SLON I
102+
HEMI I Hemisphere (char*1: 'N' or 'S' (not
103+
case-sensitive)
104+
x O Polar Stereographic X Coordinate (km)
105+
y O Polar Stereographic Y Coordinate (km)
106+
"""
107+
if units != 'm':
108+
units = 'km'
109+
110+
# print("converting lon/lat -> x/y ...")
111+
112+
# if sequence, convert to ndarray double
113+
if type(lon).__name__ in ['list', 'tuple']:
114+
lon = np.array(lon, 'f8')
115+
lat = np.array(lat, 'f8')
116+
117+
# if ndarray, convert to double if it isn't
118+
if type(lon).__name__ == 'ndarray' and lon.dtype != 'float64':
119+
lon = lon.astype(np.float64)
120+
lat = lat.astype(np.float64)
121+
122+
# convert longitude
123+
if type(lon).__name__ == 'ndarray': # is numpy array
124+
lon[lon<0] += 360. # -/+180 -> 0/360
125+
elif lon.any() < 0: # is scalar
126+
lon += 360.
127+
128+
if (str.lower(hemi) == 's'):
129+
SGN = -1
130+
else:
131+
SGN = 1
132+
if (np.abs(slat) == 90):
133+
RHO = 2. * RE / ((1 + E)**(1 + E) * (1 - E)**(1 - E))**(E/2.)
134+
else:
135+
SL = np.abs(slat) / CDR
136+
TC = np.tan(np.pi/4. - SL/2.) / ((1 - E * np.sin(SL)) \
137+
/ (1 + E * np.sin(SL)))**(E/2.)
138+
MC = np.cos(SL) / np.sqrt(1 - E2 * (np.sin(SL)**2))
139+
RHO = RE * MC / TC
140+
141+
lat = np.abs(lat) / CDR
142+
T = np.tan(np.pi/4. - lat/2.) / ((1 - E * np.sin(lat)) \
143+
/ (1 + E * np.sin(lat)))**(E/2.)
144+
145+
lon2 = -(lon - slon) / CDR
146+
x = -RHO * T * np.sin(lon2) # global vars
147+
y = RHO * T * np.cos(lon2)
148+
149+
if units == 'm': # computations are done in km
150+
x *= 1e3
151+
y *= 1e3
152+
153+
return [x, y]
154+
155+
156+
def ps2ll(x, y, slat=71, slon=0, hemi='s', units='km'):
157+
"""Polar Stereographic x/y -> Spherical lon/lat.
158+
159+
This subroutine converts from Polar Stereographic 'x,y' coordinates
160+
to geodetic longitude and latitude for the polar regions. The
161+
equations are from Snyder, J.P., 1982, Map Projections Used by the
162+
U.S. Geological Survey, Geological Survey Bulletin 1532, U.S.
163+
Government Printing Office. See JPL Technical Memorandum
164+
3349-85-101 for further details.
165+
166+
Parameters
167+
----------
168+
x, y : array-like (1d or 2d) or float
169+
Polar stereographic x and y coordinates (in 'm' or 'km').
170+
slat : float
171+
Standard latitude (e.g., 70 S), see Notes.
172+
slon : float
173+
Standard longitude (e.g., 0), see Notes.
174+
hemi : string
175+
Hemisphere: 'n' or 's' (not case-sensitive).
176+
units : string
177+
Polar Stereographic x/y units: 'm' or 'km' (not case-sensitive).
178+
179+
Returns
180+
-------
181+
lon, lat : ndarray (1d or 2d) or float
182+
Geodetic longitude and latitude (degrees, 0/360 and -/+90).
183+
184+
Notes
185+
-----
186+
SLAT is the "true" latitude in the plane projection
187+
(the map), so there is no deformation over this latitude;
188+
e.g., using the same SLON but changing SLAT from 70 to 71
189+
degrees, will move things in polar stereo. The goal is to
190+
locally preserve area and angles. Most users use 71S but
191+
the sea ice people use 70S.
192+
193+
SLON provides a "vertical" coordinate for plotting and for
194+
rectangle orientation. E.g., for Arctic sea ice, NSIDC use
195+
SLON=45 in order to make a grid that is optimized for where
196+
sea ice occurs. CATS2008a has SLON=-70 (AP roughly up), so
197+
that the grid can be long enough to include South Georgia.
198+
MOA Image Map (the GeoTIFF): SLAT=-71, SLON=0
199+
MOA mask grid (from Laurie): SLAT=-71, SLON=-70
200+
Scripps mask grid (from GeoTIFF): SLAT=-71, SLON=0
201+
History
202+
-------
203+
Written in Fortran by C.S. Morris - Apr 29, 1985
204+
Revised by C.S. Morris - Dec 11, 1985
205+
Revised by V.J. Troisi - Jan 1990
206+
SGN - provides hemisphere dependency (+/- 1)
207+
Converted to Matlab by L. Padman - Oct 25, 2006
208+
Updated for slon by L. Padman - Nov 21, 2006
209+
Converted to Python by F.S. Paolo - Mar 23, 2010
210+
Edited Susheel for is2py Sep 09, 2018
211+
212+
Example
213+
-------
214+
>>> x = [-2141.06767831 1096.06628549 1021.77465469]
215+
>>> y = [ 365.97940112 -1142.96735458 268.05756254]
216+
>>> lon, lat = xy2ll(x, y, slat=71, slon=-70, hemi='s', units='km')
217+
Original (Matlab) documentation
218+
-------------------------------
219+
ARGUMENTS:
220+
221+
Variable I/O Description
222+
223+
X I Polar Stereographic X Coordinate (km)
224+
Y I Polar Stereographic Y Coordinate (km)
225+
SLAT I Standard latitude (typ. 71, or 70)
226+
SLON I Standard longitude
227+
HEMI I Hemisphere (char*1, 'S' or 'N',
228+
not case-sensitive)
229+
lat O Geodetic Latitude (degrees, +90 to -90)
230+
lon O Geodetic Longitude (degrees, 0 to 360)
231+
"""
232+
if units != 'm':
233+
units = 'km'
234+
235+
# print("converting 'x,y' -> 'lon,lat' ...")
236+
237+
# if sequence, convert to ndarray
238+
if type(x).__name__ in ['list', 'tuple']:
239+
x = np.array(x, 'f8')
240+
y = np.array(y, 'f8')
241+
242+
# if ndarray, convert to double if it isn't
243+
if type(x).__name__ == 'ndarray' and x.dtype != 'float64':
244+
x = x.astype(np.float64)
245+
y = y.astype(np.float64)
246+
247+
if units == 'm': # computations are done in km !!!
248+
x *= 1e-3
249+
y *= 1e-3
250+
251+
if(str.lower(hemi) == 's'):
252+
SGN = -1.
253+
else:
254+
SGN = 1.
255+
slat = np.abs(slat)
256+
SL = slat / CDR
257+
RHO = np.sqrt(x**2 + y**2) # if scalar, is numpy.float64
258+
if np.alltrue(RHO < 0.1): # Don't calculate if "all points" on the equator
259+
lat = 90.0 * SGN
260+
lon = 0.0
261+
return lon, lat
262+
else:
263+
CM = np.cos(SL) / np.sqrt(1 - E2 * (np.sin(SL)**2))
264+
T = np.tan((np.pi/4.) - (SL/2.)) / ((1 - E * np.sin(SL)) \
265+
/ (1 + E * np.sin(SL)))**(E/2.)
266+
if (np.abs(slat - 90.) < 1.e-5):
267+
T = ((RHO * np.sqrt((1 + E)**(1 + E) * (1 - E)**(1 - E))) / 2.) / RE
268+
else:
269+
T = RHO * T / (RE * CM)
270+
271+
a1 = 5 * E2**2 / 24.
272+
a2 = 1 * E2**3 / 12.
273+
a3 = 7 * E2**2 / 48.
274+
a4 = 29 * E2**3 / 240.
275+
a5 = 7 * E2**3 / 120.
276+
277+
CHI = (np.pi/2.) - 2. * np.arctan(T)
278+
lat = CHI + ((E2/2.) + a1 + a2) * np.sin(2. * CHI) \
279+
+ (a3 + a4) * np.sin(4. * CHI) + a5 * np.sin(6. * CHI)
280+
lat = SGN * lat * CDR
281+
lon = -(np.arctan2(-x, y) * CDR) + slon
282+
283+
return [lon, lat]

0 commit comments

Comments
 (0)