Accessing data on the USGS LPDAAC server via OPeNDAP protocol

Author

Eli Holmes (NOAA)

Overview

The USGS LPDAAC server requires NASA Earthdata login authentication, but doesn’t require a EULA (as far as I know). However they have similar redirect issues as data that does require a EULA. I could not get the data into an xarray.

Prerequisites

I assume you have a .netrc file at ~ (home). ~/.netrc should look just like this with your username and password.

machine urs.earthdata.nasa.gov
        login yourusername
        password yourpassword

Packages

import xarray as xr
import pydap.client
import pydap

Set up a token-based session with `earthaccess`

# create an authenticated session
import earthaccess
earthaccess.login()
edl_session = earthaccess.get_requests_https_session()

# if it keeps asking for your password, run this instead to save to .netrc file
auth = earthaccess.login()
# are we authenticated?
if not auth.authenticated:
    # ask for credentials and persist them in a .netrc file
    auth.login(strategy="interactive", persist=True)

USGS LPDAAC

Their OPeNDAP server also uses NASA Earthdata login authentication.

url="https://opendap.cr.usgs.gov/opendap/hyrax/MOD13Q1.061/h09v06.ncml"

# this works
pydap_ds = pydap.client.open_url(url, session=edl_session, protocol="dap4")

# this works too
store = xr.backends.PydapDataStore(pydap_ds)
ds = xr.open_dataset(store)

# this works too
url="dap4://opendap.cr.usgs.gov/opendap/hyrax/MOD13Q1.061/h09v06.ncml"
ds = xr.open_dataset(url, engine="pydap", session=edl_session)

print(f"Dataset size: {ds.nbytes/1e9:.2f} Tb")

Dataset size: 690.46 Tb

ds

# loading data works
ds["_250m_16_days_VI_Quality"].isel({"/time": 1, "/YDim": 1}).load()

<xarray.DataArray '_250m_16_days_VI_Quality' (/XDim: 4800)> Size: 19kB
array([2116., 2116., 2116., ..., 2116., 2116., 2116.], dtype=float32)
Coordinates:
    Latitude   (/XDim) float64 38kB 30.0 30.0 30.0 30.0 ... 30.0 30.0 30.0 30.0
    Longitude  (/XDim) float64 38kB -103.9 -103.9 -103.9 ... -92.38 -92.37
Dimensions without coordinates: /XDim
Attributes:
    long_name:     250m 16 days VI Quality
    units:         bit field
    valid_range:   [0, 65534]
    Legend:        \n\t Bit Fields Description (Right to Left): \n\t[0-1] : M...
    grid_mapping:  MODIS_Grid_16DAY_250m_500m_VI_eos_cf_projection
    Maps:          ()

xarray.DataArray

'_250m_16_days_VI_Quality'

/XDim: 4800

2.116e+03 2.116e+03 2.116e+03 ... 2.116e+03 2.116e+03 2.116e+03

array([2116., 2116., 2116., ..., 2116., 2116., 2116.], dtype=float32)

Coordinates: (2)
- Latitude
  (/XDim)
  float64
  30.0 30.0 30.0 ... 30.0 30.0 30.0
  long_name :
  Latitude
  units :
  degrees_north
  Maps :
  ()
```
array([29.996875, 29.996875, 29.996875, ..., 29.996875, 29.996875,
       29.996875])
```
- Longitude
  (/XDim)
  float64
  -103.9 -103.9 ... -92.38 -92.37
  long_name :
  Longitude
  units :
  degrees_east
  Maps :
  ()
```
array([-103.91857343, -103.91616788, -103.91376233, ...,  -92.37914828,
        -92.37674273,  -92.37433718])
```
Indexes: (0)
Attributes: (6)
long_name :
250m 16 days VI Quality
units :
bit field
valid_range :
[0, 65534]
Legend :
Bit Fields Description (Right to Left): [0-1] : MODLAND_QA [2 bit range] 00: VI produced, good quality 01: VI produced, but check other QA 10: Pixel produced, but most probably cloudy 11: Pixel not produced due to other reasons than clouds [2-5] : VI usefulness [4 bit range] 0000: Highest quality 0001: Lower quality 0010..1010: Decreasing quality 1100: Lowest quality 1101: Quality so low that it is not useful 1110: L1B data faulty 1111: Not useful for any other reason/not processed [6-7] : Aerosol quantity [2 bit range] 00: Climatology 01: Low 10: Average 11: High (11) [8] : Adjacent cloud detected; [1 bit range] 1: Yes 0: No [9] : Atmosphere BRDF correction performed [1 bit range] 1: Yes 0: No [10] : Mixed clouds [1 bit range] 1: Yes 0: No [11-13] : Land/Water Flag [3 bit range] 000: Shallow ocean 001: Land (Nothing else but land) 010: Ocean coastlines and lake shorelines 011: Shallow inland water 100: Ephemeral water 101: Deep inland water 110: Moderate or continental ocean 111: Deep ocean [14] : Possible snow/ice [1 bit range] 1: Yes 0: No [15] : Possible shadow [1 bit range] 1: Yes 0: No
grid_mapping :
MODIS_Grid_16DAY_250m_500m_VI_eos_cf_projection
Maps :
()

# this is another url that works with pydap
url = "https://opendap.cr.usgs.gov/opendap/hyrax/ECOSTRESS/ECO2CLD.001/2018.07.09/ECOSTRESS_L2_CLOUD_00048_001_20180709T204901_0601_03.h5"
pydap_ds = pydap.client.open_url(url, session=edl_session, protocol="dap4")

Setting token-based authentication with pydap

If you don’t want to use earthaccess, here is how to set up token-based authentication with pydap. Go on your Earthdata login user profile to create a token.

edl_token = "put your token here"
import requests
auth_hdr="Bearer " + edl_token
token_session = requests.Session()
token_session.headers={"Authorization": auth_hdr}

Conclusion

Working with USGS data on its OPeNDAP server requires Earthdata authentication. Use token-based authentication to solve redirect errors. earthaccess helps with this.