Think of the angle $y$ at any time $t$ as the accumulation of small changes in the angle. Symbolically, when $f(t)$ is the rate of change of the angle at time $t$ and $t_0$ is the beginning of the observations,
$$y(t) = y(t_0) + \int_{t_0}^t f(t)\,\mathrm{d}t.$$
Your problem is that $y(t)$ has been recorded modulo $360$ degrees -- perhaps with some error $\epsilon(t).$ That is, you have observed only the values
$$y^{*}(t) = y(t) + \epsilon(t) \mod 360.$$
You can, however, reconstruct $y(t) + \epsilon(t)$ provided you have sufficiently frequent observations. For successive times $t \lt s,$ notice
$$y^{*}(s) - y^{*}(t) = y(s) - y(t) + \epsilon(s) - \epsilon(t) \mod 360 = \int_t^s f(t)\,\mathrm{d}t + \delta$$
where $\delta$ equals the contribution of the errors $\epsilon(s)-\epsilon(t)$ plus, perhaps, some integral multiple of $360$ whenever there has been an angular break between $y^{*}(t)$ and $y^{*}(s).$ Now, provided the size of the total error $|\epsilon(s)-\epsilon(t)|$ is less than $180$ degrees and provided the angle didn't go around more than once, we can figure out whether a break occurred: if $|\epsilon(s)-\epsilon(t)| \gt 180,$ add or subtract $360$ degrees from $\delta$ to place it into the interval from $-180$ to $+180$ degrees.
Although we can't observe these errors directly, if we are sampling frequently enough to make the increments $y(t_i) - y(t_{i-1})$ fairly small, we simply apply this adjustment to the observed differences. Thus,
Whenever $|y^{*}(s)-y^{*}(t)| \gt 180,$ add or subtract $360$ degrees from $\delta$ to place it into the interval from $-180$ to $+180$ degrees.
Equivalently, compute the differences modulo $180$ but express them in the range from $-180$ to $+180$ degrees rather than (as is conventional) the range from $0$ to $360.$
Let's call the adjusted value $\delta^{*}(t,s),$ so that
$$y^{*}(s) - y^{*}(t) = \int_t^s f(t)\,\mathrm{d}t + \delta(t,s)^{*}.$$
This is equality, not equality modulo $360.$ We may now remove the effect of recording the angles modulo $360$ by summing these adjusted differences. When the observations are made at times $t_0 \lt t_1\lt \cdots \lt t_n,$ we have
$$\begin{aligned}
y^{*}(t_i) &= y^{*}(t_0) + \left[y^{*}(t_1) - y^{*}(t_0)\right] + \cdots + \left[y^{*}(t_i) - y^{*}(t_{i-1})\right] \\
&=y(t_0) + \int_{t_0}^{t_i} f(t)\,\mathrm{d}t + \delta(t_0,t_1)^{*} + \delta(t_1,t_2)^{*} + \cdots + \delta(t_{i-1},t_i)^{*} \\
&= y(t_i) + \left[\epsilon(t_i) - \epsilon(t_0)\right].
\end{aligned}$$
The problem with computation modulo $360$ is gone: you may now use any procedure you like to model the response $y^{*}(t).$
Here is an illustration with a fairly difficult dataset. The data were generated according to the model $y(t) = 30t \mod 360$ and observed annually from 1980 through 2020 with iid Normally distributed error of standard deviation $60$ degrees (a large amount).
The trend is barely discernible in the raw data, but the angle adjustment algorithm has visibly aligned them. We may fit a least squares model to the adjusted data, for instance, producing this result:
The expanded vertical scale for the raw data shows details of the fit and their deviations from it. Incidentally, in this example the estimate of the slope is $28.0 \pm 0.74$ degrees, not remarkably different from the true value of $30$ degrees (the p-value for this comparison is $1.1\%$).
I will end by remarking that when the standard deviation of the errors $\epsilon(t)$ is large (greater than $180/2/\sqrt{2} \approx 64$ degrees, roughly), sometimes the angular adjustment will be incorrect. This will show up in the model residuals as a sudden change by a value around 360 degrees. Thus, a routine analysis of the model residuals can detect such problems, enabling you to modify the adjusted values for a better fit. The details of this will depend on your model and the fitting procedure.
This R code created the figures. At "adjust the angles" it shows how the angle adjustment can be computed efficiently.
#
# Specify the data-generation process.
#
year <- 1980:2020 # Dates to use
beta <- 30 # Annual rate of change
sigma <- 60 # Error S.D.
#
# Generate the data.
#
set.seed(17)
angle <- (year * beta + rnorm(length(year), 0, sigma)) %% 360
X <- data.frame(year, angle)
#
# Adjust the angles.
#
X$`total angle` <- with(X, {
d <- (diff(angle) + 180) %% 360 - 180
cumsum(c(angle[1], d))
})
#
# Fit a model to the adjusted angles.
#
fit <- lm(`total angle` ~ year, X)
#
# Analyze the fit.
#
b <- coefficients(fit)
y.hat <- predict(fit)
#--Compute dates the fit must wrap around from 360 to 0:
y.breaks <- seq(floor(min(y.hat) / 360)*360, max(y.hat), by=360)
year.breaks <- (y.breaks - b[1]) / b[2]
#--Make the plots:
u <- ceiling(max(X$`total angle`)/360)
par(mfcol=c(1,2))
#--The fits:
plot(X$year, X$angle, pch=19, ylim=c(0, 360), yaxp=c(0, 360, 4),
col="gray", ylab="Angle (degrees)", xlab="Year",
main="Raw Data and Fit")
for (x in year.breaks)
abline(c(-x * b[2], b[2]), col="Red", lwd=2)
plot(X$year, X$`total angle`, ylim=c(0,u*360), yaxp=c(0, u*360, u),
xlab="Year", ylab="Total angle",
main="Adjusted Data and Fit")
abline(fit, col="Red", lwd=2)
#--The raw data:
plot(X$year, X$angle, ylim=c(0,u*360), yaxp=c(0, u*360, u),
pch=19, col="gray", ylab="Angle (degrees)", xlab="Year",
main="Raw Data")
plot(X$year, X$`total angle`, ylim=c(0,u*360),
yaxp=c(0, u*360, u),
xlab="Year", ylab="Total angle",
main="Adjusted Data")
par(mfcol=c(1,1))
