This map shows the brightness1 of the night sky from Maryland to Maine, constructed using Leaflet, mostly via the {leaflet} package in R and its extension packages. Click anywhere on the map, and you'll get (by default) the closest single point that is at least 1 Exposure Value (EV) darker than the clicked point.2 You can use the "dark point properties" control to change: the number of dark points, their EV difference from the clicked point, and their maximum distance away.
My inspiration was the awesome lightpollutionmap.info,3 my deep dislike of rainbow color palettes, and my desire to take some good star photos in the Catskills. Data were downloaded from Supplement to: The New World Atlas of Artificial Night Sky Brightness.4
The map is available at: https://cgettings.github.io/Light-Pollution-Map/
The sky brightness data comes from a whole-world geotiff file of simulated zenith radiance data, given in mcd/m^2. I read the file into R as a RasterLayer using raster::raster(), which avoids loading the whole 2.8 GB file into memory. I then crop it using a bounding box that contains the states I'm mapping, and then convert the cropped raster into a stars object, which loads it into memory. This allows me to use {sf} methods to further crop the raster using state boundary data from the {tigris} package. Finally, I convert the stars object into a tbl, convert luminance (mcd/m^2) to magnitude (mag/arcsec^2),5 and then turn the tibble back into a stars object for mapping.
I create the map using {leaflet} and addProviderTiles(), with a custom tile layer drawn from the "USA_Topo_Maps" esri tile layer hosted on ArcGIS online. I then add the stars raster to the map using leafem::addGeoRaster(), specifying colors using the viridisLite::inferno() palette, with breaks computed so that qualitatively different colors roughly map onto the breaks between Crumey's sky quality bands: respectively, Pristine vs. Black (purple), Black vs. Grey (maroon), Grey vs. Bright (orange), and Bright vs. White (yellow).
I add the brightness mousover display using leafem::addImageQuery(), and the map view reset botton by modifying leaflet.extras::addResetMapButton() to simply add a position argument to the easyButton() call. Finally, I add map dependencies using the registerPlugin and leaflet.extras::addAwesomeMarkersDependencies() functions, and then pass my custom JavaScript and a tbl of raw raster data to htmlwidgets::onRender().
This custom JavaScript uses the raster data I passed to htmlwidgets::onRender(), along with the already-rendered map layer, to find the dark points.
First, the script converts the raw raster data to LatLng points, with the magnitude value saved as the alt[itude] property. It then re-reads the raster layer data from the DOM using fetch, then uses the georaster package (already loaded thanks to leafem::addGeoRaster()) to parse the data.6 The script then uses the geoblaze package to extract the magnitude value from the selected point. This point can be selected by clicking on the map, by searching for a location using the search bar (which uses Control.Geocoder), or by a user finding their location using the "locate" button (which uses Leaflet's map.locate method, implemented via an EasyButton).
Using the raster's value at the selected point, the script by default finds all points in the raw raster data that are between 1 and 27 EVs darker than the selected point. The script then uses Leaflet's built-in distanceTo function to compute the distance between the selected point and the filtered dark points, and finally selects the closest single point.
The final dark points are then displayed on the map, with tooltips giving their brightness, distance, and coordinates; the unformatted values are also sent to the console (accessible via "Developer tools" in a web browser) as LatLng objects. When using the location search, the selected location's properties are sent to the console as well.
Using the dark point control, the user can change the number of dark points returned, the EV difference between the selected and dark points (including negative difference values, which represent brighter points), and the maximum distance between the selected and dark points. These user-specified control values are read from the DOM when the map is clicked, when the dark point control's "update" button is pressed, and when a location is found.
TODO: Add options to:
Show more than 1 dark pointShow darkest point(s) within a specified radius of clicked pointChange magnitude difference between clicked point and dark points- Find dark points where specified celestial objects are visible8
Add explanation re: heat map color gradientAdd more/better search functionality- Add dynamic URL hash
1. Sky brightness values are in mag/arcsec^2. Explanation here. ↩
2. In practical terms, this means that you could e.g. find locations where a 2-second camera exposure had the same background sky brightness as a 1-second exposure at your original location. ↩
3. Shout out to Dan Jentzen for introducing me to this website. ↩
4. Falchi, Fabio; Cinzano, Pierantonio; Duriscoe, Dan; Kyba, Christopher C. M.; Elvidge, Christopher D.; Baugh, Kimberly; Portnov, Boris; Rybnikova, Nataliya A.; Furgoni, Riccardo (2016): Supplement to: The New World Atlas of Artificial Night Sky Brightness. V. 1.1. GFZ Data Services. http://doi.org/10.5880/GFZ.1.4.2016.001
Falchi F, Cinzano P, Duriscoe D, Kyba CC, Elvidge CD, Baugh K, Portnov BA, Rybnikova NA, Furgoni R. The new world atlas of artificial night sky brightness. Science Advances. 2016 Jun 1;2(6). http://dx.doi.org/10.1126/sciadv.1600377 ↩
5. Because these luminance values represent only artificial brightness, I add 0.171168465 to each value before converting it, copying the procedure of lightpollutionmap.info. This ensures that a luminance of 0 produces a magnitude of 22.0, which is the darkest value on common light pollution scales (labeled "Excellent" or "Pristine") ↩
6. This is necessary because the raster data object created by leafem::addGeoRaster() only exists within the scope of the function call that adds the georaster layer. ↩
7. Technically, it's "EV diff value + 1 EV". This is a mostly arbitrary cutoff that reduces processing demands. Note that if a user is searching for brighter points, this range will include points 1-2 EVs brighter, because the script accounts for the sign of the EV control. ↩
8. More information available here. ↩