Figuring Out Sunrise, Sunset, Dawn and Dusk from a Time Series at a Location

Example of crawl Package (Johnson et al. 2008) Applied to Culebra VPS Location Data

Step 1: Read in Culebra sample sequence

getwd()  # where am I now?
## [1] "C:/Users/jack/Dropbox/Eco691NM_GraphsAndNetworks - Working/CulebraBonefish/B30414_Seq1_07202012"
# Set the working directory to wherever the file is (change to match your
# computer):
setwd("C:\\Users\\jack\\Dropbox\\Eco691NM_GraphsAndNetworks - Working\\CulebraBonefish\\B30414_Seq1_07202012")
dir()  # look at the files here
##  [1] "CTCRW_Results_PGDB.mdb"                            
##  [2] "schema.ini"                                        
##  [3] "SunriseSunset.html"                                
##  [4] "SunriseSunset.md"                                  
##  [5] "SunriseSunset.Rmd"                                 
##  [6] "Timezones.r"                                       
##  [7] "TrueAndSimulatedTracks.png"                        
##  [8] "VPS_Positions.txt"                                 
##  [9] "VPSData_30414_.png"                                
## [10] "VPSData_30414_High_HPE_Example.png"                
## [11] "VPSData_30414_Seq1_07202012.png"                   
## [12] "VPSData_30414_Seq1_07202012_P4_All.png"            
## [13] "VPSData_30414_Seq1_07202012_TruePath.png"          
## [14] "VPSData_30414_Seq1_07202012_TruePath_P4_Detail.png"
b30414_Seq1 = read.table("VPS_Positions.txt", head = TRUE, sep = ",")
str(b30414_Seq1)
## 'data.frame':    99 obs. of  17 variables:
##  $ OBJECTID   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ TRANSMITTER: int  30414 30414 30414 30414 30414 30414 30414 30414 30414 30414 ...
##  $ DETECTEDID : Factor w/ 1 level "A69-1601-30414": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DATETIME_  : Factor w/ 99 levels "7/20/2012 2:11:46",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ X          : num  986 955 973 968 958 ...
##  $ Y          : num  1147 1010 996 981 971 ...
##  $ D          : num  1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 ...
##  $ LAT        : num  18.3 18.3 18.3 18.3 18.3 ...
##  $ LON        : num  -65.3 -65.3 -65.3 -65.3 -65.3 ...
##  $ n          : int  3 4 2 7 3 1 10 1 10 1 ...
##  $ HPE        : num  22.2 7.4 7.1 7.5 9.4 7.7 7.3 9.8 18.5 18.7 ...
##  $ HPEm       : logi  NA NA NA NA NA NA ...
##  $ TEMP       : logi  NA NA NA NA NA NA ...
##  $ DEPTH      : logi  NA NA NA NA NA NA ...
##  $ ACCEL      : logi  NA NA NA NA NA NA ...
##  $ DRX        : Factor w/ 82 levels "VPS02 VPS05 VPS06 VPS11",..: 4 62 34 21 26 75 63 82 5 82 ...
##  $ URX        : Factor w/ 82 levels "VPS02 VPS05 VPS06 VPS11 VPS12",..: 3 62 35 22 27 75 63 82 4 82 ...

We have to change DATETIME_ from factor to POSIXct

A digression on time:
First we are going to look at a single time and get the format right.

(UTCTime_ct = as.POSIXct(b30414_Seq1$DATETIME_[1], tz = "GMT", format = "%m/%d/%Y %H:%M:%S"))
## [1] "2012-07-20 02:11:46 GMT"
# There is a 4 hour offset to get local PR time PR does not have daylight
# savings time!
UTCTime_ct  # POSIXct date-time
## [1] "2012-07-20 02:11:46 GMT"
as.numeric(UTCTime_ct)  # The number of seconds since 1970/1/1 00:00:00 !
## [1] 1.343e+09
# This is how R stores the POSIXct!  POSIXlt is stored as a vector
UTCtimelt = as.POSIXlt(UTCTime_ct)  # POSIXlt date-time
UTCtimelt
## [1] "2012-07-20 02:11:46 GMT"
str(UTCTime_ct)
##  POSIXct[1:1], format: "2012-07-20 02:11:46"
str(UTCtimelt)
##  POSIXlt[1:1], format: "2012-07-20 02:11:46"
# POSIXct is a number but POSIXlt is a vector of a lot of pieces of info.
# For example:
UTCtimelt$year + 1900  # Year
## [1] 2012
UTCtimelt$mon + 1  # month
## [1] 7
UTCtimelt$yday + 1  # day of the year
## [1] 202
# etc.

