Create buffer around point location

import geopandas
import json
import shapely
import shapely.geometry
from shapely.ops import transform
from shapely.geometry import Point
from shapely.geometry import Polygon
import geoviews as gv
from pystac_client import Client  
import xarray
import rasterio as rio
import rioxarray
import os
gv.extension('bokeh', 'matplotlib')

lon = -84.0106
lat = 10.432
buffer = 0.001    # decimal degrees equal to ~222 meters

loc = Point(lon, lat).buffer(buffer)

type(loc)

shapely.geometry.polygon.Polygon

base = gv.tile_sources.EsriImagery.opts(width=650, height=500)
sample_loc = gv.Polygons(loc).opts(line_color='yellow', line_width=10, color=None)
base * sample_loc

Translate the shapely polygon to a GeoJSON Polygon Object, which is what the STAC API expects.

roi = json.loads(geopandas.GeoSeries([loc]).to_json())
roi

{'type': 'FeatureCollection',
 'features': [{'id': '0',
   'type': 'Feature',
   'properties': {},
   'geometry': {'type': 'Polygon',
    'coordinates': [[[-84.00959999999999, 10.432],
      [-84.00960481527332, 10.43190198285967],
      [-84.00961921471959, 10.431804909677984],
      [-84.00964305966427, 10.431709715322746],
      [-84.00967612046749, 10.431617316567635],
      [-84.00971807873565, 10.431528603263175],
      [-84.00976853038769, 10.43144442976698],
      [-84.00982698954664, 10.431365606715836],
      [-84.0098928932188, 10.431292893218814],
      [-84.00996560671584, 10.431226989546637],
      [-84.01004442976698, 10.431168530387698],
      [-84.01012860326317, 10.431118078735652],
      [-84.01021731656763, 10.43107612046749],
      [-84.01030971532275, 10.431043059664269],
      [-84.01040490967799, 10.431019214719598],
      [-84.01050198285967, 10.431004815273328],
      [-84.0106, 10.431000000000001],
      [-84.01069801714033, 10.431004815273328],
      [-84.010795090322, 10.431019214719598],
      [-84.01089028467725, 10.431043059664269],
      [-84.01098268343236, 10.43107612046749],
      [-84.01107139673682, 10.431118078735652],
      [-84.01115557023302, 10.431168530387698],
      [-84.01123439328416, 10.431226989546637],
      [-84.01130710678119, 10.431292893218814],
      [-84.01137301045335, 10.431365606715836],
      [-84.0114314696123, 10.43144442976698],
      [-84.01148192126435, 10.431528603263175],
      [-84.0115238795325, 10.431617316567635],
      [-84.01155694033572, 10.431709715322746],
      [-84.0115807852804, 10.431804909677984],
      [-84.01159518472667, 10.43190198285967],
      [-84.0116, 10.432],
      [-84.01159518472667, 10.43209801714033],
      [-84.0115807852804, 10.432195090322017],
      [-84.01155694033572, 10.432290284677254],
      [-84.0115238795325, 10.432382683432365],
      [-84.01148192126435, 10.432471396736826],
      [-84.0114314696123, 10.43255557023302],
      [-84.01137301045335, 10.432634393284165],
      [-84.01130710678119, 10.432707106781187],
      [-84.01123439328416, 10.432773010453364],
      [-84.01115557023302, 10.432831469612303],
      [-84.01107139673682, 10.432881921264348],
      [-84.01098268343236, 10.432923879532511],
      [-84.01089028467725, 10.432956940335732],
      [-84.010795090322, 10.432980785280403],
      [-84.01069801714033, 10.432995184726673],
      [-84.0106, 10.433],
      [-84.01050198285967, 10.432995184726673],
      [-84.01040490967799, 10.432980785280403],
      [-84.01030971532275, 10.432956940335732],
      [-84.01021731656763, 10.432923879532511],
      [-84.01012860326317, 10.432881921264348],
      [-84.01004442976698, 10.432831469612303],
      [-84.00996560671584, 10.432773010453364],
      [-84.0098928932188, 10.432707106781187],
      [-84.00982698954664, 10.432634393284165],
      [-84.00976853038769, 10.43255557023302],
      [-84.00971807873565, 10.432471396736826],
      [-84.00967612046749, 10.432382683432365],
      [-84.00964305966427, 10.432290284677254],
      [-84.00961921471959, 10.432195090322017],
      [-84.00960481527332, 10.43209801714033],
      [-84.00959999999999, 10.432]]]},
   'bbox': [-84.0116, 10.431000000000001, -84.00959999999999, 10.433]}],
 'bbox': [-84.0116, 10.431000000000001, -84.00959999999999, 10.433]}

GeoJSON Polygon Object vertices are expected to be in counterclockwise order. The current object vertices are in clockwise order. Below we reverse the order of the vertices and assign them back to our roi object.

roi['features'][0]['geometry']['coordinates'][0] = roi['features'][0]['geometry']['coordinates'][0][::-1]

Find data using the STAC API

Specify the STAC API endpoint

