This article aimed to investigate the potential of using three large‐scale precipitation products (PPs), including TRMM, ERA5‐Land, and MSWEP, and two land surface characteristics, including NDVI and Land Surface Temperature (LST), to improve the accuracy of an elevation‐based spatial non‐stationary, Locally Weighted Linear Regression (LWLR) method. A two‐step approach consisting of downscaling and merging was proposed and assessed across an orographically complex region in Iran to construct monthly precipitation maps for the 2001–2015 period. In the first step, three large‐scale PPs were downscaled to 1 km spatial resolution via the Area To Point Kriging (ATPK) method to address the spatial resolution discrepancy between the coarse pixels of PPs and pointwise gauge data. In the second step, five regression models with different combinations of predictors were developed and compared with the benchmark precipitation‐elevation regression model. The L2 regularization method was adopted to overcome the potential multicollinearity issue. The proposed framework could successfully diagnose the multicollinearity issues and dynamically estimate the regularization parameter in space and time. We utilized the holdout method to validate and compare the developed models, in which 50% of gauges were used for fitting the regression functions, and the remaining gauges were assigned to the validation dataset. The validation results showed that the overall monthly accuracy of the benchmark univariate elevation‐based model was improved where the MAE metric was decreased by 16.5%, and CC and KGE were increased by 3.4 and 5.8%, respectively. Also, the rBias was enhanced from −5.72 to 1.24%. Mean monthly precipitation maps of the 2001–2015 period generated by the model with the best combination of predictors and the benchmark model displayed the same spatial pattern consistent with the topography. Further analysis under different validation scenarios revealed that the contribution of PPs and auxiliary variables is the greatest when the training gauge network is sparse. Also, the developed multivariate model has a higher extrapolation ability and is more robust than the benchmark model.