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	proj:epsg	sentinel:utm_zone	sentinel:latitude_band	sentinel:grid_square	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	32610	10	T	EL	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	32610	10	T	EM	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	32610	10	T	EL	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	32610	10	T	EM	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	32610	10	T	EL	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

xarray.DataArray

band: 1
y: 10980
x: 10980

dask.array<chunksize=(1, 10980, 10980), meta=np.ndarray>

	Array	Chunk
Bytes	241.12 MB	241.12 MB
Shape	(1, 10980, 10980)	(1, 10980, 10980)
Count	2 Tasks	1 Chunks
Type	uint16	numpy.ndarray

Coordinates: (3)

band
(band)
int64
1
```
array([1])
```

y

(y)

float64

4.6e+06 4.6e+06 ... 4.49e+06

array([4600015., 4600005., 4599995., ..., 4490245., 4490235., 4490225.])

x

(x)

float64

5e+05 5e+05 ... 6.098e+05 6.098e+05

array([499985., 499995., 500005., ..., 609755., 609765., 609775.])

Attributes: (9)
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)

CPU times: user 2.38 s, sys: 342 ms, total: 2.72 s
Wall time: 3.29 s

Out[7]:

In [8]:

# NASA URS Authentication
import netrc, fsspec, aiohttp

(username, account, password) = netrc.netrc().authenticators("urs.earthdata.nasa.gov")
fsspec.config.conf['https'] = dict(client_kwargs={'auth': aiohttp.BasicAuth(username, password)})

In [9]:

%%time

item = items_s1['G1714627010-ASF']
xarray_opts = dict(group='/science/grids/data', engine='h5netcdf')
ds = item.data(xarray_kwargs=xarray_opts).to_dask()

subset = ds['amplitude'].sel(latitude=slice(41.53, 41.24), 
                             longitude=slice(-122.41, -121.86))

subset.hvplot.image(x='longitude', y='latitude', 
                             rasterize=True, frame_width=500,
                             cmap='gray', clim=(1,1e5), logz=True)

/Users/scott/GitHub/intake-stac/intake_stac/catalog.py:452: UserWarning: STAC Asset "type" missing, assigning None based on href suffix .nc:
{'href': 'https://grfn.asf.alaska.edu/door/download/S1-GUNW-A-R-137-tops-20200311_20200305-015951-41897N_400275N-PP-8413-v2_0_2.nc', 'type': 'application/netcdf'}
  f'STAC Asset "type" missing, assigning {entry_type} based on href suffix {suffix}:\n{entry}'  # noqa: E501

application/netcdf
CPU times: user 1.5 s, sys: 1.52 s, total: 3.02 s
Wall time: 11.1 s

Out[9]: