Land use and land cover (LULC) maps, providing crucial information for monitoring the Earth's surface, are one of the most essential products for numerous studies. Using only the spectral information in the classification process might cause poor performances in the areas with heterogeneous landscape characteristics. To overcome this problem, auxiliary and ancillary data are usually employed to improve classification accuracy. The objective of this study is to integrate auxiliary data (topographic and climatic features) and ancillary data (spectral indices and texture measures) into spectral bands of Sentinel-2A imagery and evaluate the performances of advanced feature selection methods. In this context, genetic algorithm-based random forest (GA-RF), HSIC-Lasso, and Relief-F feature selection approaches were utilized to determine the most informative features for the classification process from a high-dimensional dataset consisting of 102 features. Whilst the GA-RF algorithm selected 65 features, HSIC-Lasso chose 38 features, and Relief-F determined 51 features as ideal subsets. These feature subsets together with the whole data were inputted into a supervised classification process using the random forest (RF) classifier, whose parameters were optimized using random search algorithm. The highest overall accuracy of the produced thematic maps was estimated as 91.05% for the subset determined by the HSIC-Lasso algorithm, which was also the fastest algorithm (5.71 s). McNemar's statistical significance test confirmed the superiority of the HSIC-Lasso method over the GA-RF and Relief-F algorithms. SHapley Additive exPlanations method was also applied to analyze the relative importance of a feature according to the model output.