NASA Night Lights – Workflow

Projects & Blog
Auxiliary GIS



NASA Night Lights – Workflow

As mentioned in my previous post. I wanted to share the workflow I used to come to the end-state of a hosted web map. Anyone interested can likely use the same workflow to come to a similar outcome. This post is going to be quite a long one, but hopefully the detail will help anyone working on something similar.

For this project I used more software than I normally do, and at this time I haven’t had as much time to investigate alternatives.. which may or may not be better for this particular project. Below is what I used:

  • ArcMap (or possible open source QGIS alternative)
  • Photoshop (or possible open source GIMP or ImageMagick alternative)
  • GDAL (GDAL is open source, and I’m not aware of a better tool for the job)

Step 1: Data Download


The first step was to get the raw imagery from the NASA download page. This page has the highest resolution I could find and the important grayscale lights. On this page you’ll find five datasets. The processed color imagery with lights on a dark background for 2012 and 2016, intensity values (grayscale imagery of lights only) for 2012 and 2016, and the previous color composite for 2012 using their old technique. For this setup, I used every set but the old 2012 color composite technique. You can download the sets in full resolution (500m) by downloading the tiles for each set. I highly recommend starting with the GeoTIFFs as they will already be georeferenced, and all you have  to do is mosaic them. The colored images will provide the Basemaps for the online web map to switch between, and the intensity values will be used for actual processing and analysis.

In the thumbnails above, the bottom set (the intensities) may not look like much, but they are a critical piece to determining light change. I believe the grayscale section are the same data without background imagery. I didn’t see any provided measurements so this isn’t perfect for metrics, but essentially the higher the number the more intense light was for that pixel. Since most of the earth is not lit, it has a lowest value of 0 (or black) and the small white/gray dots you see are the few areas that are very intense and therefore closer to white.

Although a similar method of detecting change could be used with the original imagery, this set is significantly more accurate. There is no imagery in the background to “alter” the color of lights or their bloom, just raw values to work with. This means that if we were to take a particular pixel and see that it had a value of 0 for 2012, then see that same pixel had a value of 100 in 2016, it was significantly brighter. This applies vice versa and allows us to accurately determine which pixels were actually brighter, which were dimmer and which ones stayed exactly the same. My attempts at a change detection with imagery classification did yield similar results (which I unfortunately didn’t save for posting/comparison). However, the blooms seemed to deviate in some areas quite a bit. There were occasional random spots that would disappear based on my classification breaks. I just didn’t find it to be as reliable. How much red, green, and blue do you need to take away to consider light loss, especially when the background deviates in all three? I think it is best to just stick with the linear dataset.

Step 2: Data Organization


After downloading the four datasets (which ends up being several GBs), I’m ready to start some work. I end up simply dumping all of them into 4 separate Mosaic Datasets as the easiest way to process them. This leaves the original data untouched for easy use later on, and allows us to process them without creating extra data. It is perfectly fine to mosaic them into their own separate rasters if you prefer to do this, I would just recommend that you avoid applying any permanent stretch on the intensity rasters that might alter values. Here’s an example of my setup for the rasters:

Step 3: Calculating Light Intensity Values


Time to single out the intensity datasets and prep them for simple raster operations. However, they are provided in RGB format. Since the raster is grayscale, all bands (Red, Green, and Blue) have identical values (i.e., [0,0,0] for black, [255, 255, 255] for white, and [115,115,115] as a shade of gray). This means we can use any band we would like as a singular band to set as an intensity value. The quickest way I found to do this was to use the Imagery Analyst window in ArcMap. You can highlight your layer (such as BlackMarble2012) and single out one band by using a Band Arithmatic function. If you are new to functions, you can add them to any raster using the button highlighted below:

If you already have a mosaic dataset like I did, you can use this button or check out the properties of your raster and then the functions tab as an alternative. From this functions menu, right click on any of your steps and ‘Insert Function…’ > ‘Band Arithmetic’. Set the expression to your first, second, or third band. Since they all have the same values it doesn’t matter which you choose, as long as you only set one of them. In the below shots, I just went with band 1.

