Seasons are circular (the 31st of December is followed by the first of january), so tests and graphs for circularity are probably the best first line of attack on the problem.
The most basic test is Rao's test, which tests for a uniform distribution of the counts around the seasons. It is available as part of the -circular- package for Stata, written by Nick Cox.
There is also an R package, also called -circular-, which implements Rao's test.
I don't know anything about SPSS. Or, at least, nothing after 1981, when I stopped using it.
Are you taking season as categorical (typically 4 seasons) or as continuous? If continuous, then you need daily or weekly resolution, perhaps monthly. If continuous, then circular stats (interval scale), but beware of normal error model. If you are taking season as categorical, then nominal scale, one season per level of categorical variable. For counts per unit, start with Poison error model, look at residuals for homogeneity. If not homogeneous, try quasi Poisson. If residuals not homogeneous, try negative binomial error model. If still not homogeneous, then consider zero inflated model or hurdle model. There is much to be said for getting the error model right, starting with getting the parameter estimates right.
I stopped using SPSS in 1984, I used SAS from 1984 to 1992. Now I use R and sometimes SPSS (now that it has got its act together) or Minitab.
You won't need any special package. The function you need is "glm", like
model = glm(count~season, family = poisson(link = "log"))
With the data is overdispersed (variance/mean >> 1) you have the possibility to fit a quasi-poisson model with family=quasipossion or to switch to a negative binomial model.
The function to fit a negative binomial model is glm.nb and comes with the MASS package.
This is brilliant thank you very much! Actually it gives me very similar results to a straight chi-squared, but it is very useful to check through the workings.
The tests in a Poisson model *are* chi²-tests (the squared z values for individial coefficients are nothing but chi² values, like the squared t-values from standard linear models are nothing but F-values [known from ANOVA]), so the results of the tests should be the same (up to some rounding errors because of the different algorithms used). The Poisson model gives you effect estimates, it can deal with continuous predictors, can account for other variables/confounders, in any combination of continuous and categorical scales, and it can fit interactions. All this is not the case for the "simple" chi²-test.