Seasonal weights of a time series: stl function versus subtracting average values

Say that we are dealing with a time series of 25 years of daily measurements of a variable with seasonal cycles, e.g. it reaches higher values in winter than in summer.

The objective is to compute the monthly weights that produce the seasonal oscillation. Let's say that the amplitude of the cycles does not seem to increase over time, so we assume this time series as being additive.

I can identify two approaches:

1. Use R's function stl, which performs a seasonal decomposition by LOESS, access the \$time.series[,"seasonal"] object to retrieve the weights, and then subtract these weights from the daily values belonging to the relevant month;
2. Compute the monthly average of the observed variable over the entire time series, and subtract them to each daily value of the relevant month.

I know the general consensus is to use option 1 because it has long been used and it has publications backing it up, but from a theoretical point of view, why should one use option 1? And roughly speaking, what is the real difference between the two options, especially considering that they can produce hugely different weights?