Monitoring spatio-temporal changes in terrestrial gross primary productivity (GPP) of crops is key for estimating, understanding and predicting global carbon fluxes. Satellite remote sensing has been widely applied in the last decades to monitor agricultural resources, and the amount and quality of remote sensing data continuously increase. Since recently, and partly due the European Copernicus Programme, an unprecedented amount of open access data suitable for agriculture observations is now available. Benefiting from recent developments in satellite remote sensing technology, great advances in machine learning and advancements in our understanding of photosynthetic processes leading to increasingly complex and detailed photosynthesis models, we developed a hybrid approach to model GPP using satellite reflectance data by combining radiative transfer modeling and machine learning (ML). We have combined process-based model SCOPE with ML algorithms to estimate GPP of C3 crops using a variety of satellite data (Sentinel-2, Landsat and MODIS) and ancillary meteorological information. We link reflectance and meteorological data directly with crop GPP, bypassing the need of retrieving the set of input vegetation parameters needed to represent photosynthesis in an intermediate step, while still accounting for the complex processes of the original model. Several ML models, trained with the simulated data, were tested and validated using flux tower data. First, we tested our approach using Sentinel-2 data, which provide high frequency of observation, high spatial resolution of 20 m and multiple bands including red edge. Our final neural network model was able to estimate GPP at the tested flux towers with r2 of 0.92 and RMSE of 1.38 gC m$^{-2}$ d-1. Our model successfully estimated GPP across a variety of C3 crop types and environmental conditions, including periods of no vegetation, even tough it did not use any additional local information from the site. Since our learning approach is fast and efficient in the test phase and, at the same time, is based on a process-based model (and not on local empirical relationships), it can be applied globally. Furthermore, the simulated training dataset can be easily adapted to band settings of different instruments, assuring thus consistency among many sensors. However, such a global application requires high computational power and therefore we applied our approach to Landsat and MODIS data using Google Earth Engine (GEE) platform that provides cloud computing resources for processing large geospatial datasets. The results were validated using the FLUXNET2015 Dataset.