There are a few other ways that you could come up with similar results, such as the ‘Grayscale Function’ or Make Raster Layer tool. This was just my personal preference. The tool may be ideal for someone wanting to save an entirely new file, but I prefer to work with the singular mosaic dataset so I haven’t added any extra data to my workspace yet. Also functions allow you to go back and edit them if you didn’t get your parameters quite right. Totally awesome!

One other thing we’ll need to consider is the bitage of our output. Initially ArcMap is often going to want to turn our raster into a 16bit or 32bit raster. Which is unnecessary and likely to cause problems later. If you chose to utilize the function approach you can go to the second tab on any function and change the bit depth dynamically. Otherwise, when you save your raster, ensure you have a bit depth of 8-bit unsigned. I always like to think of it as unsigned has NO sign (to include a negative sign) so therefore you have no negative numbers and must be positive only. Now the values will accurately represent the intensity values we want to work with.

Once we have this setup, our raster should be a single band 8bit unsigned raster. This allows us to apply a color ramp and do some basic math to determine intensity changes. This is as simple as using the Minus Tool. With this tool a pixel value from the first raster will have the spatially relative pixel value from the next raster subtracted and the value after this operation will be applied to a new raster. In this event, if we take the 2012 intensity values and subtract the 2016 values we should get a raster of positive and negative numbers. If the numbers are negative, then the value in 2012 was lower in 2012 than 2016 since less intensity in 2012 is being subtracted from the higher intensity in 2016. Negative numbers mean the light in the area grew. If the numbers were the same, such as where there were no lights for both years, or the intensity value stayed the same, we would get 0 to indicate no change. If the value was positive the opposite would be true, the value during 2012 was higher and therefore still positive after being subtracted by a lower 2016 intensity number. This raster gives us all the information we need to quantify whether light grew or receded and by how much. All negative numbers indicate growth, positive numbers indicate receding light, zeroes indicate stagnation, and the further the value is from zero sets the intensity of that change.

Step 4: Applying Color Ramps to Intensity Values


Now to apply color to the setup, and this is where I ran into some major issues. I haven’t had the time to fully investigate QGIS’ capabilities, but I know it works with transparent values significantly better than ArcMap, but I’m not sure about its export capabilities. ArcMap (and even ArcGIS Pro) have a major issue where they simply cannot export transparency from a raster. There is no transparent options to apply in ArcMap (setting a null color on a ramp or choosing ‘No Color’ applies a gray/white in data view). In ArcGIS Pro you can setup transparency with the color ramp and it will display fine in-program, but exporting the layer always results in an off-white color instead of transparency. Kinda crappy, honestly.

At any rate the first step is to actually apply the color. I started by copying and pasting a second version of my difference layer. One for growing light and one for receding light. Since the negative numbers indicate growth, I’ll start with them. I started by applying a stretch to my raster and setting the max value to 0 and the minimum value to -255. There are a few other considerations when working here, such as choosing the highest intensity value (whether positive or negative) and making that the negative/positive limits to equalize the intensity ranges, but this will be taken care of in a later step, and I believe there are some maxed out values in this dataset anyways so it may be unnecessary in this event. Just something to consider.

You’ll also notice I applied a green to white gradient. Obviously any desired color could be chosen, but if working in ArcMap you want to choose white for the lower intensity values and not ‘No Color’. No color, looks fine in the ramp diagram, but it will not be applied to the actual raster and may not have a perfect white transfer (which is important for transforming the image later). This version will now reflect the growing light areas based off of intensity. I would include a screenshot of this step, but it is honestly unimpressive against a white background unless zoomed in to see certain areas in full detail. We can rinse and repeat for the second raster layer (or receding areas) by only altering the color and ramp direction.

Now our two rasters are all set. These can be saved to .PNG files as single file outputs to work with later once transparency is added.

Step 5: Applying Transparency to Intensity Values


We almost have a complete set with the Black Marble imagery sets and two difference layers. Of course, we now have to go through the tediousness that is adding transparency to our layers first. Again QGIS may have a much better solution for exporting transparencies to allow you to skip these next few steps. Without such an export we can still add transparency through programs such as GIMP or Photoshop. I, personally, am far more comfortable with Photoshop, so I will be using that in the next examples. Similar methods should be achievable with GIMP or ImageMagick(command-line) as well.

