#!/usr/bin/env python
# coding: utf-8

# # 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)


# 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)


# In[5]:


gf = items_s2.to_geopandas()
gf.head()


# In[6]:


get_ipython().run_cell_magic('time', '', "item = items_s2['S2A_10TEL_20201130_0_L2A']\nda = item.B04.to_dask()\nda\n")


# In[7]:


get_ipython().run_cell_magic('time', '', "import hvplot.xarray\nsubset = da.isel(x=slice(6000,10980), y=slice(0,3000))\ntitle = item.name + ' - B04'\nsubset.hvplot.image(rasterize=True, aspect='equal', cmap='magma', frame_width=500, title=title)\n")


# 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]:


get_ipython().run_cell_magic('time', '', "\nitem = items_s1['G1714627010-ASF']\nxarray_opts = dict(group='/science/grids/data', engine='h5netcdf')\nds = item.data(xarray_kwargs=xarray_opts).to_dask()\n\nsubset = ds['amplitude'].sel(latitude=slice(41.53, 41.24), \n                             longitude=slice(-122.41, -121.86))\n\nsubset.hvplot.image(x='longitude', y='latitude', \n                             rasterize=True, frame_width=500,\n                             cmap='gray', clim=(1,1e5), logz=True)\n")