## What part of the sky spends the most time in darkness?

*SUMMARY: I set out to answer the question: considering the perspective of everyone on Earth and accounting for the seasonal changes in daylight, what spots in the night sky are least and most visible? The incremental steps to get there illustrate the effect of latitude and the seasons on the night sky. The final map is shown above, but for those interested, a full explanation with code snippets is shown.*

When I moved to the London area from the US last year, I noticed that not only were the summer nights shorter than in New Jersey, but twilight seemed to go on forever. London is at about 52° N, so for a few days around the summer solstice, the sky never reaches astronomical twilight (when the sun is 15° below the horizon, which is also when it is dark enough for imaging).

So I wondered how much imaging time I get each year for each area of the sky. (Ignoring weather of course.) And beyond my view of the sky, what does the picture looks like when you consider the view from everywhere on the planet. Obviously, the summer sky isn’t in darkness as long as the winter sky, but how big is the difference? And to take the question to a larger level, what happens when you consider where people live on the earth?

I’m sure there is a way to solve this analytically for someone really talented in spherical trigonometry, but I set out to solve it numerically. I had a lot of R code that I’d written to create the maps in *The Astrophotography Sky Atlas*, so how hard could it be? It turns out that it was a good bit of effort for such an unremarkable question. (At least 30 hours, including this write up.)

Let’s start with my assumptions. A point in the sky is considered “image-able” when:

- It’s fully dark: the sun is at least 15° below the horizon
- The point is at least 20° above the horizon. Lower than that, and you’re looking through over three airmasses.

First I had to write and test the code in R to do all the positional astronomy calculations. This turned out to non-trivial, but after a couple of evenings reading and testing, I had functions that returned comparable values to the online ephemerides. The sunrise and sunset values I get differ by a few minutes from various online sources, so there are probably more precise equations available, but let’s call these good enough. I’ve posted snippets of them here for critique and transparency. These are all of the helper functions I wrote.

