Skip to content

Commit dde00b1

Browse files
jpswinskijpswinski
committed
Merge branch 'main' of github.com:ICESat2-SlideRule/sliderule-python into main
2 parents 4d3ba9c + 7ba740f commit dde00b1

File tree

4 files changed

+434
-28
lines changed

4 files changed

+434
-28
lines changed

‎examples/cmr_debug_regions.ipynb

Lines changed: 344 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,344 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"### SlideRule CMR Debug\n",
8+
"A simple tool for checking CMR queries\n",
9+
"- Creates a leaflet map for creating regions of interest\n",
10+
"- Queries CMR for granules and granule polygons\n",
11+
"- Plots granule polygons on map\n",
12+
"- Retrieves and plots granule tracks on map"
13+
]
14+
},
15+
{
16+
"cell_type": "markdown",
17+
"metadata": {},
18+
"source": [
19+
"#### Load necessary packages"
20+
]
21+
},
22+
{
23+
"cell_type": "code",
24+
"execution_count": null,
25+
"metadata": {},
26+
"outputs": [],
27+
"source": [
28+
"import re\n",
29+
"import time\n",
30+
"import logging\n",
31+
"import posixpath\n",
32+
"import numpy as np\n",
33+
"import geopandas as gpd\n",
34+
"import concurrent.futures\n",
35+
"import matplotlib.cm as cm\n",
36+
"import matplotlib.colors as colors\n",
37+
"import ipywidgets as widgets\n",
38+
"from shapely.geometry import LineString\n",
39+
"from sliderule import sliderule, icesat2, ipysliderule\n",
40+
"import sliderule.io\n",
41+
"# autoreload\n",
42+
"%load_ext autoreload\n",
43+
"%autoreload 2\n",
44+
"# create logger\n",
45+
"logging.basicConfig(level=logging.INFO)"
46+
]
47+
},
48+
{
49+
"cell_type": "markdown",
50+
"metadata": {},
51+
"source": [
52+
"#### Set ICESat-2 Product"
53+
]
54+
},
55+
{
56+
"cell_type": "code",
57+
"execution_count": null,
58+
"metadata": {},
59+
"outputs": [],
60+
"source": [
61+
"# Configure ICESat-2 API\n",
62+
"icesat2.init(\"icesat2sliderule.org\", loglevel=logging.WARNING)\n",
63+
"icesat2.get_version()\n",
64+
"\n",
65+
"# dropdown menu for ICESat-2 product\n",
66+
"product_select = widgets.Dropdown(\n",
67+
" options=['ATL03','ATL06','ATL08'],\n",
68+
" description='Product:',\n",
69+
" disabled=False\n",
70+
")\n",
71+
"display(product_select)"
72+
]
73+
},
74+
{
75+
"cell_type": "markdown",
76+
"metadata": {},
77+
"source": [
78+
"#### Select regions of interest for querying CMR"
79+
]
80+
},
81+
{
82+
"cell_type": "code",
83+
"execution_count": null,
84+
"metadata": {},
85+
"outputs": [],
86+
"source": [
87+
"# create ipyleaflet map in projection\n",
88+
"m = ipysliderule.leaflet('Global', center=(-25.3,131))\n",
89+
"m.map"
90+
]
91+
},
92+
{
93+
"cell_type": "markdown",
94+
"metadata": {},
95+
"source": [
96+
"#### Build and transmit CMR requests"
97+
]
98+
},
99+
{
100+
"cell_type": "code",
101+
"execution_count": null,
102+
"metadata": {},
103+
"outputs": [],
104+
"source": [
105+
"# for each region of interest\n",
106+
"granule_list = []\n",
107+
"granule_polygons = []\n",
108+
"for poly in m.regions:\n",
109+
" # polygon from map\n",
110+
" resources,polygons = icesat2.cmr(polygon=poly,\n",
111+
" short_name=product_select.value,\n",
112+
" return_polygons=True)\n",
113+
" granule_list.extend(resources)\n",
114+
" granule_polygons.extend(polygons)\n",
115+
"# print list of granules\n",
116+
"num_granules = len(granule_list)\n",
117+
"logging.info('Number of Granules: {0:d}'.format(num_granules))\n",
118+
"logging.debug(granule_list)"
119+
]
120+
},
121+
{
122+
"cell_type": "markdown",
123+
"metadata": {},
124+
"source": [
125+
"#### Select Granules to Plot on Map"
126+
]
127+
},
128+
{
129+
"cell_type": "code",
130+
"execution_count": null,
131+
"metadata": {},
132+
"outputs": [],
133+
"source": [
134+
"granule_select = widgets.SelectMultiple(\n",
135+
" options=granule_list,\n",
136+
" description='Granules:',\n",
137+
" layout=widgets.Layout(width='35%', height='200px'),\n",
138+
" disabled=False\n",
139+
")\n",
140+
"display(granule_select)"
141+
]
142+
},
143+
{
144+
"cell_type": "markdown",
145+
"metadata": {},
146+
"source": [
147+
"#### Add granule polygons to map"
148+
]
149+
},
150+
{
151+
"cell_type": "code",
152+
"execution_count": null,
153+
"metadata": {},
154+
"outputs": [],
155+
"source": [
156+
"granule_indices = list(granule_select.index)\n",
157+
"cmap = iter(cm.viridis(np.linspace(0,1,len(granule_indices))))\n",
158+
"for g in granule_indices:\n",
159+
" locations = [(p['lat'],p['lon']) for p in granule_polygons[g]]\n",
160+
" color = colors.to_hex(next(cmap))\n",
161+
" polygon = ipysliderule.ipyleaflet.Polygon(\n",
162+
" locations=locations,\n",
163+
" color=color,\n",
164+
" fill_color=color,\n",
165+
" opacity=0.8,\n",
166+
" weight=1, \n",
167+
" )\n",
168+
" m.map.add_layer(polygon)"
169+
]
170+
},
171+
{
172+
"cell_type": "markdown",
173+
"metadata": {},
174+
"source": [
175+
"#### Get Granules from NSIDC S3"
176+
]
177+
},
178+
{
179+
"cell_type": "code",
180+
"execution_count": null,
181+
"metadata": {},
182+
"outputs": [],
183+
"source": [
184+
"def s3_retrieve(granule, **kwargs):\n",
185+
" # set default keyword arguments\n",
186+
" kwargs.setdefault('asset','nsidc-s3')\n",
187+
" kwargs.setdefault('index_key','time')\n",
188+
" kwargs.setdefault('polygon',None)\n",
189+
" # regular expression operator for extracting information from files\n",
190+
" rx = re.compile(r'(ATL\\d{2})(-\\d{2})?_(\\d{4})(\\d{2})(\\d{2})(\\d{2})'\n",
191+
" r'(\\d{2})(\\d{2})_(\\d{4})(\\d{2})(\\d{2})_(\\d{3})_(\\d{2})(.*?).h5$')\n",
192+
" # extract parameters from ICESat-2 granule\n",
193+
" PRD,HEM,YY,MM,DD,HH,MN,SS,TRK,CYCL,GRN,RL,VRS,AUX=rx.findall(granule).pop()\n",
194+
" # variables of interest\n",
195+
" if (PRD == 'ATL03'):\n",
196+
" segment_group = \"geolocation\"\n",
197+
" segment_key = 'segment_id'\n",
198+
" lon_key = 'reference_photon_lon'\n",
199+
" lat_key = 'reference_photon_lat'\n",
200+
" vnames = ['segment_id','delta_time','reference_photon_lat',\n",
201+
" 'reference_photon_lon']\n",
202+
" elif (PRD == 'ATL06'):\n",
203+
" segment_group = \"land_ice_segments\"\n",
204+
" segment_key = 'segment_id'\n",
205+
" lon_key = 'longitude'\n",
206+
" lat_key = 'latitude'\n",
207+
" vnames = ['segment_id','delta_time','latitude','longitude']\n",
208+
" elif (PRD == 'ATL08'):\n",
209+
" segment_group = \"land_segments\"\n",
210+
" segment_key = 'segment_id_beg'\n",
211+
" lon_key = 'longitude'\n",
212+
" lat_key = 'latitude'\n",
213+
" vnames = ['segment_id_beg','segment_id_end','delta_time',\n",
214+
" 'latitude','longitude']\n",
215+
" # for each valid beam within the HDF5 file\n",
216+
" frames = []\n",
217+
" gt = dict(gt1l=10,gt1r=20,gt2l=30,gt2r=40,gt3l=50,gt3r=60)\n",
218+
" atlas_sdp_epoch = np.datetime64('2018-01-01T00:00:00')\n",
219+
" kwds = dict(startrow=0,numrows=-1)\n",
220+
" for gtx in ['gt1l','gt1r','gt2l','gt2r','gt3l','gt3r']:\n",
221+
" geodatasets = [dict(dataset=f'{gtx}/{segment_group}/{v}',**kwds) for v in vnames]\n",
222+
" try:\n",
223+
" # get datasets from s3\n",
224+
" hidatasets = icesat2.h5p(geodatasets, granule, kwargs['asset'])\n",
225+
" # copy to new \"flattened\" dictionary\n",
226+
" data = {posixpath.basename(key):var for key,var in hidatasets.items()}\n",
227+
" # Generate Time Column\n",
228+
" delta_time = (data['delta_time']*1e9).astype('timedelta64[ns]')\n",
229+
" data['time'] = gpd.pd.to_datetime(atlas_sdp_epoch + delta_time)\n",
230+
" except:\n",
231+
" pass\n",
232+
" else:\n",
233+
" # copy filename parameters\n",
234+
" data['rgt'] = [int(TRK)]*len(data['delta_time'])\n",
235+
" data['cycle'] = [int(CYCL)]*len(data['delta_time'])\n",
236+
" data['gt'] = [gt[gtx]]*len(data['delta_time'])\n",
237+
" # pandas dataframe from compiled dictionary\n",
238+
" frames.append(gpd.pd.DataFrame.from_dict(data))\n",
239+
" # concatenate pandas dataframe\n",
240+
" try:\n",
241+
" df = gpd.pd.concat(frames)\n",
242+
" except:\n",
243+
" return sliderule.icesat2.__emptyframe()\n",
244+
" # convert to a GeoDataFrame\n",
245+
" geometry = gpd.points_from_xy(df[lon_key], df[lat_key])\n",
246+
" gdf = gpd.GeoDataFrame(df.drop(columns=[lon_key,lat_key]),\n",
247+
" geometry=geometry,crs='EPSG:4326')\n",
248+
" # create global track variable\n",
249+
" track = 6*(gdf['rgt']-1) + (gdf['gt']/10)\n",
250+
" gdf = gdf.assign(track=track)\n",
251+
" # calculate global reference point\n",
252+
" global_ref_pt = 6*1387*gdf[segment_key] + track\n",
253+
" gdf = gdf.assign(global_ref_pt=global_ref_pt)\n",
254+
" # sort values for reproducible output despite async processing\n",
255+
" gdf.set_index(kwargs['index_key'], inplace=True)\n",
256+
" gdf.sort_index(inplace=True)\n",
257+
" # remove duplicate points\n",
258+
" gdf = gdf[~gdf.index.duplicated()]\n",
259+
" # intersect with geometry in projected reference system\n",
260+
" if kwargs['polygon'] is not None:\n",
261+
" gdf = gpd.overlay(gdf.to_crs(kwargs['polygon'].crs),\n",
262+
" kwargs['polygon'], how='intersection')\n",
263+
" # convert back to original coordinate reference system\n",
264+
" return gdf.to_crs('EPSG:4326')"
265+
]
266+
},
267+
{
268+
"cell_type": "code",
269+
"execution_count": null,
270+
"metadata": {},
271+
"outputs": [],
272+
"source": [
273+
"# granule resources for selected segments\n",
274+
"perf_start = time.perf_counter()\n",
275+
"gdf = sliderule.icesat2.__emptyframe()\n",
276+
"num_servers = sliderule.update_available_servers()\n",
277+
"max_workers = num_servers * icesat2.SERVER_SCALE_FACTOR\n",
278+
"with concurrent.futures.ThreadPoolExecutor(max_workers=max_workers) as executor:\n",
279+
" futures = [executor.submit(s3_retrieve, granule_list[g]) for g in granule_indices]\n",
280+
" # Wait for Results\n",
281+
" for future in concurrent.futures.as_completed(futures):\n",
282+
" # append to dataframe\n",
283+
" gdf = gdf.append(future.result())\n",
284+
"\n",
285+
"# Display Statistics\n",
286+
"perf_stop = time.perf_counter()\n",
287+
"perf_duration = perf_stop - perf_start\n",
288+
"print(\"Completed in {:.3f} seconds of wall-clock time\".format(perf_duration))\n",
289+
"print(\"Reference Ground Tracks: {}\".format(gdf[\"rgt\"].unique()))\n",
290+
"print(\"Cycles: {}\".format(gdf[\"cycle\"].unique()))\n",
291+
"print(\"Received {} elevations\".format(gdf.shape[0]))"
292+
]
293+
},
294+
{
295+
"cell_type": "markdown",
296+
"metadata": {},
297+
"source": [
298+
"#### Add Granule Track Polylines to Map"
299+
]
300+
},
301+
{
302+
"cell_type": "code",
303+
"execution_count": null,
304+
"metadata": {},
305+
"outputs": [],
306+
"source": [
307+
"# fix int columns that were converted in objects\n",
308+
"fixed = gdf.drop(columns=['geometry'])\n",
309+
"for column in fixed.select_dtypes(include='object').columns:\n",
310+
" fixed[column] = fixed[column].astype(\"int32\")\n",
311+
"fixed = gpd.GeoDataFrame(fixed,geometry=gdf.geometry,crs='EPSG:4326')\n",
312+
"# convert from points to linestrings grouping by track\n",
313+
"grouped = fixed.groupby(['track'])['geometry'].apply(\n",
314+
" lambda x: LineString(x.tolist()) if x.size > 1 else x.tolist())\n",
315+
"geodata = ipysliderule.ipyleaflet.GeoData(geo_dataframe=gpd.GeoDataFrame(grouped),\n",
316+
" style={'color': 'black', 'opacity':1, 'weight':0.1})\n",
317+
"m.map.add_layer(geodata)"
318+
]
319+
}
320+
],
321+
"metadata": {
322+
"interpreter": {
323+
"hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
324+
},
325+
"kernelspec": {
326+
"display_name": "Python 3.8.10 64-bit",
327+
"name": "python3"
328+
},
329+
"language_info": {
330+
"codemirror_mode": {
331+
"name": "ipython",
332+
"version": 3
333+
},
334+
"file_extension": ".py",
335+
"mimetype": "text/x-python",
336+
"name": "python",
337+
"nbconvert_exporter": "python",
338+
"pygments_lexer": "ipython3",
339+
"version": "3.8.10"
340+
}
341+
},
342+
"nbformat": 4,
343+
"nbformat_minor": 4
344+
}

0 commit comments

Comments
 (0)