Merge branch 'development' into main

sliderule/icesat2.py

@@ -155,6 +155,10 @@ def __todataframe(columns, time_key="time", lon_key="lon", lat_key="lat", **kwar
     # Generate Time Column
     columns['time'] = columns[time_key].astype('datetime64[ns]')
 
+    # Temporary code for backward compatibility
+    if 'delta_time' in columns:
+        del columns['delta_time']
+
     # Generate Geometry Column
     geometry = geopandas.points_from_xy(columns[lon_key], columns[lat_key])
     del columns[lon_key]

sliderule/sliderule.py

@@ -1066,7 +1066,7 @@ def toregion(source, tolerance=0.0, cellsize=0.01, n_clusters=1):
     dict
         a list of longitudes and latitudes containing the region of interest that can be used for the **poly** and **raster** parameters in a processing request to SlideRule.
 
-        region = {"poly": [{"lat": <lat1>, "lon": <lon1>, ... }], "clusters": [{"lat": <lat1>, "lon": <lon1>, ... }, {"lat": <lat1>, "lon": <lon1>, ... }, ...], "raster": {"data": <geojson file as string>, "length": <length of geojson file>, "cellsize": <parameter cellsize>}}
+        region = {"poly": [{"lat": <lat1>, "lon": <lon1> }, ...], "clusters": [[{"lat": <lat1>, "lon": <lon1>}, ...], [{"lat": <lat1>, "lon": <lon1>}, ...]], "raster": {"data": <geojson file as string>, "length": <length of geojson file>, "cellsize": <parameter cellsize>}}
 
     Examples
     --------

tests/test_arcticdem.py

@@ -9,7 +9,7 @@
 TESTDIR = Path(__file__).parent
 
 @pytest.mark.network
-class TestVrt:
+class TestMosaic:
     def test_vrt(self, domain, organization):
         icesat2.init(domain, organization=organization)
         rqst = {"samples": {"asset": "arcticdem-mosaic"}, "coordinates": [[-178.0,51.7]]}

tests/test_landsat.py

@@ -0,0 +1,41 @@
+"""Tests for sliderule-python landsat raster support."""
+
+import pytest
+from pathlib import Path
+import os.path
+from sliderule import sliderule, earthdata, icesat2
+
+TESTDIR = Path(__file__).parent
+
+@pytest.mark.network
+class TestHLS:
+    def test_samples(self, domain, organization):
+        sliderule.init(domain, organization=organization)
+        time_start = "2021-01-01T00:00:00Z"
+        time_end = "2021-02-01T23:59:59Z"
+        polygon = [ {"lon": -177.0000000001, "lat": 51.0000000001},
+                    {"lon": -179.0000000001, "lat": 51.0000000001},
+                    {"lon": -179.0000000001, "lat": 49.0000000001},
+                    {"lon": -177.0000000001, "lat": 49.0000000001},
+                    {"lon": -177.0000000001, "lat": 51.0000000001} ]
+        catalog = earthdata.stac(short_name="HLS", polygon=polygon, time_start=time_start, time_end=time_end, as_str=True)
+        rqst = {"samples": {"asset": "landsat-hls", "catalog": catalog, "bands": ["B02"]}, "coordinates": [[-178.0, 50.7]]}
+        rsps = sliderule.source("samples", rqst)
+
+    def test_ndvi(self, domain, asset, organization):
+        icesat2.init(domain, organization=organization)
+        region = sliderule.toregion(os.path.join(TESTDIR, "data/grandmesa.geojson"))
+        resource = "ATL03_20181017222812_02950102_005_01.h5"
+        time_start = "2021-01-01T00:00:00Z"
+        time_end = "2021-02-01T23:59:59Z"
+        catalog = earthdata.stac(short_name="HLS", polygon=region['poly'], time_start=time_start, time_end=time_end, as_str=True)
+        parms = { "poly": region['poly'],
+                  "raster": region['raster'],
+                  "cnf": "atl03_high",
+                  "ats": 20.0,
+                  "cnt": 10,
+                  "len": 40.0,
+                  "res": 20.0,
+                  "maxi": 1,
+                  "samples": {"ndvi": {"asset": "landsat-hls", "catalog": catalog, "bands": ["NDVI"]}} }
+        gdf = icesat2.atl06p(parms, asset=asset, resources=[resource])

utils/landsat.py

@@ -38,7 +38,7 @@ def get_temp_creds(provider):
 
     catalog = Client.open("https://cmr.earthdata.nasa.gov/stac/LPCLOUD")
 
-    timeRange = '2021-01-01/2022-01-01'
+    timeRange = '2021-01-01/2021-02-01'
     geometry  = BuildSquare(-178, 50, 1)
     mybbox    = [-179,41,-177,51]
 
@@ -48,10 +48,12 @@ def get_temp_creds(provider):
     #                          datetime=timeRange)
     # print(f"Results matched: {results.matched()}")
 
+
     print("Searching with Intersects...")
     results = catalog.search(collections=['HLSS30.v2.0', 'HLSL30.v2.0'],
                              datetime=timeRange,
-                             intersects=geometry)
+                             intersects=geometry,
+                             fields=['HORIZONTAL_CS_CODE', 'HORIZONTAL_CS_NAME'])
 
     # print(f"{results.url_with_parameters()}")
     print(f"Results matched: {results.matched()}")
