George MacKerron: code blog

GIS, software development, and other snippets

Archive for the ‘PL/R’ Category

Sunrise and sunset times with PostGIS and PL/R

with one comment

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).

Written by George

October 15th, 2012 at 5:52 pm

Posted in GIS,PL/R,PostGIS,SQL