Skip to content

Commit f464a63

Browse files
jpswinskijpswinski
committed
merging latest from development
2 parents c8cec49 + 9a5dd47 commit f464a63

File tree

11 files changed

+361
-72
lines changed

11 files changed

+361
-72
lines changed

‎environment.yml

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,9 +6,9 @@ dependencies:
66
- geopandas
77
- h5py
88
- ipykernel
9-
- ipywidgets>=7.6,<8.0
10-
- ipyleaflet>=0.15
11-
- jupyterlab=3
9+
- ipywidgets
10+
- ipyleaflet
11+
- jupyterlab
1212
- jupyterlab_widgets
1313
- matplotlib
1414
- ipympl
@@ -25,5 +25,6 @@ dependencies:
2525
- shapely
2626
- tk
2727
- xyzservices
28+
- pyarrow
2829
- pip:
2930
- -e ./

‎examples/arcticdem.ipynb

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": null,
6+
"id": "acb12a75-1636-471a-9649-48a408801d4f",
7+
"metadata": {},
8+
"outputs": [],
9+
"source": [
10+
"from sliderule import icesat2\n",
11+
"import matplotlib.pyplot as plt\n",
12+
"import matplotlib\n",
13+
"import geopandas"
14+
]
15+
},
16+
{
17+
"cell_type": "code",
18+
"execution_count": null,
19+
"id": "b8167cbe-3fe3-4dc9-a5ad-0cbba51c8a07",
20+
"metadata": {},
21+
"outputs": [],
22+
"source": [
23+
"icesat2.init(\"slideruleearth.io\", verbose=True)"
24+
]
25+
},
26+
{
27+
"cell_type": "code",
28+
"execution_count": null,
29+
"id": "98ef750a-e88b-4125-b951-d1e29ce50ce2",
30+
"metadata": {
31+
"tags": []
32+
},
33+
"outputs": [],
34+
"source": [
35+
"asset = \"nsidc-s3\"\n",
36+
"resource = \"ATL03_20190314093716_11600203_005_01.h5\"\n",
37+
"region = icesat2.toregion(\"../tests/data/dicksonfjord.geojson\")\n",
38+
"parms = { \"poly\": region['poly'],\n",
39+
" \"cnf\": \"atl03_high\",\n",
40+
" \"ats\": 5.0,\n",
41+
" \"cnt\": 5,\n",
42+
" \"len\": 20.0,\n",
43+
" \"res\": 10.0,\n",
44+
" \"maxi\": 1,\n",
45+
" \"samples\": [\"arcticdem-mosaic\"] }\n",
46+
"gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])"
47+
]
48+
},
49+
{
50+
"cell_type": "code",
51+
"execution_count": null,
52+
"id": "d0ebba64-93ab-45c8-9c53-cc5078332617",
53+
"metadata": {},
54+
"outputs": [],
55+
"source": [
56+
"gdf[\"delta\"] = gdf[\"h_mean\"] - gdf[\"arcticdem-mosaic-1980-01-06\"]\n",
57+
"gdf[\"delta\"].describe()"
58+
]
59+
},
60+
{
61+
"cell_type": "code",
62+
"execution_count": null,
63+
"id": "db8b0b39-b421-46c6-bb5f-7abc6b1168b7",
64+
"metadata": {},
65+
"outputs": [],
66+
"source": [
67+
"# Setup Plot\n",
68+
"fig,ax = plt.subplots(num=None, figsize=(10, 8))\n",
69+
"fig.set_facecolor('white')\n",
70+
"fig.canvas.header_visible = False\n",
71+
"ax.set_title(\"SlideRule vs. ArcticDEM Elevations\")\n",
72+
"ax.set_xlabel('UTC')\n",
73+
"ax.set_ylabel('height (m)')\n",
74+
"legend_elements = []\n",
75+
"\n",
76+
"# Plot SlideRule ATL06 Elevations\n",
77+
"df = gdf[(gdf['rgt'] == 1160) & (gdf['gt'] == 10) & (gdf['cycle'] == 2)]\n",
78+
"sc1 = ax.scatter(df.index.values, df[\"h_mean\"].values, c='red', s=2.5)\n",
79+
"legend_elements.append(matplotlib.lines.Line2D([0], [0], color='red', lw=6, label='ATL06-SR'))\n",
80+
"\n",
81+
"# Plot ArcticDEM Elevations\n",
82+
"sc2 = ax.scatter(df.index.values, df[\"arcticdem-mosaic-1980-01-06\"].values, c='blue', s=2.5)\n",
83+
"legend_elements.append(matplotlib.lines.Line2D([0], [0], color='blue', lw=6, label='ArcticDEM'))\n",
84+
"\n",
85+
"# Display Legend\n",
86+
"lgd = ax.legend(handles=legend_elements, loc=3, frameon=True)\n",
87+
"lgd.get_frame().set_alpha(1.0)\n",
88+
"lgd.get_frame().set_edgecolor('white')\n",
89+
"\n",
90+
"# Show Plot\n",
91+
"plt.show()"
92+
]
93+
},
94+
{
95+
"cell_type": "code",
96+
"execution_count": null,
97+
"id": "82c65e28-468e-463e-9afe-2b52064e7bae",
98+
"metadata": {},
99+
"outputs": [],
100+
"source": [
101+
"# Setup Plot\n",
102+
"fig,ax = plt.subplots(num=None, figsize=(10, 8))\n",
103+
"fig.set_facecolor('white')\n",
104+
"fig.canvas.header_visible = False\n",
105+
"ax.set_title(\"Delta Elevations between SlideRule and ArcticDEM\")\n",
106+
"ax.set_xlabel('UTC')\n",
107+
"ax.set_ylabel('height (m)')\n",
108+
"ax.yaxis.grid(True)\n",
109+
"\n",
110+
"# Plot Deltas\n",
111+
"df = gdf[(gdf['rgt'] == 1160) & (gdf['gt'] == 10) & (gdf['cycle'] == 2)]\n",
112+
"sc1 = ax.scatter(df.index.values, df[\"delta\"].values, c='green', s=2.5)\n",
113+
"\n",
114+
"# Show Plot\n",
115+
"plt.show()"
116+
]
117+
},
118+
{
119+
"cell_type": "code",
120+
"execution_count": null,
121+
"id": "eed8f243-dd0c-4473-a952-fcb2bb863e3c",
122+
"metadata": {},
123+
"outputs": [],
124+
"source": []
125+
}
126+
],
127+
"metadata": {
128+
"kernelspec": {
129+
"display_name": "Python 3 (ipykernel)",
130+
"language": "python",
131+
"name": "python3"
132+
},
133+
"language_info": {
134+
"codemirror_mode": {
135+
"name": "ipython",
136+
"version": 3
137+
},
138+
"file_extension": ".py",
139+
"mimetype": "text/x-python",
140+
"name": "python",
141+
"nbconvert_exporter": "python",
142+
"pygments_lexer": "ipython3",
143+
"version": "3.8.15"
144+
}
145+
},
146+
"nbformat": 4,
147+
"nbformat_minor": 5
148+
}

