On Wednesday I reported the results of my analysis examining the average date of first snow recorded at the Fairbanks Airport weather station. It was based on the snow_flag boolean field in the ISD database. In that post I mentioned that examining snow depth data might show the date on which permanent snow (snow that lasts all winter) first falls in Fairbanks. I’m calling this the first “true” snowfall of the season.
For this analysis I looked at the snow depth field in the ISD database for the Fairbanks station. The data was present for the years between 1973 and 1999, but isn’t in the database before that date. I’m not sure why it’s not in there after 1999, but luckily I’ve been collecting and archiving the data in the Fairbanks Daily Climate Summary (which includes a snow depth measurement) since late 2000. Combining those two data sets, I’ve got data for 27 years.
The SQL query I came up with to get the data from the data sets is a good estimate of what we’re interested in, but isn’t perfect because it only finds the date of first snow that lasts at least a week. In a place like Fairbanks where the turn to winter is so rapid and so dependent on the high albedo of snow cover, I think it’s close enough to the truth. Unfortunately, the query is brutally slow because it involves six (!) inner self-joins. The idea is to join the table containing snow depth data against itself, incrementing the date by one day at each join. The result set before the WHERE statement is the data for each date, plus the data for the six days following that date. The WHERE clause requires that snow depth on all those seven dates is above zero. This large query is a subquery of the main query which selects the earliest date found in each year.
There must be a better way to deal with conditions like this where we’re interested in the consecutive nature of the phenomenon, but I couldn’t figure out any other way to handle it in SQL, so here it is:
SELECT year, min(date) FROM ( SELECT extract(year from a.dt) AS year, to_char(extract(month from a.dt), '00') || '-' || ltrim(to_char(extract(day from a.dt), '00')) AS date FROM isd_daily AS a INNER JOIN isd_daily AS b ON a.isd_id=b.isd_id AND a.dt=b.dt - interval '1 day' INNER JOIN isd_daily AS c ON a.isd_id=c.isd_id AND a.dt=c.dt - interval '2 days' INNER JOIN isd_daily AS d ON a.isd_id=d.isd_id AND a.dt=d.dt - interval '3 day' INNER JOIN isd_daily AS e ON a.isd_id=e.isd_id AND a.dt=e.dt - interval '4 day' INNER JOIN isd_daily AS f ON a.isd_id=f.isd_id AND a.dt=f.dt - interval '5 day' INNER JOIN isd_daily AS g ON a.isd_id=g.isd_id AND a.dt=g.dt - interval '6 day' WHERE a.isd_id = '702610-26411' AND a.snow_depth > 0 AND b.snow_depth > 0 AND c.snow_depth > 0 AND d.snow_depth > 0 AND e.snow_depth > 0 AND f.snow_depth > 0 AND g.snow_depth > 0 AND extract(month from a.dt) > 7 ) AS snow_depth_conseq GROUP BY year ORDER BY year;
See what I mean? It’s pretty ugly. Running the result through the same R script as in my previous snowfall post yields this plot:
Between 1973 and 2008 we’ve gotten snow lasting the whole winter starting as early as September 12th (that was the infamous 1992), and as late as the first of November (1976). The median date is October 13th, which matches my impression. Now that the leaves have largely fallen off the trees, I’m hoping we get our first true snowfall on the early end of the distribution. We’ve still got a few things to take care of (a couple new dog houses, insulating the repaired septic line, etc.), but once those are done, I’m ready for the Creek to freeze and snow to blanket the trails.
We got our first dusting of snow last night. It stuck around until after noon, allowing me to take the photo on the right when I went for a walk with Nika around the peat bog. You can really tell where the permafrost is by the thick layer of insulating moss that keeps the ground frozen, and is keeping the snow from melting in the photo.
Every year when the first snow falls it seems like it’s earlier than the last, and there’s usually some discussion at the office about how short the summer turned out to be. The early snows of 1992 that knocked out power for days all over town are also normally mentioned. I decided to look and see if I had some data that could place this year’s first snowfall in a historical context.
One of the few free† long-term weather datasets that’s available from the National Climate Data Center is the Integrated Surface Dataset (ISD), which contains daily weather observations for more than 20,000 stations. The Fairbanks Airport station has been in operation for more than 100 years, but it moved in 1946, so I only used data from 1946–2008. In addition to a series of numerical observations (minimum and maximum temperature, pressure, wind speed, etc.), the dataset contains several fields used to indicate whether a particular phenomenon was observed during that day. One of them, snow_flag, is defined as: “True indicates there was a report of snow or ice pellets during the day.”
That’s perfect. Snow depth is another parameter I considered, but this data wasn’t collected until the mid-70s, and it doesn’t really help us answer the question because most of the time the first snowfall of the year doesn’t last long enough to be recorded as snow on the ground.
Here’s the SQL query to find the earliest snowfall date for each year for the Fairbanks Airport station:
SELECT year, min(date) FROM ( SELECT extract(year from dt) AS year, to_char(extract(month from dt), '00') || '-' || ltrim(to_char(extract(day from dt), '00')) AS date, snow_flag FROM isd_daily WHERE isd_id = '702610-26411' AND extract(month from dt) > 7 AND snow_flag = 't' ) AS snow_flag_sub GROUP BY year ORDER BY year;
Mix in a little R:
fs <- read.table("first_snow_mm-dd", header=TRUE, row.names=1) fs$date<-as.Date(fs$date, "%m-%d") png("first_snow_mm-dd.png", height=500, width=500, units="px", pointsize=12) hist(fs$date, breaks="weeks", labels=FALSE, xlab="Date of first snowfall", main="First snowfall reported, Fairbanks Airport (PAFA) station", plot=TRUE, freq=TRUE, ylim=c(0, 20), col="gray60") text(as.Date("2009-09-23"), 19, "⇦ 2009", srt=90, col="darkred") dev.off()
And you get this plot:
You can see from the plot that the first snowfall comes somewhere between August 3rd and October 26th, with the week of September 21st being the most common. So we’re right on schedule this year.
Another analysis that I’ve been meaning to do is to find the average date when the snow that falls lasts the entire winter. Since I’ve been in Fairbanks, my estimate of this date is the second week of October, but I’ve never actually looked it up to see if that’s true or not. Unfortunately, this requires good snow depth data, and the ISD dataset doesn’t have snow depth for Fairbanks prior to 1975. It’s also a bit more complicated than looking for the earliest snow_flag = 't' because you need to examine future rows to know if the snow depth observation you’re examining lasted more than a few days.
†Why isn’t all the data collected by the Weather Service freely available? Public money was used to collect, analyze, and archive it, so I think it should be made available to the public that paid for it.
In my previous post on my weather station I included a graph of the relationship between the three temperature sensors I’m collecting data from. It turns out that the plot is incorrect because I had messed up the program that converts the raw data from one of my Arduino stations. What became clear after looking at the (fixed!) data is that the enclosure I built isn’t adequate at keeping the sensors cool when in direct sunlight. I could probably improve the design of the enclosure somewhat, but I decided to aspirate the sensors instead.
The photo on the right shows the inside of the enclosure. The sensors are sitting on a platform inside the pipe, and there’s a small muffin fan (the kind you’d use to vent a computer case) on the top. It’s a 12V fan, and at the moment it’s being driven by a 9V AC/DC converter. The plan is to replace the converter with a 12V solar cell that is sufficient to drive the fan. This way it doesn’t consume any electricity, and the fan is only spinning when it’s necessary (when the sun is out). Thus far, I haven’t found a suitable solar panel. The small ones designed to charge a cell phone battery don’t operate at the correct voltage (and probably don’t produce enough current anyway), and the big ones designed to keep a car battery charged are expensive and overpowered for my needs. With winter rapidly approaching, I’ve got plenty of time to figure something out.
The pipe is a piece of 4” sewer pipe that’s been spray painted white and has a series of holes drilled into the bottom. The fan pulls air up through these holes and over the sensor array in the middle. If I had it to do over again, I’d cut the pipe a bit shorter so it’s not so difficult to get into the enclosure. But for keeping the sensors bathed in atmospheric air, it works quite well.
The plot shows the observed temperature for each of the three sensors I’ve got. The blue line is the “west” sensor that’s inside the enclosure I built and is the subject of this post. The red line is the reported temperature from the Rainwise station that sits atop a post attached to the dog yard gate. The green line is the sensor that’s behind the house and under the oil tank. I built the enclosure for the west sensors on July 12th and installed it that evening. You can immediately see the effect of the shielding during the high temperature peak on the 13th. But you can still see the two little peaks that are present in the previous plots. These peaks come from direct sun on the station, split by some trees that shade the station for an hour or two.
On the evening of the 13th I installed the pipe and fan. It was smoky on the 14th and cloudy on the 15th, but you can see the effect of the fan on the following dates. The double peaks are now gone, and the temperature from the west sensors at the high point during the day is now a few degrees cooler than the measurement from the Rainwise station. Also notice that all three sensors are virtually identical on the 15th when it was cloudy and raining.
What this demonstrates to me is that the aspirated west sensor is now the best reference sensor for our site. The Rainwise sensor is a close second, but it’s Gill multi-plate radiation shield isn’t as effective as my aspiration system at reducing the effect of solar heating on the station.
Yesterday I looked at how wind might be affecting my bicycling to and from work. Today I’ll examine the idea that Miller Hill is confounding the effect of wind on average speed by excluding this portion of the trip from the analysis. To do this, I include a bounding box comparison in the SQL statement that extracts the wind factors for track points. The additional WHERE condition looks like this:
ST_Within(point_utm, ST_SetSRID(ST_MakeBox2D(ST_Point(454861,7193973), ST_Point(458232,7199159)), 32606))
The same ST_Within test is used in the calculation of average speed for each of the trips from work to home. After compiling the wind factors and average speeds, we compare the two using R. Here are the updated results:
lm(formula = mph ~ wind, data = data) Residuals: Min 1Q Median 3Q Max -1.87808 -0.55299 0.04038 0.62790 1.19076 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 16.8544 0.2176 77.442 <2e-16 *** wind 0.3896 0.2002 1.946 0.0683 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9445 on 17 degrees of freedom Multiple R-squared: 0.1822, Adjusted R-squared: 0.1341 F-statistic: 3.788 on 1 and 17 DF, p-value: 0.06834
This time around the model and both coefficients are statistically significant (finally!), and “wind factor” is positively correlated with my average speed over the part of the route that doesn’t include Miller Hill and Railroad drive. It’s not a major contributor, but it does explain approximately 18% of the variation in average speed.
I decided to look at wind a little more deeply after yesterday’s bike ride home. It seemed clear to me that the wind was strongly at my back for much of the route. It wasn’t my fastest ride home, but it was close, and it didn’t feel like I was working all that hard.
Here’s the process. First, examine all my bicycling tracks individually, using PostGIS’s ST_Azimuth function to calculate the direction I was traveling at each point. The query uses another of the new window functions (lead) in PostgreSQL 8.4.
SELECT point_id, dt_local, ST_Azimuth( point_utm, lead(point_utm) OVER (PARTITION BY tid ORDER BY dt_local) ) / (2 * pi()) *360 FROM points WHERE tid = TID ORDER BY dt_local;
Then, for each point, find the direction the wind was blowing. This is a pretty slow query, but I haven’t found a better way to compare timestamps in the database to find the closest record. This technique, based on converting both timestamps to “epoch,” which is the number of seconds since January 1st, 1970, is faster than using an interval type of operation (like: WHERE obs_dt - POINT_DT BETWEEN interval '-3 minutes' AND interval '3 minutes').
SELECT obs_dt, wdir, wspd FROM observations WHERE abs(extract(epoch from obs_dt) - extract(epoch from POINT_DT)) < 5 * 60 AND wspd IS NOT NULL AND wdir IS NOT NULL ORDER BY abs(extract(epoch from obs_dt) - extract(epoch from POINT_DT)) LIMIT 1;
Now I’ve got the direction I was traveling and the direction the wind is coming from. I wrote a Python function that returns a value from –1 (wind is in my face) to 1 (wind is at my back). The procedure is to convert the wind directions to unit u and v vectors and get the distance between the endpoints of each vector. The distances are then scaled such that wind behind the direction traveled range from 0 – 1, and from –1 – 0 for wind blowing against the direction traveled.
def wind_effect(mydir, winddir): """ Returns a number from 1 (wind at my back) to -1 (wind in my face) based on the directions passed in. Remember that wind direction is where the wind is *from*, so a wind direction of 0 and a mydir of 0 means the wind is in my face. """ try: mydir = float(mydir) winddir = float(winddir) except: return(None) my_spd = 1.0 wind_spd = 1.0 u_mydir = -1 * my_spd * math.sin(math.radians(mydir)) v_mydir = -1 * my_spd * math.cos(math.radians(mydir)) u_winddir = -1 * wind_spd * math.sin(math.radians(winddir)) v_winddir = -1 * wind_spd * math.cos(math.radians(winddir)) distance = math.sqrt((u_mydir - u_winddir)**2 + (v_mydir - v_winddir)**2) factor = (1.41421356 - distance) if factor < 0.0: factor = factor / -0.58578644 else: factor = factor / -1.41421356 return(factor)
Finally, multiply this value by the wind speed at that time, and sum all these values for an entire bicycling track. The result is a “wind factor.” A positive wind factor means the wind was generally at my back during the ride, negative means it was blowing in my face. Yesterday’s ride home had the highest wind factor (1.07) among trips since June. So the wind really was at my back!
Can “wind factor” help predict average speed? Here’s the R and results:
$ R --save < wind_from_abr.R > data<-read.table('wind_factor_from_abr',header=TRUE) > model<-lm(speed ~ wind, data) > summary(model) Call: lm(formula = speed ~ wind, data = data) Residuals: Min 1Q Median 3Q Max -0.90395 -0.46782 -0.04334 0.40286 0.85918 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 14.7796 0.1471 100.48 <2e-16 *** wind 0.4369 0.2875 1.52 0.147 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.5522 on 17 degrees of freedom Multiple R-squared: 0.1196, Adjusted R-squared: 0.06784 F-statistic: 2.31 on 1 and 17 DF, p-value: 0.1469
Hmm. Not a whole lot of help here. The model is close to being statistically significant (although it’s not…), and it’s not very predictive (only 12% of the variation in average speed is explained by wind factor). However, the directionality of the (not quite statistically significant) wind coefficient is correct. A positive wind factor is (weakly) correlated with a higher average speed.
Thinking more about my route from work, I suspect that the route is actually two trips: the trip from ABR to the bottom of Miller Hill (4.8 miles) and the two mile trip over Miller Hill to our house. I’ll bet that wind becomes statistically significant if I only consider the first part of the trip: wind doesn’t have as much effect on a hill climb, and after making it over the top, the rest is a bumpy, gravel road where speed is determined more by safety than wind or how hard I’m pedalling. I think this might also resolve the question of why the ride home is so much easier than to work. It’s not because I’m glad to be out of work or because I’m carrying a lunchbox full of food to work, it’s because it’s downhill from ABR to the bottom of Miller Hill.