Open and plot COGs

Demonstrates opening and ploting COGs using pystac_client, xarray, and hvplot

Run this notebook

You can launch this notebook in VEDA JupyterHub by clicking the link below.

Launch in VEDA JupyterHub (requires access)

Learn more

Inside the Hub

This notebook was written on the VEDA JupyterHub and as such is designed to be run on a jupyterhub which is associated with an AWS IAM role which has been granted permissions to the VEDA data store via its bucket policy. The instance used provided 16GB of RAM.

See (VEDA Analytics JupyterHub Access)[https://nasa-impact.github.io/veda-docs/veda-jh-access.html] for information about how to gain access.

Outside the Hub

The data is in a protected bucket. Please request access by emailng aimee@developmentseed.org or alexandra@developmentseed.org and providing your affiliation, interest in or expected use of the dataset and an AWS IAM role or user Amazon Resource Name (ARN). The team will help you configure the cognito client.

You should then run:

%run -i 'cognito_login.py'

Approach

1. Use pystac_client to open the STAC catalog and retrieve the items in the collection
2. Open the collection with xarray and stackstac
3. Plot the data using hvplot

About the data

CDC’s Social Vulnerability Index (SVI) uses 15 variables at the census tract level. The data comes from the U.S. decennial census for the years 2000 & 2010, and the American Community Survey (ACS) for the years 2014, 2016, and 2018. It is a hierarchical additive index (Tate, 2013), with the component elements of CDC’s SVI including the following for 4 themes: Socioeconomic Status, Household Composition & Disability, Minority Status & Language, and Housing Type & Transportation.

SVI indicates the relative vulnerability of every U.S. Census tract–subdivisions of counties for which the Census collects statistical data. SVI ranks the tracts on 15 social factors, including unemployment, minority status, and disability, and further groups them into four related themes. Thus, each tract receives a ranking for each Census variable and for each of the four themes, as well as an overall ranking.

Scientific research

The SVI Overall Score provides the overall, summed social vulnerability score for a given tract. The Overall Score SVI Grid is part of the U.S. Census Grids collection, and displays the Center for Disease Control & Prevention (CDC) SVI score. Funding for the final development, processing and dissemination of this data set by the Socioeconomic Data and Applications Center (SEDAC) was provided under the U.S. National Aeronautics and Space Administration (NASA)¹.

The Overall SVI Score describes the vulnerability in a given county tract based on the combined percentile ranking of the four SVI scores (Socioeconomic Status, Household Composition & Disability, Minority Status & Language, and Housing Type & Transportation). The summed percentile ranking from the four themes is ordered, and then used to calculate an overall percentile ranking, ranging from 0 (less vulnerable) to 1 (more vulnerable)². Tracts with higher Overall SVI Scores typically rank high in other SVI domains, and reveal communities that may require extra support, resources, and preventative care in order to better prepare for and manage emergency situations.

Interpreting the data

The Overall SVI Score describes the vulnerability in a given county tract based on the combined percentile ranking of the four SVI scores (Socioeconomic Status, Household Composition & Disability, Minority Status & Language, and Housing Type & Transportation). The summed percentile ranking from the four themes is ordered, and then used to calculate an overall percentile ranking, ranging from 0 (less vulnerable) to 1 (more vulnerable)². Tracts with higher Overall SVI Scores typically rank high in other SVI domains, and reveal communities that may require extra support, resources, and preventative care in order to better prepare for and manage emergency situations.

Credits

• Center for International Earth Science Information Network, (CIESIN), Columbia University. 2021. Documentation for the U.S. Social Vulnerability Index Grids. Palisades, NY: NASA Socioeconomic Data and Applications Center (SEDAC). https://doi.org/10.7927/fjr9-a973. Accessed 13 May 2022.
• Centers for Disease Control and Prevention/ Agency for Toxic Substances and Disease Registry/ Geospatial Research, Analysis, and Services Program. CDC/ATSDR Social Vulnerability Index Database. https://www.atsdr.cdc.gov/placeandhealth/svi/documentation/pdf/SVI2018Documentation_01192022_1.pdf

from pystac_client import Client
import pandas as pd
import stackstac

import hvplot.xarray  # noqa

Declare your collection of interest

You can discover available collections the following ways:

• Programmatically: see example in the list-collections.ipynb notebook
• JSON API: https://staging-stac.delta-backend.com/collections
• STAC Browser: http://veda-staging-stac-browser.s3-website-us-west-2.amazonaws.com

STAC_API_URL = "https://openveda.cloud/api/stac/"
collection = "social-vulnerability-index-overall-nopop"

Find items in collection

Use pystac_client to search the STAC collection.

catalog = Client.open(STAC_API_URL)
search = catalog.search(collections=[collection])

item_collection = search.item_collection()
print(f"Found {len(item_collection)} items")

Found 5 items

Read data

Read in data using xarray using a combination of xpystac, stackstac, and rasterio.

da = stackstac.stack(item_collection, epsg=4326)
da = da.assign_coords({"time": da.start_datetime}).squeeze()
# da
da

/srv/conda/envs/notebook/lib/python3.11/site-packages/stackstac/prepare.py:408: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.
  times = pd.to_datetime(

<xarray.DataArray 'stackstac-6c45bf56d29299d91eb93a5bc48c5746' (time: 5,
                                                                y: 6298,
                                                                x: 13354)> Size: 3GB
dask.array<getitem, shape=(5, 6298, 13354), dtype=float64, chunksize=(1, 1024, 1024), chunktype=numpy.ndarray>
Coordinates: (12/17)
    id              (time) <U38 760B 'svi_2018_tract_overall_wgs84_nopop_cog'...
    band            <U11 44B 'cog_default'
  * x               (x) float64 107kB -178.2 -178.2 -178.2 ... -66.97 -66.97
  * y               (y) float64 50kB 71.38 71.37 71.37 ... 18.92 18.92 18.91
    end_datetime    (time) <U25 500B '2018-12-31T00:00:00+00:00' ... '2000-12...
    start_datetime  (time) <U25 500B '2018-01-01T00:00:00+00:00' ... '2000-01...
    ...              ...
    proj:geometry   object 8B {'type': 'Polygon', 'coordinates': [[[-178.2333...
    proj:shape      object 8B {6297, 13353}
    title           <U17 68B 'Default COG Layer'
    proj:bbox       object 8B {-178.23333334, 18.908332897999998, -66.9583337...
    epsg            int64 8B 4326
  * time            (time) <U25 500B '2018-01-01T00:00:00+00:00' ... '2000-01...
Attributes:
    spec:           RasterSpec(epsg=4326, bounds=(-178.24166595386018, 18.899...
    crs:            epsg:4326
    transform:      | 0.01, 0.00,-178.24|\n| 0.00,-0.01, 71.38|\n| 0.00, 0.00...
    resolution_xy:  (0.00833333330000749, 0.00833333329998412)

xarray.DataArray

'stackstac-6c45bf56d29299d91eb93a5bc48c5746'

• time: 5
• y: 6298
• x: 13354

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

  Array Chunk Bytes 3.13 GiB 8.00 MiB Shape (5, 6298, 13354) (1, 1024, 1024) Dask graph 490 chunks in 4 graph layers Data type float64 numpy.ndarray

• Coordinates: (17)
  • id
    (time)
    <U38
    'svi_2018_tract_overall_wgs84_no...

    array(['svi_2018_tract_overall_wgs84_nopop_cog',
           'svi_2016_tract_overall_wgs84_nopop_cog',
           'svi_2014_tract_overall_wgs84_nopop_cog',
           'svi_2010_tract_overall_wgs84_nopop_cog',
           'svi_2000_tract_overall_wgs84_nopop_cog'], dtype='<U38')

  • band
    ()
    <U11
    'cog_default'

    array('cog_default', dtype='<U11')

  • x
    (x)
    float64
    -178.2 -178.2 ... -66.97 -66.97

    array([-178.241666, -178.233333, -178.224999, ...,  -66.983333,  -66.975   ,
            -66.966666])

  • y
    (y)
    float64
    71.38 71.37 71.37 ... 18.92 18.91

    array([71.383333, 71.375   , 71.366666, ..., 18.925   , 18.916667, 18.908333])

  • end_datetime
    (time)
    <U25
    '2018-12-31T00:00:00+00:00' ... ...

    array(['2018-12-31T00:00:00+00:00', '2016-12-31T00:00:00+00:00',
           '2014-12-31T00:00:00+00:00', '2010-12-31T00:00:00+00:00',
           '2000-12-31T00:00:00+00:00'], dtype='<U25')

  • start_datetime
    (time)
    <U25
    '2018-01-01T00:00:00+00:00' ... ...

    array(['2018-01-01T00:00:00+00:00', '2016-01-01T00:00:00+00:00',
           '2014-01-01T00:00:00+00:00', '2010-01-01T00:00:00+00:00',
           '2000-01-01T00:00:00+00:00'], dtype='<U25')

  • proj:projjson
    ()
    object
    {'id': {'code': 4326, 'authority...

    array({'id': {'code': 4326, 'authority': 'EPSG'}, 'name': 'WGS 84', 'type': 'GeographicCRS', 'datum': {'name': 'World Geodetic System 1984', 'type': 'GeodeticReferenceFrame', 'ellipsoid': {'name': 'WGS 84', 'semi_major_axis': 6378137, 'inverse_flattening': 298.257223563}}, '$schema': 'https://proj.org/schemas/v0.4/projjson.schema.json', 'coordinate_system': {'axis': [{'name': 'Geodetic latitude', 'unit': 'degree', 'direction': 'north', 'abbreviation': 'Lat'}, {'name': 'Geodetic longitude', 'unit': 'degree', 'direction': 'east', 'abbreviation': 'Lon'}], 'subtype': 'ellipsoidal'}},
          dtype=object)

  • proj:epsg
    ()
    int64
    4326

    array(4326)

  • proj:wkt2
    ()
    <U277
    'GEOGCS["WGS 84",DATUM["WGS_1984...

    array('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]',
          dtype='<U277')

  • proj:transform
    ()
    object
    {0.00833333330000749, 0.0, 1.0, ...

    array({0.00833333330000749, 0.0, 1.0, -0.00833333329998412, 71.383332688, -178.23333334},
          dtype=object)

  • description
    ()
    <U47
    'Cloud optimized default layer t...

    array('Cloud optimized default layer to display on map', dtype='<U47')

  • proj:geometry
    ()
    object
    {'type': 'Polygon', 'coordinates...

    array({'type': 'Polygon', 'coordinates': [[[-178.23333334, 18.908332897999998], [-66.958333785, 18.908332897999998], [-66.958333785, 71.383332688], [-178.23333334, 71.383332688], [-178.23333334, 18.908332897999998]]]},
          dtype=object)

  • proj:shape
    ()
    object
    {6297, 13353}

    array({6297, 13353}, dtype=object)

  • title
    ()
    <U17
    'Default COG Layer'

    array('Default COG Layer', dtype='<U17')

  • proj:bbox
    ()
    object
    {-178.23333334, 18.9083328979999...

    array({-178.23333334, 18.908332897999998, -66.958333785, 71.383332688},
          dtype=object)

  • epsg
    ()
    int64
    4326

    array(4326)

  • time
    (time)
    <U25
    '2018-01-01T00:00:00+00:00' ... ...

    array(['2018-01-01T00:00:00+00:00', '2016-01-01T00:00:00+00:00',
           '2014-01-01T00:00:00+00:00', '2010-01-01T00:00:00+00:00',
           '2000-01-01T00:00:00+00:00'], dtype='<U25')

• Indexes: (3)
  • x
    PandasIndex

    PandasIndex(Index([-178.24166595386018, -178.23333262056016, -178.22499928726018,
           -178.21666595396016, -178.20833262066014, -178.19999928736013,
           -178.19166595406014, -178.18333262076013,  -178.1749992874601,
           -178.16666595416012,
           ...
            -67.04166639856027,  -67.03333306526025,  -67.02499973196025,
            -67.01666639866023,  -67.00833306536023,  -66.99999973206023,
            -66.99166639876022,  -66.98333306546022,   -66.9749997321602,
             -66.9666663988602],
          dtype='float64', name='x', length=13354))

  • y
    PandasIndex

    PandasIndex(Index([ 71.38333304766397,  71.37499971436398,    71.366666381064,
            71.35833304776402,  71.34999971446403,  71.34166638116405,
            71.33333304786406,  71.32499971456407,   71.3166663812641,
            71.30833304796411,
           ...
           18.983333257363824, 18.974999924063845, 18.966666590763857,
            18.95833325746387,  18.94999992416389, 18.941666590863903,
           18.933333257563923, 18.924999924263936, 18.916666590963956,
            18.90833325766397],
          dtype='float64', name='y', length=6298))

  • time
    PandasIndex

    PandasIndex(Index(['2018-01-01T00:00:00+00:00', '2016-01-01T00:00:00+00:00',
           '2014-01-01T00:00:00+00:00', '2010-01-01T00:00:00+00:00',
           '2000-01-01T00:00:00+00:00'],
          dtype='object', name='time'))

• Attributes: (4)

  spec :

  RasterSpec(epsg=4326, bounds=(-178.24166595386018, 18.899999924363982, -66.95833306556018, 71.38333304766397), resolutions_xy=(0.00833333330000749, 0.00833333329998412))

  crs :

  epsg:4326

  transform :

  | 0.01, 0.00,-178.24| | 0.00,-0.01, 71.38| | 0.00, 0.00, 1.00|

  resolution_xy :

  (0.00833333330000749, 0.00833333329998412)

There are 5 items representing the 5 years of data in the collection (2000, 2010, 2014, 2016, and 2018).

Plot data

Plot data using hvplot. By using rasterize=True we tell hvplot to use datashader behind the scenes to make the plot render more quickly and re-render on zoom.

%%time
da.hvplot(x="x", y="y", rasterize=True, clim=(0, 1), coastline=True, cmap="viridis", widget_location="bottom")

CPU times: user 1.33 s, sys: 87.9 ms, total: 1.42 s
Wall time: 1.61 s