Protein solubility is critical across diverse biotechnological domains, including enzyme engineering and therapeutic protein production. However, conventional experimental techniques to determine solubility remain resource-intensive, time-consuming, and financially demanding. To overcome these limitations, we developed an innovative machine learning platform integrating geometry-aware graph neural networks (GNNs) with electrostatic field modeling. Our method captures intricate amino acid-level interactions governing solubility by harnessing detailed structural geometries alongside biophysical properties. Specifically, the incorporation of electrostatic field information allows our framework to accurately represent charge distributions and spatial interactions essential for solubility prediction. We rigorously curated a comprehensive dataset derived from diverse proteins expressed in Escherichia coli, providing a robust basis for training and validation of our framework. Benchmarking our model against independently assembled test datasets demonstrated superior performance, achieving over 82.4% accuracy and surpassing existing state-of-the-art approaches, particularly for proteins exhibiting low sequence similarity to the known training data. This capability to generalize effectively to novel and under-explored proteins addresses a longstanding challenge in computational protein solubility prediction. Furthermore, our approach holds significant practical implications, offering substantial reductions in experimental overhead by minimizing trial-and-error experimentation and accelerating protein development timelines. We expect our method to become an essential tool within experimental workflows, delivering an efficient, reliable, and cost-effective alternative to traditional assays and propelling progress across industrial and pharmaceutical protein engineering efforts.