使用xarray拼接MODIS LST相邻瓦片出现无效值问题求助

问题根源 #

你拼接的两个MODIS瓦片纬度范围完全不重叠：第一个瓦片纬度覆盖30-32，第二个覆盖35-37。直接按x维度拼接时，xarray会对y坐标做外对齐——保留所有y值，其中一个瓦片在另一个的y范围外的区域自动填充无数据值（NaN），最终可视化就会出现大量垂直白线。另外你代码里把第二个瓦片标注为“tile 1”是笔误，应该改为“tile 2”。

修正步骤 #

1. 选择正确的相邻瓦片 #

调整两个瓦片的bbox，确保它们在同一纬度范围内，仅沿经度方向相邻。比如：

• 第一个瓦片bbox：[70, 30, 75, 32]（经度70-75，纬度30-32）
• 第二个瓦片bbox：[75, 30, 80, 32]（经度75-80，纬度30-32）

这样两个瓦片的y坐标（纬度）完全匹配，拼接后不会出现无数据区域。

2. 使用正确的拼接方式 #

如果两个瓦片的y坐标一致，直接用xr.concat即可；如果y坐标有部分重叠但不完全一致，建议用xr.merge按空间坐标对齐，而非强制按维度拼接：

# 替换原concat代码
ds_m = xr.merge([ds_d_t1, ds_d_t2])

xr.merge会根据x和y的坐标值自动对齐数据，保留所有有效的空间点，避免因维度坐标不匹配产生的无数据。

3. 验证坐标一致性 #

拼接前可以先检查两个数据集的y坐标是否匹配：

print("Tile1 y范围:", ds_d_t1.y.min().values, ds_d_t1.y.max().values)
print("Tile2 y范围:", ds_d_t2.y.min().values, ds_d_t2.y.max().values)

如果输出的y范围一致，再执行拼接操作。

修正后的完整代码示例 #

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 = "2000-10-01/2000-10-30"

# Tile 1：经度70-75，纬度30-32
bbox = [70, 30, 75, 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 2：经度75-80，纬度30-32（与Tile1同纬度，经度相邻）
bbox = [75, 30, 80, 32]
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")

# 处理时间坐标（可选：stackstac已生成time坐标，可验证后选择是否保留）
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

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')

# 验证y坐标一致性
print("Tile1 y范围:", ds_d_t1.y.min().values, ds_d_t1.y.max().values)
print("Tile2 y范围:", ds_d_t2.y.min().values, ds_d_t2.y.max().values)

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

# 可视化
ds_m[0].plot()

额外提示 #

• MODIS的瓦片本身是按正弦投影的分幅（如h/v编号），如果要准确获取相邻的官方瓦片，建议直接按h和v编号搜索，而非手动指定bbox，避免范围不匹配。
• stackstac生成的数据集默认已经包含正确的time坐标，你可以通过print(ds_d_t1.time)验证，无需手动提取，能减少出错概率。

内容的提问来源于stack exchange，提问作者Arindan Mandal