Load any of the difference .PNG rasters into Photoshop, it should look as expected a few red or green pixels in a sea of white. Our goal is to change any level of whiteness to transparency in a 1 for 1 match (i.e., a pure white pixel will be completely transparent, a slightly red/green pixel will be partially transparent, and a pixel that is fully green (0,255,0) or red (255,0,0) that has no “white” will not be altered and have no transparency. I was not able to find a simple tool/command to do this, but with a couple of tricky layer alterations it is actually easy to do.

When we first load the .PNG version of our difference raster, it should look something like this:

The first step I took to change this layer to transparent was to add an [Menu]Layer > Adjustment Layer > Hue/Saturation. Setting the saturation bar all the way to the left will turn our raster into grayscale. This simplifies the next few steps to define a transparent gradient.

From that same menu we add the [Menu]Layer > Adjustment Layer > Invert… which will reverse the color order, so black and white will switch. Ideally, you’ll have a mostly black raster with some gray mixed in for the intensity values. Here I apply a curve ([Menu]Layer > New Adjustment Layer > Curves…) in order to highlight the heaviest intensity values to “white” or non-transparent and scale them appropriately. This is why your choice based on Min-Max stretches earlier may have been unnecessary. If you did end up correctly calculating them then altering the curve may be unnecessary. You can also alter the curve to apply a desired amount of transparency. Here is how my finalized adjusted layer looks:

This gradient allows us to apply it to the original layer and utilize it as a transparent gradient. The black pixels will be 100% transparent and white pixels will be 0% transparent, the shade of gray will determine how transparent the original pixel value will be. Awesome!

First highlight the “Curves” layer and select all pixels: [Menu]Select > All (or ctrl+A). Then use [Menu]Edit > Copy Merged (or shift+ctrl+C). The regular copy command will not work, ensure you are copying the merged result. In my experience, this sometimes took a while, but it just may have been my poor laptop overworking itself. To apply this to your original layer, highlight the layer and add a layer mask via [Menu]Layer > Layer Mask > Reveal All. You’ll notice a new frame to the right of your original layer that indicates the mask. To edit this mask alt+click on it. Your background should go white showing you the raw mask. Here is where you can paste your previously copied merged layer to apply it.

After pasting the layer, turn off all of your other adjustment layers and highlight your original layer again to check out the result. You should see a checkered pattern where the black pixels used to be to indicate the pixels are transparent. You should be all set, save your .PNG as it now has appropriate transparency! Here is what my final output looked like:

You should be able to repeat this process with only slight modifications to the Curve Adjustment Layer (since the values are inverted) for your alternate difference raster.

Step 6: Re-georeferencing Intensity Values


Now here’s the bad part. We took one step forward by getting appropriate transparency, but we took one step back because Photoshop strips the georeferencing information from our .PNG. Doh! If there is a way/plugin to avoid this, I have not been aware at the time I created this. However, it is fairly easy to re-georeference this image. The best way I found is through GDAL. This is because any georeferencing in ArcMap requires saving, and I don’t think ArcMap saves with transparency… betcha saw that one coming. At any rate, I used the command:

in order to save to a new .PNG file and define the coordinate system and extent. Since we already know that this covers the entire world, we can set the extent to the extreme values for the latitudes (90, -90) and longitudes (-180, 180). The -of PNG options tells GDAL to save as a .PNG. The -a_ullr expects four values with the upper-left and lower-right coordinates in the format: -a_ullr <left_longitude> <top_latitude> <right_longitude> <bottom_latitude>. Seconds later, and you’ll have a clean version with transparency and spatial information.

The only change that is required for your alternative raster is updating the input and outputs. Afterwards you’ll be left with two .PNG files ready for web processing.

Step 7: Preparing Data for Web Mapping


At this step, there are numerous ways that you can go from here. My typical goal for such projects is to have a version of a web map that can

  • Be used in an “offline” environment as a simple structure.
  • Can be hosted without the need for server architecture setups (i.e., ArcServer, GeoServer).

You can use open source libraries and specifications to turn this data into an offline web map you can utilize from your computer’s local files. You can even upload these files to any raw storage you would like and they will continue to work with minimal (if any) changes. You can get a web map up and running using web technologies and libraries easily, you only need to worry about cost once you decide to host, which is typically not expensive either. I believe this project would cost me $0.20 a month if it ever got used heavily.

In order to use this data, I utilize the TMS tile specification which should work flawlessly for web libraries such as OpenLayers, Leaflet, etc. This is a genius way that data can be programmatically cut up into MANY very small tiles. They are ordered in folders with numbers and given numbered names. They will be non-georeferenced images that mapping libraries will know how to appropriately load based on the numbering system they are saved in. Essentially, the mapping libraries will load only a few tiles at any given time and since most are small (think 5-20kb). Even on very large monitors, you may only be loading 100 tiles at any given time. This will be true at any scale, so load/render times are consistent and you can load as much detail as you want to process/host. There are many other alternatives that may be more efficient and/or powerful, but this setup is simple to create, understand, and work with.

There are several simple ways to create tiles. If you are brand new and not a big fan of command-line, I would recommend checking out QTiles. It is an amazing plugin that will take your data as rendered in QGIS and create tiles out of the project. The user interface is fairly easy to understand, and the progress bar is awesome. This plugin does falter when your data is so large that QGIS has difficulty rendering it (shouldn’t be an issue with the relatively small size of the data in this project), and that it only uses one core (although you can process several layers from separate QGIS instances). However, I decided to go with GDAL for this one again. GDAL has gdal2tiles.py which is a simple command-line based option that will process your data quickly with quite a few nice options. It is actually the foundation of the QTiles plugin.

gdal2tiles.py uses a format similar to this: gdal2tiles.py -r near -z 3-8 <source file> <destination folder>. With this command, gdal2tiles takes your input file, and does all the work to cut up the tiles and save them in the supplied destination folder. So really all we need is something like our intensity .PNG files, and an empty folder to dump the tiles into. A command like this should start creating tiles for us:

Note the quotes, these are only required if your path has spaces in it, if you don’t have spaces they won’t hurt though. In this example GDAL will take our RecedingLights.png and save the TMS structure in a folder right next to it with the name RecedingLights_Tiles. -z 3-8 refers to the zoom levels. 0-2 were left out because I don’t expect users to zoom out that far, and since it provides no useful resolution it just cleans things up a bit. The resolution of this data (at 500m) is fully detailed by level 8. Anything further and I’m wasting tons of time processing and storing the extra data for no additional resolution. In fact, I could likely get away with 7 and interpolating level 8 if I wanted to save some space. You can build any levels you might like, but keep in mind that every additional level approximately quadruples in the number of tiles from the previous level. For the entire world, that adds up quick. Either way, here’s an example of how my output turned out.

Rinse and repeat for your alternate intensity raster, then on to the Black Marble Imagery!
Initially, if you were using Mosaic Datasets such as in this workflow, you’ll need to take a quick additional step as they are not supported by GDAL. If you happened to export them to a singular .TIF file (as an example), then you can skip this step.

GDAL can mosaic rasters on the fly similar to a mosaic dataset using virtual rasters (.VRTs). These don’t have all the power that a mosaic dataset has, but will quickly allow all of our .TIF files to act as one file which is perfect for use in things like QGIS or GDAL.

You can create a .VRT in seconds with a GDAL command like this:

Note the asterisk. Normally GDAL requires a list of files that you want to add to the .VRT, but I chose the wildcard (*) to tell GDAL to accept all .TIF files within that folder. This can be VERY useful when you have more than 8 files! Either way, it will take all supplied files and output the .VRT. In my case it saves in the same folder right next to my .TIF files.

From here, the process is the same as our previous gdal2tiles.py script except changing the names appropriately and using .vrt instead of .png. Here is another example:

Once again, rinse and repeat for your alternate Black Marble layer to get the final tileset, woo!

Note 7.1

In my exact instance, I utilized an amazing tool, gdal2tiles_parallel.py. It allows a system to utilize all of its cores in order to process tiles, and can save you immense amounts of time if you ever decide to process tiles through GDAL. It does have a slightly different setup however.
In my instance, I utilized a command like this:

This allows GDAL to maximize my computer’s resources with the applied necessary options to run the command in a similar manner to those above.

Step 8: Connecting the Data to a Webmap.


Again, numerous ways that this can be tackled. I personally use Leaflet as my preferred mapping library. I find it the lightest and easiest solution to work with. It will be what this post covers here on out.

Leaflet has many tutorials to get anyone started in web mapping, but even if you are not terribly familiar with it. Here is a download link for a bare-bones browser setup so you can see how all of the files (based on my personal storage preference) work together. You can load the leaflet map without any internet connection. Of note here, is that this template only contains levels 3-5 in order to vastly reduce the number of files to download/extract/work with. If you made your own tiles simply replace them.

If you take a look at the local files in this template there are only a few components. The .HTML file to load in a browser, a CSS folder with our style files, a JS folder for scripts, and the data folder for… well… data. The data folder is where we are going to store the tiles. In this project, the tiles are the only data we will be using, so this setup allows us to easily point to them. The template’s HTML code in its entirety looks this:

Note 8.1

A quirk with this code that is not common among other coders is that I personally use the variable MainMap as opposed to map. This has helped me personally with differentiating between ‘map’ when it is referring to the element ID or the variable within Javascript, especially when utilizing multiple maps (like with a reference map). It makes a big difference to me, but if you are copying code from other areas, it won’t be applied unless you ensure the MainMap variable is being matched or changed appropriately. It is significantly more common to see map as the variable when referencing a Leaflet code.

We have all of the basics of an HTML setup to include pulling local Leaflet dependencies: js/leaflet.js  and css/leaflet.css.

css/nightlight.css  is simply a short css setup that sets this styling:

These options are required for Leaflet to cover the entire browser window (without being in fullscreen mode). Just a convenience.

In the above HTML code lines 25-28 are highlighted since they show where we are linking the tile data. By using relative paths in our .HTML file we can ensure that the data is being pulled even from a local directory. As long as the tiles are in the same location relative to the .HTML file it will be able to load them, whether on your desktop or on a server. Leaflet understands this statement: L.tileLayer ('data/tiles/BlackMarble2012/{z}/{x}/{y}.png', {tms:true}); as a new layer where {z}, {x}, and {y} are variables that will search for the available tiles currently in Leaflet’s map extent. {tms:true} informs Leaflet that this is a TMS layer (the default output of gdal2tiles.py and its parallel branch). This is simply a slight specification difference from similar XYZ tiles where the ‘y’ value is inverted. One starts y’s value from top to bottom and the other from bottom to top. At this point we just repeat the process by updating to the appropriate folder path per layer.

The next section involves calling the actual map built by Leaflet’s library and rendering it in the HMTL:

Options such as center, zoom, minZoom, maxZoom, maxBounds are ones that I commonly use to keep users in place in case their panning/scrolling goes haywire and keeps the user in the area of interest. maxBounds is interesting in that a user can normally scroll indefinitely left or right and the map will loop back to where it was before. Tiles will load successfully if moving several “worlds” over, but many types of vector data will not load if you are not on the original “world.” I add some extra wiggle room on both sides so there isn’t a hard limit if the user is interested in viewing the edges in the middle of the screen. I cut off the northern and southern poles because the default Mercator projection has massive amounts of space for these areas, and in this instance it is a lot of unnecessary space.

So now all tiles are being connected and being called by specific variables and we utilize one of them to add to the map right away with layers:[BlackMarble2012]. Next we can add our difference layers with some simple addTo(MainMap); commands. If this seems odd see Note 8.1.
This is done by adding the highlighted lines below:

Now your map should have the new/receding lights on it as well, but the next major step is to add some user control to this. We need to be able to switch between 2012 and 2016 imagery as well as turn difference layers on and off so the user has as much control as possible.

This can easily be done by using Leaflet’s Layer Control. You can think of controls as widgets for a map that Leaflet provides some functionality for. This one allows users to turn off/switch between layers. In Leaflet there are a few types of layers and we will be working with basemaps and layers. Basemaps are typically listed first, and are designed to only have one on at a time. Layers are designed to be loaded as the user desires. For instance, there is no use loading 2012 imagery behind the 2016 imagery. Once both are loaded the 2012 imagery will not be visible, but will spend time loading and making the map slower. I prefer to separate them out using the additional lines of code below:

In this code I separate my basemap and layers into groups, specifically I use the variables LegendBasemaps  and LegendOverlays. These groups have the visible label in the control (I call it a legend for now), and the variable after the colon. Similarly to the GDAL command the quotes are necessary if your labels have spaces, but don’t hurt otherwise. Afterwords we add the control to the map with options: position: 'bottomright', collapsed:false to attach it to the bottom right of our map and to leave it open. If you refresh the map after adding these lines you can alter the layers all you like now. 😀

The next goal is to buff up the layer control a bit more by making it bit easier to understand, we can edit the font color and add a color ramp to make it more apparent how the the light intensities have changed. This can be easily added by adding some html into the label text. As an example we can be turn LegendOverlays into:

Here we use a <font>  tag to change the color and turn it bold, and a CSS gradient to give it a color bar to goes from fully transparent to the relative color. The various versions of the gradients exist to support multiple browsers.

At this point we have all of the critical pieces built for the project as displayed in my previous post. However there is a bit more polish we can add for purely aesthetic reasons.

Step 9: Polishing the Webmap.


I’ll be brief with these as this has already been a rather long post. Below I added a few convenience controls and the necessary stylings. As well as some third party layers, that help bring more information and aesthetics to the map if desired. Here is my completed code example with changes highlighted:

I also add a few styles to nightlights.css to make the default Leaflet elements darker and style other map elements. Here are the additional items added:

Now that the elements have a darker theme, we can add some useful content as well.

  • Lines 21-22 simply add a title banner in a futuristic style I’m fond of.
  • We add some plugins for convenience and smoother experience. In particular I often use the Default Extent plugin and Fade Basemap plugin. These guys did a great job making it easy to install these. You can see some styling for the Default Extent control in line 13, with the control actually added in line 80. The Fade Basemap plugin just simply needed to be appended to the document in line 30. You can grab the necessary CSS/JS files from these areas or check them out in the “finished template” at the bottom of this post.
  • I added two new layers from third parties (ESRI and CartoDB). ESRILabels and Dark Basemap instantly give the map more info and the ability to go for cleaner look. Both are called in lines 39 and 40 and added to their respective legend layers in lines 66 and 71. One trick about the ESRILabels is that since they aren’t loaded as vector data they, by default, will fall behind other tiles that are loaded based on order. They will always, by default, load behind vector data. To offset this, newer versions of Leaflet allow us to define custom panes with adjustable z-values. That is where lines 55-56 come into play (and the end of line 40). At the end of 40 we add it to a custom pane and in lines 55-56 we create that pane and set its z-index to be above other layers, but still under map elements.

I believe that covers all of the polishing changes. My hosted project has a few other things added that I feel are necessary for explaination and identification/branding but are of no use to anyone else, so I’ll leave those out. If you wanted to compare my “finished setup” you can download it here to mess around with it: Finished Template. (Remember for storage/processing sake these only have levels 3-5 and won’t load tiles when zoomed in further.)

That’s all there was to it. I think it seems like a lot after almost 6,000 words, but in actuality if you are familiar with html/css/js/gdal even a little bit this could be done from scratch in one or two days… and that’s mostly just downloading and processing time.

If you happened to have read through all that, I hope you get up and stretch. It was good to write about and I’m glad I got a chance to showcase something fun. I have a few older projects I’m dusting off and a fun desktop idea related to Westworld I’m entertaining. So my next project is likely to be considerably simpler instead of a wall of text. Till then, thank you for reading!


Alan Clack

Leave a Reply

Your email address will not be published. Required fields are marked *

Alan Clack © 2016