Basics of nighttime lights data

NOAA provides tif files of the nightlights images. Images in the package are represented as 2D arrays with floating-point values. Images of different months are stacked together to form 3D arrays. In this package, such 3D arrays are called datacubes. The following examples demonstrate how array indices work in Julia in the context of this package.

image[1, 2]

Returns the value of the image at location [1, 2]. 1st row and 2nd column.

datacube[:, :, 3]

Returns the image of the 3rd month.

datacube[1, 2, :]

Returns the time series values of the pixel at location [1, 2].

datacube[1, 2, 3]

Returns the value of the image at location [1, 2] of the 3rd month.

Loading and saving files

Images and datacubes can be be loaded and saved using the following functions.

NighttimeLights.load_img — Method

NOAA provides nighttime lights images as tif files. They can be opened as 2D arrays using the load_img function.

load_img("example.tif")

NighttimeLights.save_img — Method

Images in the form of 2D arrays can be saved as tif files.

save_img("example.tif", img)

NighttimeLights.load_datacube — Method

NOAA provides images for each month since April 2012. Images of the same place taken at different times can be stacked together to make a 3D array representating a panel data. In this package, we refer to such arrays as datacubes. JLD files containing 3D arrays can be loaded using the load_datacube function.

load_datacube("example.jld")

NighttimeLights.make_datacube — Method

Loads all images (tif files) in a folder and generates a datacube. The function prints the file names to you the order in which they are loaded.

make_datacube("~/Downloads/ntl_images")

NighttimeLights.save_datacube — Method

3D arrays can be saved as JLD files.

save_datacube("example.jld", datacube)

Mapping between earth and arrays

Suppose you want to find which location has the maximum value of light in an image. You can use the findmax function in julia.

findmax(my_image)

Suppose, the answer is [2000, 3000]. If you want to find the coordinates of this location, you'll need a mapping between the earth's coordinates and the image's indices.

Similarly, if you were given a pair of latitude and longitude, for example, (19.05, 73.01), and you need to find the radiance of that coordinate (approximately), you'll need to convert these number to your image's indices.

NighttimeLights.CoordinateSystem — Type

The mapping between the coordinates of earth and indices of the image is needed to convert from one to another. In this package such mapping is referred to as a coordinate system. To define a coordinate system, you need to provide coordinates of the top-left and bottom-right pixels, and the height and the width of the image.

Suppose the image we plan to load had 8000 rows and 7100 columns. The top-left coordinate (37.5, 67.91666) is and bottom-right (4.166, 97.5). We can use these to create a the mapping.

top_left     = Coordinate(37.5, 67.91666)
bottom_right = Coordinate(4.166, 97.5)
height = 8000
width  = 7100
my_coordinate_system = CoordinateSystem(top_left, bottom_right, height, width)

Now functions can use this mapping to go from one to another.

NighttimeLights.lat_to_row — Method

To convert from latitude to row number of an image, use the lat_to_row function.

Example:

lat_to_row(my_coordinate_system, 19.6)

NighttimeLights.long_to_column — Method

To convert from longitude to column number of an image use the long_to_column function.

Example:

long_to_column(my_coordinate_system, 73.33)

NighttimeLights.row_to_lat — Method

To convert from row number of an image to latitude use the row_to_lat function.

Example:

row_to_lat(my_coordinate_system, 100)

NighttimeLights.column_to_long — Method

To convert from column number of an image to longitude use the column_to_long function.

Example:

column_to_long(my_coordinate_system, 100)

NighttimeLights.coordinate_to_image — Method

To convert a coordinate from earth's coordinate system to image's row number and column number, use the coordinate_to_image function.

Example:

coordinate_to_image(my_coordinate_system, Coordinate(19.49907, 72.721252))

NighttimeLights.image_to_coordinate — Method

To convert a coordinate from image's row number and column number to earth's coordinate system use the image_to_coordinate function.

Example:

image_to_coordinate(my_coordinate_system, [4320, 1153])

Masks

Masks are 2D arrays consisting of 0s and 1s. The 1s determine the region of interest. The pixels with the value of 1 are referred to as lit and the ones with the value of 0 are referred to as dark. Masks are useful because they can be easily combined with one another.

To demonstrate the concept of mask, here are 3 examples:

If the region of interest is India, the points inside the boundary of India will be 1, while the remaining will be 0.
If all pixels in an image below a threshold are considered background noise, such pixels can be marked as zero and the remaining can be marked as 1 to produce a background noise mask. For following code demonstrates this example.

image = rand(0:10.0, 10, 10)
noise_threshold = function(image, threshold)
    if x < threshold
        return 0
    else 
        return 1
    end
end
threshold = 0.3
mask_mask = noise_threshold.(image, threshold)

The standard deviation of each pixel in a datacube can be computed. Suppose those pixels with standard deviation greater than a threshold are considered outliers. In a 2D array, these pixels can be marked as 0 and the remaining can be marked as 1 to generate an outlier mask.

datacube = rand(0:10.0, 10, 10, 10)
mask = zeros(10, 10)
std_threshold = 1
for i in 1:10
    for j in 1:10
        if std(datacube[i, j, :]) > std_threshold
            mask[i, j] = 1
        end
    end
end

If you multiply (elementwise) the 3 masks, you get a mask that represents the pixels above the threshold, which are inside India and which are not outliers.

To find the number of pixels in a mask, one simply needs to do:

sum(mask)

To find the area of the mask, the mask_area function of the package can be used.

NighttimeLights.mask_area — Method

The area of each pixel is added to determine the total area of a mask. The res parameter is needed for the resolution in arc seconds. By default the resolution is 30 arc seconds.

mask = rand(0:1,8000,7100) 
mask_area(INDIA_COORDINATE_SYSTEM, mask)

An image can be multiplied with a mask (elementwise) so that only the pixels lit in the mask are lit in the images. For example:

image .* noise_mask

Returns as image with 0s wherever there is background noise and the remaining values are the same as the original image.

Aggregating over masks

While using nighttime lights, you may need to find the total lit of an image over a mask.

For example, if you have a background noise mask (where pixels considered noise are marked as 0 and the remaining at marked as 1), you may need to find the the total of an image over the lit pixels of the mask.

NighttimeLights.aggregate — Method

The aggregate value of an image over a mask can be computed by the aggregate function. The function multiplies the image and the mask (elementwise) and then performs the summation.

rand_image = rand(10, 10)
rand_mask = rand(0:1, 10, 10)
aggregate(rand_image, rand_mask)

The aggregate function is equivalent to

sum(image .* mask)

NighttimeLights.aggregate_timeseries — Method

To find the time series of aggregate values of a datacube over a mask, using the aggregate_timeseries function

rand_datacube = rand(10, 10, 10)
rand_mask = rand(0:1, 10, 10)
aggregate_timeseries(rand_datacube, rand_mask)