Sunrise and sunset times with PostGIS and PL/R
Since light can affect happiness, two important pieces of environmental data I add to the Mappiness data set during analysis are: (1) whether a response was made in daylight; and (2) day length when and where the response was made.
To derive these, I need the date, time and location of the response, and sunrise and sunset times for that date and location. R’s StreamMetabolism library provides sunrise/sunset calculations based on NOAA routines. And since my data is in PostGIS, it’s handy to use PL/R to access these.
I set this up as follows on my Ubuntu 12.04 GIS box.
In bash:
sudo aptitude install postgresql-9.1-plr sudo R |
In R:
install.packages('StreamMetabolism', dependencies = TRUE) |
In Postgres (9.1):
create extension plr; create table plr_modules (modseq int4, modsrc text); insert into plr_modules values (0, 'library("StreamMetabolism")'); create or replace function _r_daylight_period(lat double precision, lon double precision, datestring text, timezone text) returns setof integer as $$ return(as.integer(sunrise.set(lat, lon, datestring, timezone))) $$ language plr immutable strict; create or replace function sunrise(location geometry, moment timestamp without time zone, timezone text) returns timestamp without time zone as $$ select to_timestamp(min(s)) at time zone $3 from _r_daylight_period( st_y(st_transform($1, 4326)), -- 4326 = WGS84 st_x(st_transform($1, 4326)), to_char($2, 'YYYY/MM/DD'), $3 ) s $$ language sql immutable strict; create or replace function sunset(location geometry, moment timestamp without time zone, timezone text) returns timestamp without time zone as $$ select to_timestamp(max(s)) at time zone $3 from _r_daylight_period( st_y(st_transform($1, 4326)), -- 4326 = WGS84 st_x(st_transform($1, 4326)), to_char($2, 'YYYY/MM/DD'), $3 ) s $$ language sql immutable strict; create or replace function is_daylight(location geometry, moment timestamp without time zone, timezone text) returns boolean as $$ select (sunrise($1, $2, $3), sunset($1, $2, $3)) overlaps ($2, cast('0' as interval)); $$ language sql immutable strict; create or replace function hours_in_the_day(location geometry, moment timestamp without time zone, timezone text) returns double precision as $$ select extract(epoch from sunset($1, $2, $3) - sunrise($1, $2, $3)) / 3600.0; $$ language sql immutable strict; |
These new functions can be used like so:
select is_daylight(geom, moment, 'Europe/London'), hours_in_the_day(geom, moment, 'Europe/London') from mytable; |
(Note: what I actually do with the Mappiness data is to use a single call to _r_daylight_period, and calculate the other quantities as a second step. This speeds things up a lot, because I have millions of rows to process and Postgres doesn’t appear to do as much caching of immutable function results as one would like here).
-
Mike