STAC_URL = 'https://cmr.earthdata.nasa.gov/stac'

Connect to the LPCLOUD catalog using PySTAC Client

catalog = Client.open(f'{STAC_URL}/LPCLOUD/')

Submit a query for: - HLSL30 and HLSS30 - The roi’s geometry - Dates between 05-01-2021 and 06-30-2022

search = catalog.search(
    collections = ['HLSL30.v2.0', 'HLSS30.v2.0'],
    intersects = roi['features'][0]['geometry'],
    datetime = '2021-05/2022-06'
)

search.matched()

Get all the STAC Items/granules

items = search.get_all_items()

len(items)

Get the first STAC Item from the list and print off the available assets

items[0].assets

{'B07': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B07.tif>,
 'B01': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B01.tif>,
 'B04': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B04.tif>,
 'B06': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B06.tif>,
 'VZA': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.VZA.tif>,
 'Fmask': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.Fmask.tif>,
 'B02': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B02.tif>,
 'SAA': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.SAA.tif>,
 'B11': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B11.tif>,
 'B03': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B03.tif>,
 'VAA': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.VAA.tif>,
 'B10': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B10.tif>,
 'B05': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B05.tif>,
 'B09': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B09.tif>,
 'SZA': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.SZA.tif>,
 'browse': <Asset href=https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-public/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.jpg>,
 'metadata': <Asset href=https://cmr.earthdata.nasa.gov/search/concepts/G2144745416-LPCLOUD.xml>}

Grab the URL for the B04 asset

data_url = items[0].assets['B04'].href
data_url

'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T16PHS.2021133T155359.v2.0/HLS.L30.T16PHS.2021133T155359.v2.0.B04.tif'

NOTE: The above returns the HTTPS URL. For direct S3 access of cloud assets the S3 URL is need. See the following resource on how to generate the S3 URL: https://nasa-openscapes.github.io/2021-Cloud-Hackathon/tutorials/02_Data_Discovery_CMR-STAC_API.html

Read in the URL using `rioxarray`

NOTE: The example below show how to access the asset via HTTPS. An example for performing direct S3 access can be found here: https://nasa-openscapes.github.io/2021-Cloud-Hackathon/tutorials/05_Data_Access_Direct_S3.html

First we need to set our gdal configs. Make sure a .netrc file has been created in your home directory

rio_env = rio.Env(GDAL_DISABLE_READDIR_ON_OPEN='EMPTY',
                  GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),
                  GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt'))
rio_env.__enter__()

<rasterio.env.Env at 0x7f846be57f70>

da = rioxarray.open_rasterio(data_url)

da

Drop the ‘band’ coordinate

da.squeeze('band', drop=True)

<xarray.DataArray (y: 3660, x: 3660)>
[13395600 values with dtype=int16]
Coordinates:
  * x            (x) float64 8e+05 8e+05 8.001e+05 ... 9.097e+05 9.098e+05
  * y            (y) float64 1.2e+06 1.2e+06 1.2e+06 ... 1.09e+06 1.09e+06
    spatial_ref  int64 0
Attributes:
    _FillValue:    -9999.0
    scale_factor:  0.0001
    add_offset:    0.0
    long_name:     Red

xarray.DataArray

y: 3660
x: 3660

...
```
[13395600 values with dtype=int16]
```
Coordinates: (3)
- x
  (x)
  float64
  8e+05 8e+05 ... 9.097e+05 9.098e+05
```
array([799995., 800025., 800055., ..., 909705., 909735., 909765.])
```
- y
  (y)
  float64
  1.2e+06 1.2e+06 ... 1.09e+06
```
array([1199985., 1199955., 1199925., ..., 1090275., 1090245., 1090215.])
```
- spatial_ref
  ()
  int64
  0
  crs_wkt :
  PROJCS["UTM Zone 16, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
  semi_major_axis :
  6378137.0
  semi_minor_axis :
  6356752.314245179
  inverse_flattening :
  298.257223563
  reference_ellipsoid_name :
  WGS 84
  longitude_of_prime_meridian :
  0.0
  prime_meridian_name :
  Greenwich
  geographic_crs_name :
  Unknown datum based upon the WGS 84 ellipsoid
  horizontal_datum_name :
  Not_specified_based_on_WGS_84_spheroid
  projected_crs_name :
  UTM Zone 16, Northern Hemisphere
  grid_mapping_name :
  transverse_mercator
  latitude_of_projection_origin :
  0.0
  longitude_of_central_meridian :
  -87.0
  false_easting :
  500000.0
  false_northing :
  0.0
  scale_factor_at_central_meridian :
  0.9996
  spatial_ref :
  PROJCS["UTM Zone 16, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
  GeoTransform :
  799980.0 30.0 0.0 1200000.0 0.0 -30.0
```
array(0)
```
Attributes: (4)
_FillValue :
-9999.0
scale_factor :
0.0001
add_offset :
0.0
long_name :
Red

Find data using the STAC API

Read in the URL using rioxarray

Read in the URL using `rioxarray`