Merge branch 'development' of github.com:ICESat2-SlideRule/sliderule-python into development

Author: elidwa

examples/arcticdem_strip_boundaries.ipynb

@@ -34,7 +34,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "8ac0d0cc-a65d-4368-baa9-7d5dfb4936a7",
+   "id": "78db0877-2ba9-4e8f-ba2e-a07e4de491ca",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -71,7 +71,6 @@
     "          \"len\": 40.0,\n",
     "          \"res\": 120.0,\n",
     "          \"maxi\": 5,\n",
-    "          \"cycle\": 15,\n",
     "          \"rgt\": 658,\n",
     "          \"time_start\":'2020-01-01',\n",
     "          \"time_end\":'2021-01-01',\n",
@@ -82,14 +81,15 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "970e0770-abe9-4e2a-b1af-bc6e330fb566",
-   "metadata": {},
+   "id": "b815e8d8-c71a-4305-8649-fa4e95b64325",
+   "metadata": {
+    "tags": []
+   },
    "outputs": [],
    "source": [
     "#\n",
     "# Set DEM of Interest\n",
     "#\n",
-    "id = 17\n",
     "gdf.attrs['file_directory']"
    ]
   },
@@ -143,8 +143,11 @@
     "#\n",
     "# Get Boundaries for each Raster\n",
     "#\n",
+    "p0 = 132\n",
+    "p1 = 148\n",
     "raster_of_interest = {}\n",
-    "for i in gdf.attrs['file_directory']:\n",
+    "for i in list(gdf.attrs['file_directory'].keys())[p0:p1]:\n",
+    "    print(\"Retrieving raster info for:\", gdf.attrs['file_directory'][i])\n",
     "    rlist = getBB(gdf.attrs['file_directory'][i])\n",
     "    raster_of_interest[\"dem\"+str(i)] = sliderule.toregion(rlist)"
    ]
@@ -168,39 +171,36 @@
     "    else:\n",
     "        return None\n",
     "sampled_data = gdf[gdf['strips.time'].notnull()]\n",
-    "for i in gdf.attrs['file_directory']:\n",
+    "for i in list(gdf.attrs['file_directory'].keys())[p0:p1]:\n",
     "    sampled_data[\"dem\"+str(i)] = sampled_data.apply(lambda x: getValue(x, i), axis=1)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "24e1c0c3-b4e5-4bcd-af34-26bfa6885c91",
+   "id": "f23a3827-93b6-47cb-9842-8f8269e0a871",
    "metadata": {},
    "outputs": [],
    "source": [
-    "\n",
-    "#\n",
-    "# Create Boundaries for Region and Raster\n",
-    "#\n",
-    "\n",
     "#\n",
-    "# Plot\n",
+    "# Plot Overlays of Boundaries and Returns\n",
     "#\n",
     "fig = plt.figure(num=None, figsize=(24, 24))\n",
     "region_lons = [p[\"lon\"] for p in region_of_interest[\"poly\"]]\n",
     "region_lats = [p[\"lat\"] for p in region_of_interest[\"poly\"]]\n",
     "ax = {}\n",
-    "for i in gdf.attrs['file_directory']:\n",
+    "k = 0\n",
+    "for i in list(gdf.attrs['file_directory'].keys())[p0:p1]:\n",
     "    raster_lons = [p[\"lon\"] for p in raster_of_interest[\"dem\"+str(i)][\"poly\"]]\n",
     "    raster_lats = [p[\"lat\"] for p in raster_of_interest[\"dem\"+str(i)][\"poly\"]]\n",
     "    plot_data = sampled_data[sampled_data[\"dem\"+str(i)].notnull()]\n",
     "    plot_data = sampled_data[sampled_data[\"dem\"+str(i)] > -9990]\n",
-    "    ax[i] = plt.subplot(5,4,i+1)\n",
-    "    gdf.plot(ax=ax[i], column='h_mean', color='y', markersize=0.5)\n",
-    "    plot_data.plot(ax=ax[i], column='h_mean', color='b', markersize=0.5)\n",
-    "    ax[i].plot(region_lons, region_lats, linewidth=1.5, color='r', zorder=2)\n",
-    "    ax[i].plot(raster_lons, raster_lats, linewidth=1.5, color='g', zorder=2)\n",
+    "    ax[k] = plt.subplot(5,4,k+1)\n",
+    "    gdf.plot(ax=ax[k], column='h_mean', color='y', markersize=0.5)\n",
+    "    plot_data.plot(ax=ax[k], column='h_mean', color='b', markersize=0.5)\n",
+    "    ax[k].plot(region_lons, region_lats, linewidth=1.5, color='r', zorder=2)\n",
+    "    ax[k].plot(raster_lons, raster_lats, linewidth=1.5, color='g', zorder=2)\n",
+    "    k += 1\n",
     "plt.tight_layout()"
    ]
   },