SunCoords<-read.csv("SunRADec2017.csv") #get sun coordinates for any year, expressed in decimal degrees CalcSunRiseSetTime <- function(lat, day, altitude, sunset=TRUE){ #calculates time sun is at altitude. Default is sunset=descending sun. d<-day+6208.5 #6208.5 days from Jan 1 2000 noon to Jan 1 2017 0:00, which is the year I'm using here ST0<-GMT2GMSThours(d)*15 #calculate sidereal time at long=0 at midnight, expressed in degrees #use ifelse so function is vectorized... #modulus handles discontinuity near equinox EoT<-((SunCoords[day,]$RADegrees - ST0 - 0) / 15)%%24 #"Equation of time" at long=0, expressed in hours #semi-diurnal arc cosH<-(sin(altitude/57.2957795) - sin(lat/57.2957795)*sin(SunCoords[day,]$DecDegrees/57.2957795)) / (cos(lat/57.2957795) * cos(SunCoords[day,]$DecDegrees/57.2957795)) #handle cases near poles: if cosH >1, return ~11:59pm, as sun never rises. if <-1, return ~12:01am, as sun never sets if(sunset){ r<-ifelse(cosH > 1, 0.001, ifelse(cosH < -1, 23.999, EoT + (acos(cosH)*57.2957795/15))) } else{ r<-ifelse(cosH > 1, 23.999, ifelse(cosH < -1, 0.001, EoT - (acos(cosH)*57.2957795/15))) } r #time in decimal hours } CalcObjectRiseSetTime <- function(lat, day, RA, Dec, altitude, set=TRUE){ #RA and Dec of object in degrees cosLHA<-(sin(altitude/57.2957795) - sin(lat/57.2957795)*sin(Dec/57.2957795)) / (cos(lat/57.2957795) * cos(Dec/57.2957795)) if(set){ if(cosLHA > 1) return(0.001) if(cosLHA < -1) return(23.999) SiderealAnswer<-RA + (acos(cosLHA)*57.2957795) #in degrees } else{ if(cosLHA > 1) return(23.999) if(cosLHA < -1) return(0.001) SiderealAnswer<-RA - (acos(cosLHA)*57.2957795) #in degrees } #convert Greenwich sidereal time to GMT GMST2GMT(day, SiderealAnswer/360*24) } GMST2GMT<- function(day, ST){#input day of 2017 and sidereal time in hours ***why is it up to 8 minutes off? #return value in hours ((ST/24-0.277-day*0.00273042)*24)%%24 #2017 starts at 0.277 day sidereal time, declines by 3m56s each day } GMT2GMSThours<- function(GMTday){ #input days and fraction since J2000.0 (18.697374558 + 24.06570982441908 *GMTday)%%24 #return sidereal time at long=0, expressed in hours } TimeInNightSky <- function(lat, day, RA, Dec, MinObjectAltitude, SunAltitude){ # ObjectRise<-CalcObjectRiseSetTime(lat, day, RA, Dec, MinObjectAltitude, FALSE) ObjectSet<-CalcObjectRiseSetTime(lat, day, RA, Dec, MinObjectAltitude, TRUE) SunRise<-CalcSunRiseSetTime(lat, day, SunAltitude, sunset=FALSE) SunSet<-CalcSunRiseSetTime(lat, day, SunAltitude, sunset=TRUE) if(ObjectRise > ObjectSet){ time<-min(ObjectSet, SunRise) + 24 - max(ObjectRise, SunSet) } else{ time<-max(0, SunRise-ObjectRise) + max(0, ObjectSet-SunSet) } time } CalcEoT <- function(day){ d<-day+6208.5 #6208.5 days from Jan 1 2000 noon to Jan 1 2017 0:00, which is the year I'm using here ST0<-GMT2GMSThours(d)*15 #calculate sidereal time at long=0 at midnight, expressed in degrees #use ifelse so function is vectorized... EoT<-((SunCoords[day,]$RADegrees - ST0 - 0) / 15)%%24 #"Equation of time" at long=0, expressed in hours #EoT<-ifelse(abs(EoT) > 12, EoT+24, EoT) #handle discontinuity near vernal equinox EoT } #calculate map data for one night at one location lati<- 25 #latitude to use minAlt<-20 #lowest imageable altitude lowDec<-max(-89, lati-90+minAlt) #only compute values for Dec within range of imageability for latitude highDec<-min(89, 90+lati-minAlt) decSequence<-seq(lowDec,highDec,2) RASequence<-seq(1,365,2) timegrid=data.frame(matrix(0, nrow=180, ncol=365)) for (r in RASequence){ for(d in decSequence){ timegrid[d+90, r]<-timegrid[d+90, r] + TimeInNightSky(lati, 105, r, d, minAlt, -15) #-15 is position of sun at astronomical twilight } } #now convert the dataframe above into a format readable by the plotting code timemap2 <- data.frame(RA=numeric(), Dec=numeric(), time=numeric() ) for (x in RASequence){ #only look where values were loaded into the matrix .. old version was "1:ncol(timegrid)){" for(y in decSequence+90){ #1:nrow(timegrid)){ temp<-data.frame(x, y-90, timegrid[y,x]) #the +/-90 is to handle the negative Dec values, which are invalid indexes in a matrix cols <- c("RA","Dec","time") colnames(temp) <- cols timemap2<-rbind(timemap2, temp) } } #flip the RA axis timemap2$RA<-365-timemap2$RA

With these functions, R can easily plot a few interesting things to verify the code makes sense. One is the Equation of Time.

Sunrise and sunsets (astronomical twilight times, actually) for a year look pretty good, too. This is sunrise and sunset for every day of the year at 25° N, not accounting for Daylight Savings Time of course.

The functions appear to work well, so I iterated those functions to calculate a single night sky (April 15^{th}) from London, to see what parts get the most image-able hours.