# Local time?
format(UTCTime_ct, tz = "America/Puerto_Rico")  # What's the local time? 
## [1] "2012-07-19 22:11:46"
# (Atlantic Standard Time or in R 'America/Puerto_Rico')

Now do it for the whole dataset

b30414_Seq1$UTCTime_ct = as.POSIXct(b30414_Seq1$DATETIME_, tz = "GMT", format = "%m/%d/%Y %H:%M:%S")
b30414_Seq1$AST_chr = format(b30414_Seq1$UTCTime_ct, tz = "America/Puerto_Rico")
b30414_Seq1$AST_ct = as.POSIXct(b30414_Seq1$AST_chr, tz = "America/Puerto_Rico", 
    format = "%F %T")
b30414_Seq1$AST_lt = as.POSIXlt(b30414_Seq1$AST_chr, tz = "America/Puerto_Rico", 
    format = "%F %T")

# Note that AST_ct is also a POSIXct
str(b30414_Seq1)
## 'data.frame':    99 obs. of  21 variables:
##  $ OBJECTID   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ TRANSMITTER: int  30414 30414 30414 30414 30414 30414 30414 30414 30414 30414 ...
##  $ DETECTEDID : Factor w/ 1 level "A69-1601-30414": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DATETIME_  : Factor w/ 99 levels "7/20/2012 2:11:46",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ X          : num  986 955 973 968 958 ...
##  $ Y          : num  1147 1010 996 981 971 ...
##  $ D          : num  1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 ...
##  $ LAT        : num  18.3 18.3 18.3 18.3 18.3 ...
##  $ LON        : num  -65.3 -65.3 -65.3 -65.3 -65.3 ...
##  $ n          : int  3 4 2 7 3 1 10 1 10 1 ...
##  $ HPE        : num  22.2 7.4 7.1 7.5 9.4 7.7 7.3 9.8 18.5 18.7 ...
##  $ HPEm       : logi  NA NA NA NA NA NA ...
##  $ TEMP       : logi  NA NA NA NA NA NA ...
##  $ DEPTH      : logi  NA NA NA NA NA NA ...
##  $ ACCEL      : logi  NA NA NA NA NA NA ...
##  $ DRX        : Factor w/ 82 levels "VPS02 VPS05 VPS06 VPS11",..: 4 62 34 21 26 75 63 82 5 82 ...
##  $ URX        : Factor w/ 82 levels "VPS02 VPS05 VPS06 VPS11 VPS12",..: 3 62 35 22 27 75 63 82 4 82 ...
##  $ UTCTime_ct : POSIXct, format: "2012-07-20 02:11:46" "2012-07-20 02:17:54" ...
##  $ AST_chr    : chr  "2012-07-19 22:11:46" "2012-07-19 22:17:54" "2012-07-19 22:20:26" "2012-07-19 22:21:27" ...
##  $ AST_ct     : POSIXct, format: "2012-07-19 22:11:46" "2012-07-19 22:17:54" ...
##  $ AST_lt     : POSIXlt, format: "2012-07-19 22:11:46" "2012-07-19 22:17:54" ...
head(b30414_Seq1)
##   OBJECTID TRANSMITTER     DETECTEDID         DATETIME_     X      Y   D
## 1        1       30414 A69-1601-30414 7/20/2012 2:11:46 986.4 1147.5 1.5
## 2        2       30414 A69-1601-30414 7/20/2012 2:17:54 954.8 1009.6 1.5
## 3        3       30414 A69-1601-30414 7/20/2012 2:20:26 972.6  996.1 1.5
## 4        4       30414 A69-1601-30414 7/20/2012 2:21:27 968.0  981.2 1.5
## 5        5       30414 A69-1601-30414 7/20/2012 2:22:39 957.9  970.7 1.5
## 6        6       30414 A69-1601-30414 7/20/2012 2:24:15 945.0  944.6 1.5
##    LAT    LON n  HPE HPEm TEMP DEPTH ACCEL
## 1 18.3 -65.25 3 22.2   NA   NA    NA    NA
## 2 18.3 -65.25 4  7.4   NA   NA    NA    NA
## 3 18.3 -65.25 2  7.1   NA   NA    NA    NA
## 4 18.3 -65.25 7  7.5   NA   NA    NA    NA
## 5 18.3 -65.25 3  9.4   NA   NA    NA    NA
## 6 18.3 -65.25 1  7.7   NA   NA    NA    NA
##                                   DRX                                 URX
## 1       VPS02 VPS05 VPS08 VPS09 VPS12       VPS02 VPS05 VPS08 VPS09 VPS12
## 2       VPS08 VPS09 VPS12 VPS15 VPS16       VPS08 VPS09 VPS12 VPS15 VPS16
## 3             VPS02 VPS12 VPS15 VPS16             VPS02 VPS12 VPS15 VPS16
## 4 VPS02 VPS08 VPS09 VPS12 VPS15 VPS16 VPS02 VPS08 VPS09 VPS12 VPS15 VPS16
## 5       VPS02 VPS08 VPS11 VPS12 VPS16       VPS02 VPS08 VPS11 VPS12 VPS16
## 6                   VPS09 VPS12 VPS17                   VPS09 VPS12 VPS17
##            UTCTime_ct             AST_chr              AST_ct
## 1 2012-07-20 02:11:46 2012-07-19 22:11:46 2012-07-19 22:11:46
## 2 2012-07-20 02:17:54 2012-07-19 22:17:54 2012-07-19 22:17:54
## 3 2012-07-20 02:20:26 2012-07-19 22:20:26 2012-07-19 22:20:26
## 4 2012-07-20 02:21:27 2012-07-19 22:21:27 2012-07-19 22:21:27
## 5 2012-07-20 02:22:39 2012-07-19 22:22:39 2012-07-19 22:22:39
## 6 2012-07-20 02:24:15 2012-07-19 22:24:15 2012-07-19 22:24:15
##                AST_lt
## 1 2012-07-19 22:11:46
## 2 2012-07-19 22:17:54
## 3 2012-07-19 22:20:26
## 4 2012-07-19 22:21:27
## 5 2012-07-19 22:22:39
## 6 2012-07-19 22:24:15
Variables of interest to us are:
    Transmitter - the tag we are listening to  
    DATETIME_   - date and time as a factor  
    X, Y        - Position in meters using a custom World Azimuthal Equidistant projection  
                  centered on the array. (Center of array is at X=1000, Y=1000).  
    D           - depth in meters.  (Set at 1.5m, not calculated, for this dataset)  
    LAT, LON    - Position in Polar coordinates ion decimal degrees (geographic)  
    HPE         - Horizontal Position Error in meters  
    DRX, URX    - the list of receivers that heard this ping  
    UTCTime_ct  - Same as DATETIME_ but converted to POSIXct format  
    AST_ct      - Local time.  Easier to compare with sunrise/sunset calculations  