@@ -218,7 +218,7 @@
     "######################\n",
     "# Select DEM File ID #\n",
     "######################\n",
-    "file_id = 5\n",
+    "file_id = list(gdf.attrs['file_directory'].keys())[p0]\n",
     "\n",
     "# Setup Plot\n",
     "fig,ax = plt.subplots(num=None, figsize=(10, 8))\n",

sliderule/gedi.py

@@ -51,9 +51,6 @@
 # default GEDI standard data product version
 DEFAULT_GEDI_SDP_VERSION = '2'
 
-# gps-based epoch for delta times
-GEDI_SDP_EPOCH = datetime.datetime(2018, 1, 1)
-
 # gedi parameters
 ALL_BEAMS = -1
 
@@ -64,7 +61,7 @@
 #
 #  Dictionary to GeoDataFrame
 #
-def __todataframe(columns, delta_time_key="delta_time", lon_key="longitude", lat_key="latitude", **kwargs):
+def __todataframe(columns, time_key="time", lon_key="longitude", lat_key="latitude", **kwargs):
 
     # Latch Start Time
     tstart = time.perf_counter()
@@ -78,9 +75,7 @@ def __todataframe(columns, delta_time_key="delta_time", lon_key="longitude", lat
         return sliderule.emptyframe(**kwargs)
 
     # Generate Time Column
-    delta_time = (columns[delta_time_key]*1e9).astype('timedelta64[ns]')
-    gedi_sdp_epoch = numpy.datetime64(GEDI_SDP_EPOCH)
-    columns['time'] = geopandas.pd.to_datetime(gedi_sdp_epoch + delta_time)
+    columns['time'] = columns[time_key].astype('datetime64[ns]')
 
     # Generate Geometry Column
     geometry = geopandas.points_from_xy(columns[lon_key], columns[lat_key])

sliderule/icesat2.py

@@ -88,9 +88,6 @@
 P = { '5':   0, '10':  1, '15':  2, '20':  3, '25':  4, '30':  5, '35':  6, '40':  7, '45':  8, '50': 9,
       '55': 10, '60': 11, '65': 12, '70': 13, '75': 14, '80': 15, '85': 16, '90': 17, '95': 18 }
 
-# gps-based epoch for delta times
-ATLAS_SDP_EPOCH = datetime.datetime(2018, 1, 1)
-
 ###############################################################################
 # LOCAL FUNCTIONS
 ###############################################################################
@@ -142,7 +139,7 @@ def __calcspot(sc_orient, track, pair):
 #
 #  Dictionary to GeoDataFrame
 #
-def __todataframe(columns, delta_time_key="delta_time", lon_key="lon", lat_key="lat", **kwargs):
+def __todataframe(columns, time_key="time", lon_key="lon", lat_key="lat", **kwargs):
 
     # Latch Start Time
     tstart = time.perf_counter()
@@ -156,9 +153,7 @@ def __todataframe(columns, delta_time_key="delta_time", lon_key="lon", lat_key="
         return sliderule.emptyframe(**kwargs)
 
     # Generate Time Column
-    delta_time = (columns[delta_time_key]*1e9).astype('timedelta64[ns]')
-    atlas_sdp_epoch = numpy.datetime64(ATLAS_SDP_EPOCH)
-    columns['time'] = geopandas.pd.to_datetime(atlas_sdp_epoch + delta_time)
+    columns['time'] = columns[time_key].astype('datetime64[ns]')
 
     # Generate Geometry Column
     geometry = geopandas.points_from_xy(columns[lon_key], columns[lat_key])

sliderule/sliderule.py

@@ -119,7 +119,7 @@
     "BITFIELD": { "fmt": 'x', "size": 0, "nptype": numpy.byte   },  # unsupported
     "FLOAT":    { "fmt": 'f', "size": 4, "nptype": numpy.single },
     "DOUBLE":   { "fmt": 'd', "size": 8, "nptype": numpy.double },
-    "TIME8":    { "fmt": 'Q', "size": 8, "nptype": numpy.byte   },
+    "TIME8":    { "fmt": 'q', "size": 8, "nptype": numpy.int64  }, # numpy.datetime64
     "STRING":   { "fmt": 's', "size": 1, "nptype": numpy.byte   }
 }
 
@@ -917,7 +917,7 @@ def authenticate (ps_organization, ps_username=None, ps_password=None):
 #
 # gps2utc
 #
-def gps2utc (gps_time, as_str=True, epoch=gps_epoch):
+def gps2utc (gps_time, as_str=True):
     '''
     Convert a GPS based time returned from SlideRule into a UTC time.
 
@@ -927,8 +927,6 @@ def gps2utc (gps_time, as_str=True, epoch=gps_epoch):
                     number of seconds since GPS epoch (January 6, 1980)
         as_str:     bool
                     if True, returns the time as a string; if False, returns the time as datatime object
-        epoch:      datetime
-                    the epoch used in the conversion, defaults to GPS epoch (Jan 6, 1980)
 
     Returns
     -------
@@ -941,15 +939,11 @@ def gps2utc (gps_time, as_str=True, epoch=gps_epoch):
         >>> sliderule.gps2utc(1235331234)
         '2019-02-27 19:34:03'
     '''
-    gps_time = epoch + timedelta(seconds=gps_time)
-    tai_time = gps_time + timedelta(seconds=19)
-    tai_timestamp = (tai_time - tai_epoch).total_seconds()
-    utc_timestamp = datetime.utcfromtimestamp(tai_timestamp)
+    rsps = source("time", {"time": gps_time, "input": "GPS", "output": "DATE"})
     if as_str:
-        return str(utc_timestamp)
+        return rsps["time"]
     else:
-        return utc_timestamp
-
+        return datetime.strptime(rsps["time"], '%Y-%m-%dT%H:%M:%SZ')
 #
 # get_definition
 #

tests/test_api.py

@@ -118,7 +118,7 @@ def test_definition(self, domain, organization):
             "rectype": "atl06rec.elevation",
         }
         d = sliderule.source("definition", rqst)