#calculate map data for one night at one location lati<- 25 #latitude to use minAlt<-20 #lowest imageable altitude lowDec<-max(-89, lati-90+minAlt) #only compute values for Dec within range of imageability for latitude highDec<-min(89, 90+lati-minAlt) decSequence<-seq(lowDec,highDec,2) RASequence<-seq(1,365,2) timegrid=data.frame(matrix(0, nrow=180, ncol=365)) for (r in RASequence){ for(d in decSequence){ timegrid[d+90, r]<-timegrid[d+90, r] + TimeInNightSky(lati, 105, r, d, minAlt, -15) #-15 is position of sun at astronomical twilight } } #now convert the dataframe above into a format readable by the plotting code timemap2 <- data.frame(RA=numeric(), Dec=numeric(), time=numeric() ) for (x in RASequence){ #only look where values were loaded into the matrix .. old version was "1:ncol(timegrid)){" for(y in decSequence+90){ #1:nrow(timegrid)){ temp<-data.frame(x, y-90, timegrid[y,x]) #the +/-90 is to handle the negative Dec values, which are invalid indexes in a matrix cols <- c("RA","Dec","time") colnames(temp) <- cols timemap2<-rbind(timemap2, temp) } } #flip the RA axis timemap2$RA<-365-timemap2$RA

Then I plotted the resulting data (using code not shown):

For comparison, I also plotted London’s sky at the winter and summer solstices:

You can see that in the winter, London provides plenty of dark sky time (though it’s rarely clear), but near the summer solstice the sun never gets 15° below the horizon!

For comparison, let’s calculate the same chart for Miami (at 25° N) on April 15^{th} to see how things change. (Because the Earth is rotating, longitude is irrelevant to our calculations.)

You can see how the longer spring nights in Miami, compared to London, give you more dark sky time this time of year. The opposite is true for winter, where the nights are looooong in London.

Now that I had one night, it was time to sum all of the nights for a year. Even on a fairly fast computer, this took some time, so I had to go back and see if I could make the code more efficient. R is an interpreted language, but highly optimized for vectors. Unfortunately, this code required several “for” loops, which R does not handle quickly. There is probably an elegant and fast solution with dplyr or something, but I couldn’t make it work. Using a resolution of 2×2° squares, it still took 40-80 minutes to crunch the numbers for each map, depending on latitude.

Here is the map from 50° N.

How do the maps look near the equator? And what about near the poles? Here are maps at 70°N, 20°N, and the equator.

You can see that the seasonal changes are greater the further you move away from the equator.

Now that everything was working, I generated full-sky data for latitudes of 80° N down to 80° S in 2.5° increments. This took about 60 hours of computer time (mostly while I was sleeping or at work). I’m sure there is a more elegant and faster way to do it in R, but again, I used for loops. This gave me the data I needed to calculate the “average” view from everywhere on Earth.

###avoid rbind by using a matrix, then re-shape the matrix back into a dataframe the plotting code can read for(lati in seq(-80, 80, 2.5)) { minAlt<-20 #lowest imageable altitude lowDec<-max(-89, lati-90+minAlt) #only compute values for Dec within range of imageability for latitude highDec<-min(89, 90+lati-minAlt) decSequence<-seq(lowDec,highDec,2) RASequence<-seq(1,365,2) timegrid=data.frame(matrix(0, nrow=180, ncol=365)) for (r in RASequence){ for(d in decSequence){ for(date in 1:365){ timegrid[d+90, r]<-timegrid[d+90, r] + TimeInNightSky(lati, date, r, d, minAlt, -15) #-15 is position of sun at astronomical twilight } } } #get the average timegrid<-timegrid/365 #covert the grid into a dataframe of the right format to be read by mapping code below timemap2 <- data.frame(RA=numeric(), Dec=numeric(), time=numeric() ) for (x in RASequence){ #only look where values were loaded into the matrix .. old version was "1:ncol(timegrid)){" for(y in decSequence+90){ #1:nrow(timegrid)){ temp<-data.frame(x, y-90, timegrid[y,x]) #the +/-90 is to handle the negative Dec values, which are invalid indexes in a matrix cols <- c("RA","Dec","time") colnames(temp) <- cols timemap2<-rbind(timemap2, temp) } } #flip the RA axis, since map is printed from 24 to 0 hrs timemap2$RA<-365-timemap2$RA #save as a separate file fname=paste("NightMap Lat=", lati, ".Rda", sep="") save(timemap2, file=fname) }

Bringing all of my maps together, this is the view from everywhere on earth (technically, everywhere between latitudes 80° N and 80° S):

Two things stand out:

