(An EXPLAIN ANALYSE would be better here). Look at the expected number of stations
"Nested Loop (cost=0.00..994.94 rows=4046 width=4) (actual time=0.053..41.173 rows=78 loops=1)"
" Join Filter: ((6371.009::double precision * sqrt((pow(radians(((c.latitude_decimal - s.latitude_decimal))::double precision), 2::double precision) + (cos((radians(((c.latitude_decimal + s.latitude_decimal))::double precision) / 2::double precision)) * pow(radians(((c.longitude_decimal - s.longitude_decimal))::double precision), 2::double precision))))) <= 25::double precision)"
" -> Index Scan using city_pkey1 on city c (cost=0.00..6.27 rows=1 width=16) (actual time=0.014..0.016 rows=1 loops=1)"
" Index Cond: (id = 5182)"
" -> Seq Scan on station s (cost=0.00..321.08 rows=12138 width=20) (actual time=0.007..5.256 rows=12139 loops=1)"
" Filter: ((s.elevation >= 0) AND (s.elevation <= 3000))"
"Total runtime: 41.235 ms"
expects to have to touch a large proportion of the measurement table, therefore it thinks that it will be fastest to do a seq scan. In actual fact, for 78 stations, the index would be faster, but for 4046 it wouldn't.
This is rather unexpected. I'd have figured it would use the actual number.
If you will be querying by season quite regularly, had you considered partitioning by season?
Dave