-        assert d["delta_time"]["offset"] == 256
+        assert d["time"]["offset"] == 192
 
     def test_version(self, domain, organization):
         icesat2.init(domain, organization=organization)

tests/test_arcticdem.py

@@ -32,7 +32,7 @@ def test_nearestneighbour(self, domain, asset, organization):
                   "samples": {"mosaic": {"asset": "arcticdem-mosaic"}} }
         gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])
         assert len(gdf) == 964
-        assert len(gdf.keys()) == 20
+        assert len(gdf.keys()) == 19
         assert gdf["rgt"][0] == 1160
         assert gdf["cycle"][0] == 2
         assert gdf['segment_id'].describe()["min"] == 405240
@@ -54,7 +54,7 @@ def test_zonal_stats(self, domain, asset, organization):
                   "samples": {"mosaic": {"asset": "arcticdem-mosaic", "radius": 10.0, "zonal_stats": True}} }
         gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])
         assert len(gdf) == 964
-        assert len(gdf.keys()) == 27
+        assert len(gdf.keys()) == 26
         assert gdf["rgt"][0] == 1160
         assert gdf["cycle"][0] == 2
         assert gdf['segment_id'].describe()["min"] == 405240

tests/test_client.py

@@ -13,9 +13,10 @@ def test_seturl_empty(self):
         with pytest.raises(TypeError, match=('url')):
             sliderule.set_url()
 
