Soil nutrient estimation can be used as a key input to increase crop yield and agriculture fertilization. Due to the shortage of the ground measured spectrum technology and costliness to obtain hyperspectral images, multispectral remote sensing data is used to explore the soil nutrient content estimation. In this paper, we have used optical remote sensing data (Landsat-8 and Sentinel-2), terrain/climate data (precipitation, radiation, slope etc.) and ground truth value to estimate four nutrients: N, K, P, and OC for two districts of Maharashtra, India. We compared four linear and non-linear regression models: multiple linear regression (MLR), random forest regression (RFR), support vector machine for regression (SVR) and gradient boosting (GB) for estimation of NPK and OC. Comparative results suggest that, GB and RFR performed better than other models with sMAPE in range of 0.125-0.377 for all nutrients, which is better or comparable with literature reported accuracy . Therefore, the approach has potential to generate high resolution (< ha) soil nutrients map and can reduce soil sampling effort/cost.