@@ -66,17 +68,24 @@ def get_temp_creds(provider):
 
     urlList = []
 
+
     for i in reversed(range(len(itemsDict["features"]))):
         del itemsDict["features"][i]["links"]
         del itemsDict["features"][i]["stac_version"]
         del itemsDict["features"][i]["stac_extensions"]
         del itemsDict["features"][i]["collection"]
         del itemsDict["features"][i]["bbox"]
         del itemsDict["features"][i]["assets"]["browse"]
-        del itemsDict["features"][i]["assets"]["metadata"]
 
         propertiesDict = itemsDict["features"][i]["properties"]
         assetsDict     = itemsDict["features"][i]["assets"]
+        metaFileUrl    = assetsDict["metadata"]["href"]
+        # we may have to read metaFile, get some algo related attributes from it
+        # and add them to properties ie:
+        # propertiesDict['algo1'] = someAlgo1value
+
+        del itemsDict["features"][i]["assets"]["metadata"]
+
         for val in assetsDict:
             if "href" in assetsDict[val]:
                 # add raster url to properties
@@ -88,7 +97,7 @@ def get_temp_creds(provider):
 
 
 
-    # Dumped trimmed dictionary as geojson file, for testing
+    # Dump trimmed dictionary as geojson file, for testing
     file = 'hls_trimmed.geojson'
     print(f"Writing reults to file {file}")
     with open(file, 'w') as fp:

utils/query_cmr.py

@@ -18,7 +18,6 @@
         "tolerance": 0.0,
         "dataset": "ATL03",
         "version": "005",
-        "provider": "NSIDC_ECS"
     }
 
     # Command line parameters
@@ -28,7 +27,7 @@
     region = sliderule.toregion(cfg["region"], cfg["tolerance"])
 
     # Query CMR for list of resources
-    resources = earthdata.cmr(provider=cfg["provider"], short_name=cfg["dataset"], version=cfg["version"], polygon=region["poly"])
+    resources = earthdata.cmr(short_name=cfg["dataset"], version=cfg["version"], polygon=region["poly"])
     print("Region: {} points, {} files".format(len(region["poly"]), len(resources)))
     for resource in resources:
         print(resource)

utils/query_stac.py

@@ -0,0 +1,38 @@
+#
+# Imports
+#
+import sys
+import json
+import sliderule
+from utils import parse_command_line
+from sliderule import earthdata
+
+###############################################################################
+# MAIN
+###############################################################################
+
+if __name__ == '__main__':
+
+    # Defaults
+    cfg = {
+        "dataset": "HLS",
+        "region": "examples/grandmesa.geojson",
+        "time_start": "2021-01-01T00:00:00Z",
+        "time_end": "2022-01-01T23:59:59Z",
+        "display_all": False
+    }
+
+    # Command line parameters
+    parse_command_line(sys.argv, cfg)
+
+    # Region of interest
+    polygon = sliderule.toregion(cfg["region"])["poly"]
+
+    # Query CMR for list of resources
+    geojson = earthdata.stac(short_name=cfg["dataset"], polygon=polygon, time_start=cfg["time_start"], time_end=cfg["time_end"])
+
+    # Display Results
+    print("Returned {} features".format(geojson["context"]["returned"]))
+    if cfg["display_all"]:
+        print(json.dumps(geojson, indent=2))
+

utils/stac.py

@@ -0,0 +1,61 @@
+from pystac_client import Client
+import sliderule
+import requests
+import json
+import sys
+
+# parameters
+debug = False
+client = False
+manual = True
+page = None
+url = 'https://cmr.earthdata.nasa.gov/stac/LPCLOUD'
+headers = {'Content-Type': 'application/json'}
+limit = 200
+collections = ["HLSS30.v2.0", "HLSL30.v2.0"]
+time_start = "2021-01-01T00:00:00Z"
+time_end = "2022-01-01T23:59:59Z"
+coordinates = [[[coord["lon"], coord["lat"]] for coord in polygon] for polygon in [sliderule.toregion("examples/grandmesa.geojson")["poly"]]]
+
+# debug output
+if debug:
+    import logging
+    import http.client as http_client
+    http_client.HTTPConnection.debuglevel = 1
+    logging.basicConfig()
+    logging.getLogger().setLevel(logging.DEBUG)
+    requests_log = logging.getLogger("requests.packages.urllib3")
+    requests_log.setLevel(logging.DEBUG)
+    requests_log.propagate = True
+
+
+# client stac request
+if client:
+    catalog = Client.open(url)
+    timeRange = '{}/{}'.format(time_start, time_end)
+    geometry = {"type": "Polygon", "coordinates": coordinates}
+    results = catalog.search(collections=collections, datetime=timeRange, intersects=geometry)
+    items = results.get_all_items_as_dict()
+    print(f"Returned {results.matched()} features")
+
+
+# manual stac request
+if manual:
+    rqst = {
+        "limit": 100,
+        "datetime": '{}/{}'.format(time_start, time_end),
+        "collections": collections,
+        "intersects": {
+            "type": "Polygon",
+            "coordinates": coordinates
+        }
+    }
+    if page:
+        rqst["page"] = page
+    context = requests.Session()
+    context.trust_env = False
+    data = context.post(url+"/search", data=json.dumps(rqst), headers=headers)
+    data.raise_for_status()
+    geojson = data.json()
+    print("Returned {} features".format(geojson["context"]["returned"]))
+