-    def test_gps2utc(self):
-        utc = sliderule.gps2utc(1235331234)
-        assert utc == '2019-02-27 19:34:03'
+    def test_gps2utc(self, domain, organization):
+        sliderule.init(domain, organization=organization)
+        utc = sliderule.gps2utc(1235331234000)
+        assert utc == '2019-02-27T19:33:36Z'
 
 @pytest.mark.network
 class TestRemote:

tests/test_geojson.py

@@ -10,7 +10,7 @@
 
 @pytest.mark.network
 class TestGeoJson:
-    def test_atl06(self, domain, asset, organization):
+    def test_atl03(self, domain, asset, organization):
         icesat2.init(domain, organization=organization)
         for testfile in ["data/grandmesa.geojson", "data/grandmesa.shp"]:
             region = sliderule.toregion(os.path.join(TESTDIR, testfile))

tests/test_parquet.py

@@ -2,6 +2,7 @@
 
 import pytest
 from pathlib import Path
+import numpy
 import os
 import os.path
 import sliderule
@@ -26,7 +27,7 @@ def test_atl06(self, domain, asset, organization):
                   "output": { "path": "testfile1.parquet", "format": "parquet", "open_on_complete": True } }
         gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])
         assert len(gdf) == 964
-        assert len(gdf.keys()) == 17
+        assert len(gdf.keys()) == 16
         assert gdf["rgt"][0] == 1160
         assert gdf["cycle"][0] == 2
         assert gdf['segment_id'].describe()["min"] == 405240
@@ -48,9 +49,49 @@ def test_atl03(self, domain, asset, organization):
                   "output": { "path": "testfile2.parquet", "format": "parquet", "open_on_complete": True } }
         gdf = icesat2.atl03sp(parms, asset=asset, resources=[resource])
         assert len(gdf) == 194696
-        assert len(gdf.keys()) == 18
+        assert len(gdf.keys()) == 17
         assert gdf["rgt"][0] == 1160
         assert gdf["cycle"][0] == 2
         assert gdf['segment_id'].describe()["min"] == 405240
         assert gdf['segment_id'].describe()["max"] == 405915
         os.remove("testfile2.parquet")
+
+    def test_atl06_index(self, domain, asset, organization):
+        icesat2.init(domain, organization=organization)
+        resource = "ATL03_20181017222812_02950102_005_01.h5"
+        region = sliderule.toregion(os.path.join(TESTDIR, "data/grandmesa.geojson"))
+        parms = {
+            "poly": region["poly"],
+            "raster": region["raster"],
+            "srt": icesat2.SRT_LAND,
+            "cnf": icesat2.CNF_SURFACE_HIGH,
+            "ats": 10.0,
+            "cnt": 10,
+            "len": 40.0,
+            "res": 20.0,
+            "output": { "path": "testfile3.parquet", "format": "parquet", "open_on_complete": True } }
+        gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])
+        assert len(gdf) == 275
+        assert gdf.index.values.min() == numpy.datetime64('2018-10-17T22:31:35.330187264')
+        assert gdf.index.values.max() == numpy.datetime64('2018-10-17T22:31:37.693747456')
+        os.remove("testfile3.parquet")
+
+    def test_atl03_index(self, domain, asset, organization):
+        icesat2.init(domain, organization=organization)
+        resource = "ATL03_20181017222812_02950102_005_01.h5"
+        region = sliderule.toregion(os.path.join(TESTDIR, "data/grandmesa.geojson"))
+        parms = {
+            "poly": region["poly"],
+            "raster": region["raster"],
+            "srt": icesat2.SRT_LAND,
+            "cnf": icesat2.CNF_SURFACE_HIGH,
+            "ats": 10.0,
+            "cnt": 10,
+            "len": 40.0,
+            "res": 20.0,
+            "output": { "path": "testfile4.parquet", "format": "parquet", "open_on_complete": True } }
+        gdf = icesat2.atl03sp(parms, asset=asset, resources=[resource])
+        assert len(gdf) == 21029
+        assert gdf.index.values.min() == numpy.datetime64('2018-10-17T22:31:35.330047488')
+        assert gdf.index.values.max() == numpy.datetime64('2018-10-17T22:31:37.695347456')
+        os.remove("testfile4.parquet")