Example Notebooks

medaprep.skim and medaprep.visualize modules to investigate Sentinel-2 data¶

In this notebook, we will show how using medaprep.skim.features,medaprep.visualize.query, and medaprep.visualize.distributions can help prepare and visualize large sets of data. We will load a set of Sentinel-2 imagery over Austin, TX and quickly retreive statistics about the features within to assist with downstream processing. We will use the following process:

1. Initialize Dask cluster to load xarray Dataset from a public S3 bucket.
2. visualize.query bounding boxes to check geographic coverage.
3. skim.features from the dataset for a quick set of summary statistics.
4. visualize.distributions of variables within the dataset.

In [1]:

from bokeh.io import show, output_notebook
from bokeh.layouts import layout
import dask.distributed
import folium
import geopandas as gpd
from IPython.display import HTML, display
from odc.stac import configure_rio, stac_load
import pandas as pd
from medaprep import skim, visualize
from pystac_client import Client
import shapely.geometry
import xarray

client = dask.distributed.Client()
configure_rio(cloud_defaults=True, aws={"aws_unsigned": True}, client=client)

/Users/dylan/opt/anaconda3/envs/odc-stac/lib/python3.10/site-packages/distributed/node.py:179: UserWarning: Port 8787 is already in use.
Perhaps you already have a cluster running?
Hosting the HTTP server on port 49607 instead
  warnings.warn(

Load Sentinel-2 data from public STAC catalog¶

The following cell loads an xarray Dataset containing Sentinel-2 imagery from a STAC catalog. We initialize the center point of our query to be the center of Austin, TX and search for imagery from July, 2022 within a 100 km radius.

In [2]:

config = {
    "sentinel-s2-l2a-cogs": {
        "assets": {
            "*": {"data_type": "uint16", "nodata": 0},
            "SCL": {"data_type": "uint8", "nodata": 0},
            "visual": {"data_type": "uint8", "nodata": 0},
        },
        "aliases": {"red": "B04", "green": "B03", "blue": "B02"},
    },
    "*": {"warnings": "ignore"},
}

km2deg = 1.0 / 111

x, y = (-97.744, 30.266)

r = 100 * km2deg

bbox = (x - r, y - r, x + r, y + r)

catalog = Client.open("https://earth-search.aws.element84.com/v0")

query = catalog.search(
    collections=["sentinel-s2-l2a-cogs"], datetime="2022-07-01/2022-07-31", limit=100, bbox=bbox
)

items = list(query.get_items())
stac_json = query.get_all_items_as_dict()

crs = "epsg:3857"
zoom = 2**5

data = stac_load(
    items,
    crs=crs,
    resolution=10*zoom,
    chunks={},
    groupby="solar_day",
    stac_cfg=config,
)
data

Out[2]:

<xarray.Dataset>
Dimensions:      (y: 1142, x: 1137, time: 12)
Coordinates:
  * y            (y) float64 3.716e+06 3.715e+06 ... 3.351e+06 3.351e+06
  * x            (x) float64 -1.102e+07 -1.102e+07 ... -1.066e+07 -1.066e+07
    spatial_ref  int32 3857
  * time         (time) datetime64[ns] 2022-07-03T17:25:22 ... 2022-07-30T17:...
Data variables: (12/16)
    visual       (time, y, x) uint8 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    B01          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    B02          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    B03          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    B04          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    B05          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    ...           ...
    B09          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    B11          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    B12          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    AOT          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    WVP          (time, y, x) uint16 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>
    SCL          (time, y, x) uint8 dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • y: 1142
  • x: 1137
  • time: 12
• Coordinates: (4)
  • y
    (y)
    float64
    3.716e+06 3.715e+06 ... 3.351e+06

    units :

    metre

    resolution :

    -320.0

    crs :

    EPSG:3857

    array([3715680., 3715360., 3715040., ..., 3351200., 3350880., 3350560.])

  • x
    (x)
    float64
    -1.102e+07 ... -1.066e+07

    units :

    metre

    resolution :

    320.0

    crs :

    EPSG:3857

    array([-11020640., -11020320., -11020000., ..., -10657760., -10657440.,
           -10657120.])

  • spatial_ref
    ()
    int32
    3857

    spatial_ref :

    PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["Popular Visualisation Pseudo-Mercator",METHOD["Popular Visualisation Pseudo Mercator",ID["EPSG",1024]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Web mapping and visualisation."],AREA["World between 85.06°S and 85.06°N."],BBOX[-85.06,-180,85.06,180]],ID["EPSG",3857]]

    crs_wkt :

    PROJCRS["WGS 84 / Pseudo-Mercator",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["Popular Visualisation Pseudo-Mercator",METHOD["Popular Visualisation Pseudo Mercator",ID["EPSG",1024]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Web mapping and visualisation."],AREA["World between 85.06°S and 85.06°N."],BBOX[-85.06,-180,85.06,180]],ID["EPSG",3857]]

    array(3857, dtype=int32)

  • time
    (time)
    datetime64[ns]
    2022-07-03T17:25:22 ... 2022-07-...

    array(['2022-07-03T17:25:22.000000000', '2022-07-05T17:15:15.000000000',
           '2022-07-08T17:25:13.000000000', '2022-07-10T17:15:23.000000000',
           '2022-07-13T17:25:21.000000000', '2022-07-15T17:15:14.000000000',
           '2022-07-18T17:25:13.000000000', '2022-07-20T17:15:21.000000000',
           '2022-07-23T17:25:19.000000000', '2022-07-25T17:15:15.000000000',
           '2022-07-28T17:25:14.000000000', '2022-07-30T17:15:20.000000000'],
          dtype='datetime64[ns]')

• Data variables: (16)
  • visual
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 14.86 MiB 1.24 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint8 numpy.ndarray

  • B01
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B02
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B03
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B04
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B05
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B06
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B07
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B08
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B8A
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B09
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B11
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • B12
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • AOT
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • WVP
    (time, y, x)
    uint16
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 29.72 MiB 2.48 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint16 numpy.ndarray

  • SCL
    (time, y, x)
    uint8
    dask.array<chunksize=(1, 1142, 1137), meta=np.ndarray>

    nodata :

    0

    Array Chunk Bytes 14.86 MiB 1.24 MiB Shape (12, 1142, 1137) (1, 1142, 1137) Count 98 Tasks 12 Chunks Type uint8 numpy.ndarray

• Attributes: (0)

visualize the bounding boxes of the input query and the returned data using visualize.query¶

visualize.query returns a visualization of the input and returned bounding boxes.

In [3]:

gdf = gpd.GeoDataFrame.from_features(stac_json, 3857)

bbox2 = (gdf.bounds['minx'].min(),gdf.bounds['miny'].min(),gdf.bounds['maxx'].max(),gdf.bounds['maxy'].max())
bbox3 = tuple([b+2 for b in bbox2])

visualize.query(bbox=[bbox,bbox2],
                    name=["Query","Returned"],
                    folium_map=folium.Map(),
                    color=["red","green"])

Out[3]:

Make this Notebook Trusted to load map: File -> Trust Notebook

skim through features of the data using skim.features¶

skim.features returns basic statistics and properties of an xarray Dataset.

In [4]:

small_data = data.sel(time='2022-07-03').compute()

df = skim.features(small_data)
df

Out[4]:

	variables	data_types	NaNs	mean	std	maximums	minimums
0	visual	uint8	False	83.206896	85.090139	255	0
1	B01	uint16	False	694.111545	1057.890177	12671	0
2	B02	uint16	False	715.116638	1019.690505	11615	0
3	B03	uint16	False	831.311709	1033.465046	10787	0
4	B04	uint16	False	897.211163	1064.686861	10284	0
5	B05	uint16	False	1144.481014	1258.273113	11399	0
6	B06	uint16	False	1453.429606	1447.058901	10691	0
7	B07	uint16	False	1585.919234	1543.215470	10361	0
8	B08	uint16	False	1592.048660	1549.807719	10519	0
9	B8A	uint16	False	1737.786114	1662.727114	10352	0
10	B09	uint16	False	2445.116097	3098.090223	15918	0
11	B11	uint16	False	1869.718284	1790.813024	8431	0
12	B12	uint16	False	1358.351765	1357.018681	8239	0
13	AOT	uint16	False	59.987095	57.874945	173	0
14	WVP	uint16	False	1847.768434	1681.527584	6201	0
15	SCL	uint8	False	3.550229	3.664158	10	0

visualize the distributions of the variables using visualize.distributions¶

visualize.distributions returns a list of bokeh figures containing estimated distributions of the variables.

In [5]:

output_notebook()
show(layout(visualize.distributions(small_data, df, 500)))

Loading BokehJS ...