[Python] MODIS LST tiles concatenation or merging issue with xarray

I am having an issue while concatenating two xarray data files. The two files (MODIS LST, accessed through the STAC API) are from the same date/time but correspond to adjacent geographical tiles (side-by-side, by longitude or x). When I use the xr.concat() function, everything works fine without any errors, but when I visualize the data, I notice many vertical white lines, figure below (which appear to be no-data values). When I plot the individual files (from the same time step), the data exists in the regions where the vertical lines appear (these figures not shown here).

A reproducible example is provided below if anyone would like to give it a try:

import pystac_client
import stackstac
import planetary_computer
import xarray as xr
import pandas as pd

catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1",
modifier=planetary_computer.sign_inplace,
)

# time range
time_range = "2000-10-01/2000-10-30"

# tile 1
bbox = [75, 30, 80, 32]
search = catalog.search(collections=["modis-11A2-061"], bbox=bbox, datetime=time_range)
items = search.get_all_items()
ds_d_t1 = stackstac.stack(items, epsg=4326, chunksize='auto', assets=["LST_Day_1km"], properties=False)
ds_d_t1 = ds_d_t1.rename("LST")

# tile 1
bbox = [66, 35, 67, 37]
search = catalog.search(collections=["modis-11A2-061"], bbox=bbox, datetime=time_range)
items = search.get_all_items()
ds_d_t2 = stackstac.stack(items, epsg=4326, chunksize='auto', assets=["LST_Day_1km"], properties=False)
ds_d_t2 = ds_d_t2.rename("LST")

# lets handle the time, because there is no proper time coordinate
# extract year and day of the year from MODIS ID
def get_date(modis_id):
year = int(modis_id.split('.')[1][1:5])
day_of_year = int(modis_id.split('.')[1][5:])
date = pd.to_datetime(str(year) + str(day_of_year), format='%Y%j')
return date

# create a new time coordinate variable using the get_date function
ds_d_t1['time'] = xr.DataArray([get_date(modis_id) for modis_id in ds_d_t1['id'].values],
dims='time', coords={'time': ds_d_t1['time']})

ds_d_t2['time'] = xr.DataArray([get_date(modis_id) for modis_id in ds_d_t2['id'].values],
dims='time', coords={'time': ds_d_t2['time']})

ds_d_t1 = ds_d_t1.sortby('time')
ds_d_t2 = ds_d_t2.sortby('time')

# concat
ds_m = xr.concat([ds_d_t1, ds_d_t2], dim="x")
ds_m = ds_m.sortby('x')

# visualise
ds_m[0].plot()


Here is the output figure, with many vertical white lines:

Any help in concatenating these files properly would be greatly appreciated.

Continue reading...