Graph learning aims to infer a network structure directly from observed data, enabling the analysis of complex dependencies in irregular domains. Traditional methods focus on scalar signals at each node, ignoring dependencies along additional dimensions such as time, configurations of the observation device, or populations. In this work, we propose a graph signal processing framework for learning graphs from two-dimensional signals, modeled as matrix graph signals generated by joint filtering along both dimensions. This formulation leverages the concept of graph stationarity across the two dimensions and leverages product graph representations to capture structured dependencies. Based on this model, we design an optimization problem that can be solved efficiently and provably recovers the optimal underlying Kronecker/Cartesian/strong product graphs. Experiments on synthetic data demonstrate that our approach achieves higher estimation accuracy and reduced computational cost compared to existing methods.