Sunrise/Sunset Calculations:

An interactive Calculator is at http://www.esrl.noaa.gov/gmd/grad/solcalc/sunrise.html or a new one using Google Maps is at http://www.esrl.noaa.gov/gmd/grad/solcalc/ .
On 7/20/2012, apparent sunrise was at 5:55 local time (9:55 UTC) &
apparent sunset was at 19:00 local time (23:00 UTC).

To calculate for lots of points, the maptools package has several functions including sunriset (for sunrise or sunset) and crepuscule (for various definitions of dawn/dusk). The angle of the sun below the horizon (solarDep) determines which definition we are talking about:

         Definition name            solarDep (in degrees)
       Astronomical dawn/dusk:       18
       Nautical dawn/dusk:           12
       Civil dawn/dusk:               6
See the example in the help for these functions to see how to specify a
   location with Lat/Long and find the event you want for the day you want.

Make a variable called Light that is “dawn”, “day”, “dusk”, or “night” to use as a covariate in the movement model. Use the Civil dawn/dusk definition

require(rgeos)
## Loading required package: rgeos
## rgeos version: 0.3-4, (SVN revision 438)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE
require(maptools)
## Loading required package: maptools
## Loading required package: sp
## Checking rgeos availability: TRUE
######################################## The example from maptools Location of Helsinki, Finland, in decimal
######################################## degrees, as listed in NOAA's website
hels <- matrix(c(24.97, 60.17), nrow = 1)
Hels <- SpatialPoints(hels, proj4string = CRS("+proj=longlat +datum=WGS84"))
d041224 <- as.POSIXct("2004-12-24", tz = "EET")
## Civil dawn
crepuscule(Hels, d041224, solarDep = 6, direction = "dawn", POSIXct.out = TRUE)
##        day_frac                time
## newlon   0.3519 2004-12-24 08:26:46
## Civil dusk
crepuscule(Hels, d041224, solarDep = 6, direction = "dusk", POSIXct.out = TRUE)
##        day_frac                time
## newlon   0.6757 2004-12-24 16:13:00
## Sunrise
sunriset(Hels, d041224, direction = "sunrise", POSIXct.out = TRUE)
##        day_frac                time
## newlon   0.3924 2004-12-24 09:25:05
## Sunset
sunriset(Hels, d041224, direction = "sunset", POSIXct.out = TRUE)
##        day_frac                time
## newlon   0.6352 2004-12-24 15:14:42
######################################################### 

