Accessing tiled NetCDFs at xarray (MSLSP30NA)

[mitchbon]

Hello,

I am interested in working with phenology data over my study area and want to use the MSLSP30NA product (https://lpdaac.usgs.gov/products/mslsp30nav011/) for that.

My goal is to have a gridded product over my area that I can use for easy visualization and as input into modelling. I am familiar with STAC and xarray and have had no issues building similar products using HLS (e.g., https://lpdaac.usgs.gov/products/hlsl30v002/), for example. I am trying to convert the NetCDF tiles that make up this dataset into this xarray format over a larger area. This would involve merging all the tiles along lat/lon and maybe time (e.g., get the mean value for variables between 2018/2019 years for each pixel).

I have tried a couple different ways to accomplish this but without success.
1. Earthaccess > xarray.open_mfdataset. If I narrow the earthaccess search to a single granule this works. However, if increase to more than one (e.g., a second time dimension) it will error out (e.g., ValueError: Could not find any dimension coordinates to use to order the datasets for concatenation). It is possible some combination of parameters will allow me to specify dimension coordinates but I am unsure.
2. Accessing the STAC directly using pystac_client and stackstac. The main issue here appears to be the lack of an EPSG metadata attribute 'out_epsg' (the x/y variables looks like lat/lon EPSG:4326).

I am not experienced with using NetCDF files, especially in an xarray environment. So I am not sure if the issue with how to load them or if there is an issue with this dataset (e.g., lack of set EPSG) that prevents it being manipulated in this way.

Thanks for your help.
Mitchell B.

[LP DAAC - dgolon]

Hi @mitchbon Let me talk to our STAC developers and one of us will get back to you. Please note: most of them are at AGU this week, so response might be delayed. Thanks --Danielle

[LP DAAC - dgolon]

Hi @mitchbon Would you mind emailing (lpdaac@usgs.gov) , or privately posting here, us a sample of your code so our developers can take a look at it?

[mitchbon]

I will send an email, thanks!

[LP DAAC - afriesz]

@mitchbon,

I'm not sure if I have an exact answer for you, but here are a few thoughts to consider.

- The MSLSP30NA data product is currently archived and distributed from the LP DAAC on-prem archive. In order to read these data remotely, you’ll need to pass your Earthdata Login credentials to authenticate and access these data. I have done this using Earthaccess and its use of fsspec. However, it does seem like remotely accessing these data is slow. You may be better off downloading the data first before working with them in Python/Xarray.

- pystac_client can certainly be used to find assets, but I'm uncertain if stackstac can read netCDF4 files as STAC assets. I’m also uncertain whether stackstac can pass authentication credentials via anything other than GDAL env variables.

- Finally, the MSLSP30NA data product is in a Universal Transverse Mercator (UTM) projection on the Military Grid Reference System (MGRS). When querying for granules using CMR and a region of interest, you may end up with granules that are from two different UTM zones. It’s recommended to only mosaic granules that are within a common zone. Mosaicking across zones will likely have negative impacts on your analysis. So, when trying to read granules into a Python object (i.e., xarray dataset) be sure that those granules are from within the same UTM zone. If you’re intent is to have a final mosaiced map, do the mosaicing after you perform your analysis by zone.

[mitchbon]

No worries for the delay, thanks for the tips! I will keep them in mind.