Hello,
New user here, perhaps with a beginner question 
I am trying to collect Daymet data and MODIS LST imagery in one xarray.Dataset from Jan 1 to July 1, 2020 for an arbitrary 0.5 x 0.5 degree area. Both of these sources have a native resolution of 1km, and I’d like the result to be at that resolution as well.
I’m working through microsoft planetary computer, my approach is coming from their tutorials. Daymet comes as single Dataset indexed by time and the x/y coordinate of the Daymet projection (Lambert conformal conic). Since Daymet is quite large, I think it will be easier to use Dataset.sel on Daymet with its projected coordinates, then reproject MODIS to align them. So, I’m creating a proj4 string from the parameters provided in Daymet, and then passing that string to odc.stac.load when loading MODIS. The projection information in Daymet is:
false_easting :
0.0
false_northing :
0.0
grid_mapping_name :
lambert_conformal_conic
inverse_flattening :
298.257223563
latitude_of_projection_origin :
42.5
longitude_of_central_meridian :
-100.0
semi_major_axis :
6378137.0
standard_parallel :
[25.0, 60.0]
And the resulting proj4 string is:
+proj=lcc +lat_0=42.5 +lon_0=-100 +lat_1=25 +lat_2=60 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
(note that inverse_flattening and semi_major_axis are the defaults for a WGS84 ellipse)
I’m converting my lat/lon bounding box into Daymet coordinates via pyproj.Transformer and rounding to the nearest 1000m as below.
import odc.stac
import planetary_computer as pc
import pystac_client
import pystac
import pyproj
import fsspec
# Setup view location/size
latitude, longitude = 47.387, -119.232
buffer = 0.25
bbox = [longitude - buffer, latitude - buffer, longitude + buffer, latitude + buffer]
start = "2020-01-01"
end = "2020-07-01"
daymet_crs = pyproj.CRS.from_proj4("""
+proj=lcc +lat_0=42.5 +lon_0=-100 +lat_1=25 +lat_2=60 +x_0=0 +y_0=0
+ellps=WGS84 +units=m +no_defs +type=crs
""".replace("\n", " "))
latlon_daymet_transform = pyproj.Transformer.from_crs("EPSG:4326", daymet_crs, always_xy=True)
# Snap to the nearest km
xmin, ymin = tuple(map(lambda x: round(x, -3), latlon_daymet_transform.transform(bbox[0], bbox[1])))
xmax, ymax = tuple(map(lambda x: round(x, -3), latlon_daymet_transform.transform(bbox[2], bbox[3])))
Then I’m slicing both datasets using those bounds. First, Daymet:
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=pc.sign_inplace,
)
daymet_collection = pystac.read_file(
"https://planetarycomputer.microsoft.com/api/stac/v1/collections/daymet-daily-na"
)
asset = daymet_collection.assets['zarr-https']
store = fsspec.get_mapper(asset.href)
daymet = xr.open_zarr(store, **asset.extra_fields["xarray:open_kwargs"])
daymet_subset = ds.sel(
time=slice(start, end),
x = slice(xmin, xmax),
y = slice(ymax, ymin)
)
Then, MODIS:
modis_search = catalog.search(
collections=["modis-11A1-061"],
bbox=bbox,
datetime=[start,end],
# afternoon overpasses only
query={"platform": {"eq": "aqua"}}
)
modis_items = modis_search.get_all_items()
modis_data = odc.stac.load(
modis_items,
crs=daymet_crs,
bands=["LST_Day_1km", "QC_Day", "Day_view_time", "Day_view_angl"],
resolution=1000,
x=(xmin, xmax),
y=(ymin, ymax)
)
This all works, but the coordinates are misaligned by 250m.
I’ve tried adjusting the bounding box and modifying the Daymet proj4 string, but to no avail. I know that stackstac has a RasterSpec that can set the bounding box/resolution explicitly, but it only accepts projections with an EPSG code. Daymet’s projection does not have an EPSG code, so I cannot use that approach without reprojecting Daymet (very slow).
I am feeling stuck, and would appreciate any guidance on combining these two datasets.