Here is the code for Culebra

culebra = matrix(c(-65.252, 18.297), nrow = 1)
culebra
##        [,1] [,2]
## [1,] -65.25 18.3
Culebra <- SpatialPoints(culebra, proj4string = CRS("+proj=longlat +datum=WGS84"))
Culebra
## SpatialPoints:
##      coords.x1 coords.x2
## [1,]    -65.25      18.3
## Coordinate Reference System (CRS) arguments: +proj=longlat
## +datum=WGS84

Calculate a Dawn and Dusk time

Nautical dawn

dawn = crepuscule(Culebra, b30414_Seq1$AST_ct, solarDep = 12, direction = "dawn", 
    POSIXct.out = TRUE)
str(dawn)
## 'data.frame':    99 obs. of  2 variables:
##  $ day_frac: num  0.21 0.21 0.21 0.21 0.21 ...
##  $ time    : POSIXct, format: "2012-07-19 05:03:06" "2012-07-19 05:03:06" ...
b30414_Seq1$Dawn = dawn$time

Nautical dusk

dusk = crepuscule(Culebra, b30414_Seq1$AST_ct, solarDep = 12, direction = "dusk", 
    POSIXct.out = TRUE)
b30414_Seq1$Dusk = dusk$time

Calculate sunrise and sunset time

Sunrise

sunrise = sunriset(Culebra, b30414_Seq1$AST_ct, direction = "sunrise", POSIXct.out = TRUE)
b30414_Seq1$Sunrise = sunrise$time

Sunset

sunset = sunriset(Culebra, b30414_Seq1$AST_ct, direction = "sunset", POSIXct.out = TRUE)
b30414_Seq1$Sunset = sunset$time
str(b30414_Seq1)
## 'data.frame':    99 obs. of  25 variables:
##  $ OBJECTID   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ TRANSMITTER: int  30414 30414 30414 30414 30414 30414 30414 30414 30414 30414 ...
##  $ DETECTEDID : Factor w/ 1 level "A69-1601-30414": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DATETIME_  : Factor w/ 99 levels "7/20/2012 2:11:46",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ X          : num  986 955 973 968 958 ...
##  $ Y          : num  1147 1010 996 981 971 ...
##  $ D          : num  1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 ...
##  $ LAT        : num  18.3 18.3 18.3 18.3 18.3 ...
##  $ LON        : num  -65.3 -65.3 -65.3 -65.3 -65.3 ...
##  $ n          : int  3 4 2 7 3 1 10 1 10 1 ...
##  $ HPE        : num  22.2 7.4 7.1 7.5 9.4 7.7 7.3 9.8 18.5 18.7 ...
##  $ HPEm       : logi  NA NA NA NA NA NA ...
##  $ TEMP       : logi  NA NA NA NA NA NA ...
##  $ DEPTH      : logi  NA NA NA NA NA NA ...
##  $ ACCEL      : logi  NA NA NA NA NA NA ...
##  $ DRX        : Factor w/ 82 levels "VPS02 VPS05 VPS06 VPS11",..: 4 62 34 21 26 75 63 82 5 82 ...
##  $ URX        : Factor w/ 82 levels "VPS02 VPS05 VPS06 VPS11 VPS12",..: 3 62 35 22 27 75 63 82 4 82 ...
##  $ UTCTime_ct : POSIXct, format: "2012-07-20 02:11:46" "2012-07-20 02:17:54" ...
##  $ AST_chr    : chr  "2012-07-19 22:11:46" "2012-07-19 22:17:54" "2012-07-19 22:20:26" "2012-07-19 22:21:27" ...
##  $ AST_ct     : POSIXct, format: "2012-07-19 22:11:46" "2012-07-19 22:17:54" ...
##  $ AST_lt     : POSIXlt, format: "2012-07-19 22:11:46" "2012-07-19 22:17:54" ...
##  $ Dawn       : POSIXct, format: "2012-07-19 05:03:06" "2012-07-19 05:03:06" ...
##  $ Dusk       : POSIXct, format: "2012-07-19 19:51:23" "2012-07-19 19:51:23" ...
##  $ Sunrise    : POSIXct, format: "2012-07-19 05:54:47" "2012-07-19 05:54:47" ...
##  $ Sunset     : POSIXct, format: "2012-07-19 18:59:45" "2012-07-19 18:59:45" ...

