google / Xee

An Xarray extension for Google Earth Engine
Apache License 2.0
240 stars 28 forks source link

Transposing MODIS data introduces data artifacts #129

Open nielja opened 7 months ago

nielja commented 7 months ago

Dear developers,

thanks for this great package! I've been trying to work with MODIS MOD13Q1 data and ran into some inconsistencies when trying to mask values based on the quality band when transposing the data. This is my first time posting here, so I hope I'm following all guidelines!

It seems when loading the original MODIS data, the longitude and latitude variable are switched, so I transposed them. However, this leads to unexpected results when computing masks.

#Set dates
start_date = ee.Date.fromYMD(2001, 1, 1)
end_date = ee.Date.fromYMD(2024, 1, 1)

#Load MODIS NDVI data
ds = ee.ImageCollection("MODIS/061/MOD13Q1").filterDate(start_date, end_date)
dx = xr.open_dataset(ds, engine = "ee", crs = "EPSG:4326")
dt = dx.transpose('time','lat','lon')

Transposing the values introduces more of the shading artifacts already mentioned in #119.

#Plot transposed and original Summary QA bands

fig, ax = plt.subplots(1, 2, figsize = (8, 4))
dx.SummaryQA[0].plot(ax = ax[0])
ax[0].set_title("not transposed, full filter, then first timestep", fontsize = 8)
ax[0].set_xlim([-90, 90])
dt.SummaryQA[0].plot(ax = ax[1])
ax[1].set_title("not transposed, filter on first time step", fontsize = 8)
plt.tight_layout()

image

This leads to further issues when conducting masking. In the first case, I filtered the data for all time steps and then took the first time slice, and in the second case I computed the filter only on the first time slice. Theoretically, these should be the same, and they are in case of the non-transposed data (first row of the plots below), but not for the transposed data (second row).

#Do masking on full global data

fig, ax = plt.subplots(2, 2)
dx.SummaryQA.isin([0,1])[0].plot(ax = ax[0, 0])
ax[0,0].set_title("not transposed, full filter, then first timestep", fontsize = 8)
dx.SummaryQA[0].isin([0, 1]).plot(ax = ax[0, 1])
ax[0,1].set_title("not transposed, filter on first time step", fontsize = 8)
dt.SummaryQA.isin([0,1])[0].plot(ax = ax[1, 0])
ax[1,0].set_title("transposed, full filter, then first time step", fontsize = 8)
dt.SummaryQA[0].isin([0, 1]).plot(ax = ax[1, 1])
ax[1,1].set_title("transposed, filter on first time step", fontsize = 8)
plt.tight_layout()

image

This issue somehow gets worse when looking only at a subarea.

#Masking on subarea

dxa = dx.sel(lon = slice(-100, -30), lat = slice(9, -40))
dta = dt.sel(lon = slice(-100, -30), lat = slice(9, -40))

#Plot masks on subarea

fig, ax = plt.subplots(2, 2)
dxa.SummaryQA.isin([0,1])[0].plot(ax = ax[0, 0])
ax[0,0].set_title("not transposed, full filter, then first timestep", fontsize = 8)
dxa.SummaryQA[0].isin([0, 1]).plot(ax = ax[0, 1])
ax[0,1].set_title("not transposed, filter on first time step", fontsize = 8)
dta.SummaryQA.isin([0,1])[0].plot(ax = ax[1, 0])
ax[1,0].set_title("transposed, full filter, then first time step", fontsize = 8)
dta.SummaryQA[0].isin([0, 1]).plot(ax = ax[1, 1])
ax[1,1].set_title("transposed, filter on first time step", fontsize = 8)
plt.tight_layout()

image

Bottom left is the mask when transposing, then filtering on all values and taking the first time slice. It should look like the bottom right plot and seems to have some of the correct features (e.g. that little bump), but most values are incorrectly assigned.

My guess was that this has something to do with both the shading issue mentioned in #119 and some other artefacts introduced by the transposing.

naschmitz commented 7 months ago

Sounds like this is related to #119. We're looking into it.