- The contours are a little jagged. This is a result of running the calculations for latitudes at 2.5° increments. Increasing the “latitude resolution” would further smooth them.
- The hemispheres aren’t completely symmetrical. Summer in the southern hemisphere seems to have a slightly greater impact on the southern night sky than the northern hemisphere’s summer does on the northern sky. At first I thought it could be an error in my equations (and it could indeed be). But I think it’s the effect of Earth’s elliptical orbit. The changes in sunlight on the days near the solstices are not symmetrical (see the shape of the analemma for proof). Someone more skilled in positional astronomy could surely offer a clearer explanation or show me I’m wrong.

All in all, I think it’s a pretty neat map. The winter sky clearly gets more time in darkness, and the southern hemisphere has more extreme differences between the summer and winter. Surprisingly, the area of sky with the least good imaging time on average is just south of the most-imaged object in the night sky: M42 in Orion. And the area with most imaging time available on average is the stinger of Scorpius, where there are a slew of beautiful objects.

But we can take this even further. You may object that people are not evenly distributed by latitude. And indeed, that is why you see so many images of M42 compared to the Cat’s Paw Nebula, which gets more hours in darkness. A proper map of “image-able” time for the night sky would incorporate where people actually live by latitude. So let’s account for where people live. First, I had to get the data, which had to be distilled from a much larger data set.

The data on population from SEDAC* come in a raster format, so it took a little more code to get a simple percentage of population at each latitude. (Note that in the code here, I cut the data slightly differently, so the bins would be centered on 5° increments. )

library(rgdal) library(sp) library(raster) GPW<-raster("gpw-v4-population-count-adjusted-to-2015-unwpp-country-totals_2015.tif") popPerLatDegree<-aggregate(GPW, c(43200,1200), fun=sum) #aggregate data across longitudes and into 2.5 degree groups by latitude popPerLatDegree<-as.data.frame(popPerLatDegree, xy=TRUE) names(popPerLatDegree)[3] <- "pop" #shorten column name popPerLatDegree$x<-NULL #we don't need the longitude, since there's only one value now popPerLatDegree[is.na(popPerLatDegree$pop),]$pop<-0 #replace the NA's with zeros popPerLatDegreeTrimmed<-popPerLatDegree[3:142,] #trim top two degrees off so center of each group of five degrees (approximately) aligns with 0,5,10,etc. #sum every five degrees into one value popPerFiveDegree <- data.frame(centralLat=numeric(), pop=numeric() ) for(n in seq(1, nrow(popPerLatDegreeTrimmed), 5)) { total<-0 for(y in 0:4) { total<-total + popPerLatDegreeTrimmed[n+y,]$pop } tmp<-data.frame(popPerLatDegreeTrimmed[n+2,]$y, total) cols <- c("centralLat","pop") colnames(tmp) <- cols popPerFiveDegree<-rbind(popPerFiveDegree, tmp) } #calculate proportion in each band worldPop<-sum(popPerFiveDegree$pop) popPerFiveDegree$pop<-popPerFiveDegree$pop/worldPop

The following table shows the percent of the world’s population by latitude:

A couple of things are immediately apparent:

- Only 13% of the world’s population lives in the southern hemisphere.
- 49% of the world’s people live in the range from 20 to 40° N.

Okay, now that we know where people live, we can re-create the sky map weighted by where people live.

Sadly, it turned out to be pretty uninteresting. The results are simpler than I first thought because the distribution of human population is so skewed that its effect swamps any seasonality. The final plot mostly consists of data from the mid-northern latitudes, with only a tiny contribution from other latitudes. The further south you go in declination, the fewer people there are to see it. The short answer to my original question is that the most and least image-able parts of the sky is basically, “near the north and south celestial poles, respectively.”

* Population data from: Center for International Earth Science Information Network – CIESIN – Columbia University. 2016. Gridded Population of the World, Version 4 (GPWv4): Population Count Adjusted to Match 2015 Revision of UN WPP Country Totals. Palisades, NY: NASA Socioeconomic Data and Applications Center (SEDAC). http://dx.doi.org/10.7927/H4SF2T42. Accessed 20 April 2017.