Okay, so it seems that this might quickly become a performance nightmare. I am now looking into doing something like https://arxiv.org/pdf/1703.09112.pdf because this nicely handles imbalanced data in the different dimensions.
The linear model for coregionalization does involve computing a Kronecker product though and this seems to be quite inefficient as not hard-coded in STAN (A homemade kronecker product function).
Any ideas on how to cope with this efficiently?