Calculate Light (can be “Dawn”, “Day”, “Dusk”, or “Night”)

b30414_Seq1$Light = ifelse(b30414_Seq1$AST_ct < b30414_Seq1$Dawn, "Night", ifelse(b30414_Seq1$AST_ct < 
    b30414_Seq1$Sunrise, "Dawn", ifelse(b30414_Seq1$AST_ct < b30414_Seq1$Sunset, 
    "Day", ifelse(b30414_Seq1$AST_ct < b30414_Seq1$Dusk, "Dusk", "Night"))))
b30414_Seq1$Light = factor(b30414_Seq1$Light)
str(b30414_Seq1)
## 'data.frame':    99 obs. of  26 variables:
##  $ OBJECTID   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ TRANSMITTER: int  30414 30414 30414 30414 30414 30414 30414 30414 30414 30414 ...
##  $ DETECTEDID : Factor w/ 1 level "A69-1601-30414": 1 1 1 1 1 1 1 1 1 1 ...
##  $ DATETIME_  : Factor w/ 99 levels "7/20/2012 2:11:46",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ X          : num  986 955 973 968 958 ...
##  $ Y          : num  1147 1010 996 981 971 ...
##  $ D          : num  1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 ...
##  $ LAT        : num  18.3 18.3 18.3 18.3 18.3 ...
##  $ LON        : num  -65.3 -65.3 -65.3 -65.3 -65.3 ...
##  $ n          : int  3 4 2 7 3 1 10 1 10 1 ...
##  $ HPE        : num  22.2 7.4 7.1 7.5 9.4 7.7 7.3 9.8 18.5 18.7 ...
##  $ HPEm       : logi  NA NA NA NA NA NA ...
##  $ TEMP       : logi  NA NA NA NA NA NA ...
##  $ DEPTH      : logi  NA NA NA NA NA NA ...
##  $ ACCEL      : logi  NA NA NA NA NA NA ...
##  $ DRX        : Factor w/ 82 levels "VPS02 VPS05 VPS06 VPS11",..: 4 62 34 21 26 75 63 82 5 82 ...
##  $ URX        : Factor w/ 82 levels "VPS02 VPS05 VPS06 VPS11 VPS12",..: 3 62 35 22 27 75 63 82 4 82 ...
##  $ UTCTime_ct : POSIXct, format: "2012-07-20 02:11:46" "2012-07-20 02:17:54" ...
##  $ AST_chr    : chr  "2012-07-19 22:11:46" "2012-07-19 22:17:54" "2012-07-19 22:20:26" "2012-07-19 22:21:27" ...
##  $ AST_ct     : POSIXct, format: "2012-07-19 22:11:46" "2012-07-19 22:17:54" ...
##  $ AST_lt     : POSIXlt, format: "2012-07-19 22:11:46" "2012-07-19 22:17:54" ...
##  $ Dawn       : POSIXct, format: "2012-07-19 05:03:06" "2012-07-19 05:03:06" ...
##  $ Dusk       : POSIXct, format: "2012-07-19 19:51:23" "2012-07-19 19:51:23" ...
##  $ Sunrise    : POSIXct, format: "2012-07-19 05:54:47" "2012-07-19 05:54:47" ...
##  $ Sunset     : POSIXct, format: "2012-07-19 18:59:45" "2012-07-19 18:59:45" ...
##  $ Light      : Factor w/ 1 level "Night": 1 1 1 1 1 1 1 1 1 1 ...

It appears that this sequence happened entirely at night!