Data-driven discovery of governing equations for fluid systems is challenged by sparse, noisy measurements and unobserved variables like pressure. This study addresses these challenges in a decaying, nearly incompressible fluid system, characterized by localized sharp spatial gradients, using a dataset of three-dimensional velocity and density fields on a periodic domain across only 10 time slices. To generate robust spatio-temporal features, temporal derivatives were estimated using second-order local polynomial regression for each spatial point, mitigating noise from sparse data. Spatial derivatives were computed by benchmarking spectral methods against a 5th-order Weighted Essentially Non-Oscillatory (WENO5) scheme, with the selection guided by adherence to the incompressibility constraint (). A comprehensive feature library was then constructed, including advection, viscous, and various proxy terms for pressure and buoyancy, motivated by a preliminary momentum residual analysis. All primary fields and generated features were independently standardized. Temporal derivatives were accurately extracted, exhibiting low Root Mean Square Errors (e.g., 0.0031 for in smooth regions and 0.0066 in sharp gradient regions). Benchmarking confirmed the spectral method's marginal superiority in satisfying the incompressibility constraint (RMS error 0.5914 vs. WENO5's 0.5934), leading to its selection for all spatial derivative computations. A preliminary kinematic viscosity estimate was very small (). Crucially, momentum residual analysis revealed significant, spatially coherent residuals, unequivocally confirming the dominant role of an unmeasured pressure gradient and validating the necessity of its proxy terms in the feature library. The resulting 26-term, independently standardized feature matrix and target temporal derivative vectors are now optimally conditioned for subsequent sparse regression-based equation discovery.