‎requirements.txt

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,5 +2,5 @@ requests
22
numpy
33
fiona
44
geopandas
5-
shapely<2
5+
shapely
66
scikit-learn

‎sliderule/icesat2.py

Lines changed: 43 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,6 @@
4141
import geopandas
4242
from shapely.geometry.multipolygon import MultiPolygon
4343
from shapely.geometry import Polygon
44-
from sklearn.cluster import KMeans
4544
import sliderule
4645

4746
###############################################################################
@@ -51,6 +50,14 @@
5150
# create logger
5251
logger = logging.getLogger(__name__)
5352

53+
# import cluster support
54+
clustering_enabled = False
55+
try:
56+
from sklearn.cluster import KMeans
57+
clustering_enabled = True
58+
except:
59+
logger.warning("Unable to import sklearn... clustering support disabled")
60+
5461
# profiling times for each major function
5562
profiles = {}
5663

@@ -503,12 +510,8 @@ def __gdf2poly(gdf):
503510
#
504511
def __procoutputfile(parm, lon_key, lat_key):
505512
if "open_on_complete" in parm["output"] and parm["output"]["open_on_complete"]:
506-
# Read Parquet File as DataFrame
507-
df = geopandas.pd.read_parquet(parm["output"]["path"])
508-
# Build GeoDataFrame
509-
gdf = __todataframe(df, lon_key=lon_key, lat_key=lat_key)
510-
# Return Results
511-
return gdf
513+
# Return GeoParquet File as GeoDataFrame
514+
return geopandas.read_parquet(parm["output"]["path"])
512515
else:
513516
# Return Parquet Filename
514517
return parm["output"]["path"]
@@ -804,7 +807,7 @@ def atl06p(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, callb
804807
columns = {}
805808
elevation_records = []
806809
num_elevations = 0
807-
field_dictionary = {} # ['field_name'] = {"extent_id": [], field_name: []}
810+
field_dictionary = {} # [<field_name>] = {"extent_id": [], <field_name>: []}
808811
if len(rsps) > 0:
809812
# Sort Records
810813
for rsp in rsps:
@@ -814,15 +817,23 @@ def atl06p(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, callb
814817
elif 'extrec' == rsp['__rectype']:
815818
field_name = parm['atl03_geo_fields'][rsp['field_index']]
816819
if field_name not in field_dictionary:
817-
field_dictionary[field_name] = {"extent_id": [], field_name: []}
820+
field_dictionary[field_name] = {'extent_id': [], field_name: []}
818821
# Parse Ancillary Data
819-
data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
822+
data = __get_values(rsp['data'], rsp['datatype'], len(rsp['data']))
820823
# Add Left Pair Track Entry
821824
field_dictionary[field_name]['extent_id'] += rsp['extent_id'] | 0x2,
822825
field_dictionary[field_name][field_name] += data[LEFT_PAIR],
823826
# Add Right Pair Track Entry
824827
field_dictionary[field_name]['extent_id'] += rsp['extent_id'] | 0x3,
825828
field_dictionary[field_name][field_name] += data[RIGHT_PAIR],
829+
elif 'rsrec' == rsp['__rectype']:
830+
for sample in rsp["samples"]:
831+
time_str = sliderule.gps2utc(sample["time"])
832+
field_name = parm['samples'][rsp['raster_index']] + "-" + time_str.split(" ")[0].strip()
833+
if field_name not in field_dictionary:
834+
field_dictionary[field_name] = {'extent_id': [], field_name: []}
835+
field_dictionary[field_name]['extent_id'] += rsp['extent_id'],
836+
field_dictionary[field_name][field_name] += sample['value'],
826837
# Build Elevation Columns
827838
if num_elevations > 0:
828839
# Initialize Columns
@@ -967,23 +978,23 @@ def atl03sp(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, call
967978
# Get Field Type
968979
field_name = parm['atl03_geo_fields'][rsp['field_index']]
969980
if field_name not in extent_field_types:
970-
extent_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['data_type']]]["nptype"]
981+
extent_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['datatype']]]["nptype"]
971982
# Initialize Extent Dictionary Entry
972983
if extent_id not in extent_dictionary:
973984
extent_dictionary[extent_id] = {}
974985
# Save of Values per Extent ID per Field Name
975-
data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
986+
data = __get_values(rsp['data'], rsp['datatype'], len(rsp['data']))
976987
extent_dictionary[extent_id][field_name] = data
977988
elif 'phrec' == rsp['__rectype']:
978989
# Get Field Type
979990
field_name = parm['atl03_ph_fields'][rsp['field_index']]
980991
if field_name not in photon_field_types:
981-
photon_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['data_type']]]["nptype"]
992+
photon_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['datatype']]]["nptype"]
982993
# Initialize Extent Dictionary Entry
983994
if extent_id not in photon_dictionary:
984995
photon_dictionary[extent_id] = {}
985996
# Save of Values per Extent ID per Field Name
986-
data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
997+
data = __get_values(rsp['data'], rsp['datatype'], len(rsp['data']))
987998
photon_dictionary[extent_id][field_name] = data
988999
# Build Elevation Columns
9891000
if num_photons > 0:
@@ -1331,22 +1342,25 @@ def toregion(source, tolerance=0.0, cellsize=0.01, n_clusters=1):
13311342
# generate clusters
13321343
clusters = []
13331344
if n_clusters > 1:
1334-
# pull out centroids of each geometry object
1335-
if "CenLon" in gdf and "CenLat" in gdf:
1336-
X = numpy.column_stack((gdf["CenLon"], gdf["CenLat"]))
1345+
if clustering_enabled:
1346+
# pull out centroids of each geometry object
1347+
if "CenLon" in gdf and "CenLat" in gdf:
1348+
X = numpy.column_stack((gdf["CenLon"], gdf["CenLat"]))
1349+
else:
1350+
s = gdf.centroid
1351+
X = numpy.column_stack((s.x, s.y))
1352+
# run k means clustering algorithm against polygons in gdf
1353+
kmeans = KMeans(n_clusters=n_clusters, init='k-means++', random_state=5, max_iter=400)
1354+
y_kmeans = kmeans.fit_predict(X)
1355+
k = geopandas.pd.DataFrame(y_kmeans, columns=['cluster'])
1356+
gdf = gdf.join(k)
1357+
# build polygon for each cluster
1358+
for n in range(n_clusters):
1359+
c_gdf = gdf[gdf["cluster"] == n]
1360+
c_poly = __gdf2poly(c_gdf)
1361+
clusters.append(c_poly)
13371362
else:
1338-
s = gdf.centroid
1339-
X = numpy.column_stack((s.x, s.y))
1340-
# run k means clustering algorithm against polygons in gdf
1341-
kmeans = KMeans(n_clusters=n_clusters, init='k-means++', random_state=5, max_iter=400)
1342-
y_kmeans = kmeans.fit_predict(X)
1343-
k = geopandas.pd.DataFrame(y_kmeans, columns=['cluster'])
1344-
gdf = gdf.join(k)
1345-
# build polygon for each cluster
1346-
for n in range(n_clusters):
1347-
c_gdf = gdf[gdf["cluster"] == n]
1348-
c_poly = __gdf2poly(c_gdf)
1349-
clusters.append(c_poly)
1363+
raise sliderule.FatalError("Clustering support not enabled; unable to import sklearn package")
13501364

13511365
# update timing profiles
13521366
profiles[toregion.__name__] = time.perf_counter() - tstart

0 commit comments

Comments
 (0)