jgomezdans / KaFKA

Kascading Fast Kalman Assimilation
GNU General Public License v3.0
5 stars 8 forks source link

Kernels from do_MCD43_prior.py #14

Open NPounder opened 7 years ago

NPounder commented 7 years ago

priorkernels

Here are the kernels generated. The labels match the labels in the code: Top row should be kernels and middle row should be uncertainties. Bottom row is the BHR (and unc) calculated from these kernels. This is all for modis band 1. The BHR is higher than the BHR product I'm typically seeing - the uncertainty values are closer to typical BHR values but they're clearly not swapped in the code.

NPounder commented 7 years ago

priorbhr BHR plots

jgomezdans commented 7 years ago

You sure you go the BHR transformation OK? Third weight is -ve (-1.37762)... (but might be labelling)

NPounder commented 7 years ago

Good spot! That has improved things, although it is still quite high. priorbhr I'll fix in the KaFKA code to.

My plotting code:

band = 1
file_prior = '/data/selene/ucfajlg/Aurade_MODIS/MCD4/MCD43_average_2009_300_367_b{}.tif'.format(band)
dst_prior = gdal.Open(file_prior)
P_diag = dst_prior.ReadAsArray()

fig, axs = plt.subplots(2,3, figsize=(21,14))
labels = ['P1', 'P2', 'P3', 'S1', 'S2', 'S3']
for i, ax in enumerate(axs.flatten()):
    im = ax.imshow(P_diag[i,...], vmin=0, vmax = 0.5)
    ax.set_title(labels[i])
    plt.colorbar(im, ax=ax,fraction=0.046, pad=0.04)

fig, ax2 = plt.subplots(1,2, figsize=(14,7))
ax2=ax2.flatten()
bhr = P_diag[0,...]+P_diag[1,...]*0.189184 - P_diag[2,...]*1.377622
im2 = ax2[0].imshow(bhr, vmin=0, vmax=0.5)
ax2[0].set_title('BHR = P1 + 0.189184*P2 - 1.377622*P3')
plt.colorbar(im2, ax=ax2[0], fraction=0.046, pad=0.04)
bhr = P_diag[3,...]+P_diag[4,...]*0.189184 - P_diag[5,...]*1.377622
im2 = ax2[1].imshow(bhr, vmin=0, vmax=0.5)
ax2[1].set_title('BHR unc = S1 + 0.189184*S2 - 1.377622*S3')
plt.colorbar(im2, ax=ax2[1], fraction=0.046, pad=0.04)
jgomezdans commented 7 years ago

I'm thinking that it might be better to show a comparison between "synergy BHR/iso/vol/geo" and "MCD43 BHR/iso/vol/geo" based on a z-score (e.g. ((iso_synergy - iso_mcd3) /iso_synergy_std)). Also, check timing, as this can be shifted between products (we talked about this, just putting it here for references).

NPounder commented 7 years ago

I can do that.

Just to be clear, the plots above aren't comparison, it's BHR and unc. It's just that when I plot any MCD43 bhr for band 1 the values (ignoring bright spots) are typically less than 0.15 (by which I mean the image is blue) but here the image is mostly turquoise. So either I'm not getting the right bands/products or I have another bug in the calcs or there is a bug in the prior estimate or MODIS because unless 2009 was very different I'd expect the plots to be of the same order here.

jgomezdans commented 7 years ago

Ah, OK I see... The current version in master uses data from southern Spain + Morocco, and starts on DoY 45 to avoid snow and clouds. I have also tightened the model uncertainty a lot (but that's just eyeballing things), meaning that things shouldn't deviate much from the calculated prior for a while, but haven't tested against MCD43 official

NPounder commented 7 years ago

OK. I'll set something running here too and see 1) if it runs and 2) how it looks against MCD43.

Just to check - have you already made the prior files or do I need to make them (for the new site)?

jgomezdans commented 7 years ago

I think it should all work (prior files do exist already). You will need to change the output directory (2nd argument of create_output_file) to somewhere you can write to. Also, creating the output is quite wasteful as it creates 366 bands for what's effectively a 45 day example!!!