Successful hydraulic infrastructure design necessitates information about future precipitation conditions, i.e., the return period for a specific event magnitude, which is why we need to treat rainfall as a random variable controlled by a distribution law. Our research addresses the sampling uncertainty associated with parameter estimation of the GEV distribution, providing a comparison between at-site and regional approaches. Regional frequency analysis was performed in the region of Thessaly, Greece, using annual maximum rainfall data from 55 rain gauges. To identify statistically homogeneous regions, four spatial covariates (rain gauge elevation, mean annual precipitation, latitude, and longitude) were used as inputs to the principal component analysis followed by the k-means clustering algorithm. We then calculated the parameters and their confidence intervals using two distinct methods. For the regional estimates, we used Bayesian MCMC sampling using the metropolis algorithm with non-informative priors, while for the at-site estimates, we used a typical L-moments procedure along with bootstrap resampling. Lastly, we compared the results to identify relationships with the aforementioned covariates. The results indicate that the gauge elevation and the mean annual precipitation are the best covariates. The location and scale parameters show a strong correlation with these covariates, while the shape parameter does not. The study procedure demonstrates the sampling uncertainty of the regional frequency analysis method for precipitation maxima. The results will be useful for hydrologic studies in the area.