For solving large-scale inconsistent linear systems with full column-rank coefficient matrices, we propose a maximum residual coordinate descent method integrated with k-means clustering, abbreviated as MRCD(k). The method begins by partitioning all equations (or equivalently, hyperplanes) into k clusters via k-means clustering, where the cosine distance between normal vectors of hyperplanes serves as the similarity metric. In each iteration of MRCD(k), a cluster is selected uniformly at random, and the approximate solution is updated using the equation with the maximum residual within the selected cluster. To enhance computational efficiency, we further extend this framework to a block variant, termed Block MRCD(k) or BMRCD(k), which updates the solution by simultaneously using the maximum-residual equation from each of the k clusters in every iteration. The convergence properties of both iteration methods are carefully studied. Numerical results demonstrate the robustness and the efficiency of the MRCD(k) and the BMRCD(k) iteration methods.