Medium Blog post GIST

In [1]:
import os
os.environ['GDAL_DISABLE_READDIR_ON_OPEN']='EMPTY_DIR'
In [2]:
import intake
from satsearch import Search

def get_STAC_items(url, collection, dates, bbox):
    results = Search.search(url=url,
                        collections=[collection], 
                        datetime=dates,
                        bbox=bbox,    
                        sortby=['-properties.datetime'])

    items = results.items()
    print(f'Found {len(items)} Items')
    
    return intake.open_stac_item_collection(items)
In [3]:
url='https://cmr.earthdata.nasa.gov/stac/ASF'
bbox = [-122.4, 41.3, -122.1, 41.5] # (min lon, min lat, max lon, max lat)
collection = 'C1595422627-ASF' # sentinel-1 beta interferograms
dates = '2020-01-01/2020-12-31' 
items_s1 = get_STAC_items(url, collection, dates, bbox)

Found 46 Items
In [4]:
url='https://earth-search.aws.element84.com/v0'
collection = 'sentinel-s2-l2a-cogs' # sentinel-2 w/ atmospheric corrections
items_s2 = get_STAC_items(url, collection, dates, bbox)

Found 275 Items
In [5]:
gf = items_s2.to_geopandas()
gf.head()

/Users/scott/miniconda3/envs/intake-stac-gui/lib/python3.7/site-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))
Out[5]:
geometry datetime platform constellation instruments gsd view:off_nadir proj:epsg sentinel:utm_zone sentinel:latitude_band sentinel:grid_square sentinel:sequence sentinel:product_id sentinel:data_coverage eo:cloud_cover sentinel:valid_cloud_cover created updated data_coverage
0 POLYGON ((-121.70337 40.55548, -122.71338 40.5... 2020-12-05T19:03:09Z sentinel-2b sentinel-2 [msi] 10 0 32610 10 T EL 0 S2B_MSIL2A_20201205T185749_N0214_R113_T10TEL_2... 66.35 17.15 True 2020-12-06T00:12:45.530Z 2020-12-06T00:12:45.530Z NaN
1 POLYGON ((-121.68559 41.45625, -122.43644 41.4... 2020-12-05T19:02:54Z sentinel-2b sentinel-2 [msi] 10 0 32610 10 T EM 0 S2B_MSIL2A_20201205T185749_N0214_R113_T10TEM_2... 45.77 0.00 True 2020-12-06T00:10:16.786Z 2020-12-06T00:10:16.786Z NaN
2 POLYGON ((-121.70337 40.55548, -123.00022 40.5... 2020-12-03T19:13:06Z sentinel-2a sentinel-2 [msi] 10 0 32610 10 T EL 0 S2A_MSIL2A_20201203T190751_N0214_R013_T10TEL_2... 100.00 0.26 True 2020-12-04T00:28:15.281Z 2020-12-04T00:28:15.281Z NaN
3 POLYGON ((-121.68559 41.45625, -123.00023 41.4... 2020-12-03T19:12:52Z sentinel-2a sentinel-2 [msi] 10 0 32610 10 T EM 0 S2A_MSIL2A_20201203T190751_N0214_R013_T10TEM_2... 100.00 20.43 True 2020-12-04T00:22:22.583Z 2020-12-04T00:22:22.583Z NaN
4 POLYGON ((-121.70337 40.55548, -122.71551 40.5... 2020-11-30T19:03:11Z sentinel-2a sentinel-2 [msi] 10 0 32610 10 T EL 0 S2A_MSIL2A_20201130T185741_N0214_R113_T10TEL_2... 66.50 2.05 True 2020-11-30T23:41:45.912Z 2020-11-30T23:41:45.912Z NaN
In [6]:
%%time
item = items_s2['S2A_10TEL_20201130_0_L2A']
da = item.B04.to_dask()
da

CPU times: user 386 ms, sys: 92.2 ms, total: 478 ms
Wall time: 955 ms
Out[6]:
<xarray.DataArray (band: 1, y: 10980, x: 10980)>
dask.array<open_rasterio-99e229fd7b208acb72b593a6303cbf08<this-array>, shape=(1, 10980, 10980), dtype=uint16, chunksize=(1, 10980, 10980), chunktype=numpy.ndarray>
Coordinates:
  * band     (band) int64 1
  * y        (y) float64 4.6e+06 4.6e+06 4.6e+06 ... 4.49e+06 4.49e+06 4.49e+06
  * x        (x) float64 5e+05 5e+05 5e+05 ... 6.098e+05 6.098e+05 6.098e+05
Attributes:
    transform:           (10.0, 0.0, 499980.0, 0.0, -10.0, 4600020.0)
    crs:                 +init=epsg:32610
    res:                 (10.0, 10.0)
    is_tiled:            1
    nodatavals:          (0.0,)
    scales:              (1.0,)
    offsets:             (0.0,)
    AREA_OR_POINT:       Area
    OVR_RESAMPLING_ALG:  AVERAGE
In [7]:
%%time 
import hvplot.xarray
subset = da.isel(x=slice(6000,10980), y=slice(0,3000))
title = item.name + ' - B04'
subset.hvplot.image(rasterize=True, aspect='equal', cmap='magma', frame_width=500, title=title)