Open In Colab

1. DMSP-OLS annual composites in Google Earth Engine (5 min)

This module looks at the various ways to transform and aggregate nighttime lights data.

Two common processes are compositing and mosaicing.

In the context of remote sensing, the process of making a composite generally refers to taking a collection of spatially overlapping images and applying an aggregate function, also known as a "reducer" function, to these images to create a single composite image. For example, you may want to create an annual composite for a series of images by applying a reducer function, such as calculating the median, to the images to get a single image for the year that contains the median values per pixel.
Mosaicing refers to the process of assembling different spatially located images from the same time period together to create a seamless single image.

In this tutorial, we’ll identify and visualize annual composites in the DMSP-OLS nighttime lights series.

The DMSP-OLS image collection in Google Earth Engine, sourced by the Earth Observation Group, Payne Institute for Public Policy, Colorado School of Mines, has been processed as a series of annual composites per year, per satellite (some years include two satellite sensors), so we dont have to compose them ourselves. See Elvidget et al. (1997) [EBK+97] for the methdology.

Our tasks in this exercise:

  1. Load and inspect the DMSP-OLS nighttime lights Image Collection in GEE using geemap

  2. Search the reference table of DMSP-OLS satellites by year for particular annual composite (from 1996).

  3. Add the selected annual composite to a map.

  4. Add a 2nd annual composite from another year (2010) and create a slider panel to view and compare both.

1.1. Load and inspect the DMSP-OLS nighttime lights Image Collection

1.1.1. Initialize map object

First, let’s import geemap and initialize a geemap object centered on the greater Washington, DC area. We’ll also add the default satellite (daytime) basemap.

# reminder that if you are installing libraries in a Google Colab instance you will be prompted to restart your kernal

try:
    import geemap, ee
except ModuleNotFoundError:
    if 'google.colab' in str(get_ipython()):
        print("package not found, installing w/ pip in Google Colab...")
        !pip install geemap
    else:
        print("package not found, installing w/ conda...")
        !conda install mamba -c conda-forge -y
        !mamba install geemap -c conda-forge -y
    import geemap, ee
try:
        ee.Initialize()
except Exception as e:
        ee.Authenticate()
        ee.Initialize()

# set our initial map parameters for Washington, DC
center_lat = 38.9072
center_lon = -77.0369
zoomlevel=10

# initialize our map
map1 = geemap.Map(center=[center_lat,center_lon], zoom=zoomlevel)
map1.add_basemap('SATELLITE')

map1.addLayerControl()
map1

1.1.2. Get DMSP-OLS Image Collection and review meta-data and description

Recall the image collection for DMSP-OLS is located at: NOAA/DMSP-OLS/NIGHTTIME_LIGHTS

dmsp = ee.ImageCollection("NOAA/DMSP-OLS/NIGHTTIME_LIGHTS")

What is the total number of images in this collection? (It spans 1992-2013, but some years contain multiple images due to satellite overlap.)

In the Google Earth Engine Editor, we call the .size() method on our collection.

With geemap we can do the same thing, however, .size() will produce a size “object”, so we have to add the extra step of using the .getInfo() method so it prints out to our notebook.

In the GEE editor (in JavaScript):

print(dmsp.size());

In our Python notebook:

print(dmsp.size().getInfo())
35

TIP: We can use Python’s functional string method by adding f and brackets to dynamically print our collection size in a sentence:

print(f"There are {dmsp.size().getInfo()} images in this collection.")
There are 35 images in this collection.

What is the date range of our collection?

GEE has a set of methods called “Reducers,” which do a range of functions, such as get the sum or avg value of a collection. They are quite handy. The function Reducer.minMax() can be used to get a date range.

In the GEE editor (in JavaScript):

var imgrange = dmsp.reduceColumns(ee.Reducer.minMax(), ["system:time_start"]);
var start = ee.Date(imgrange.get('min'));
var end = ee.Date(imgrange.get('max'));
print('Date range: ', start, end);

Let’s try this in Python…

imgrange = dmsp.reduceColumns(ee.Reducer.minMax(), ["system:time_start"])
start = ee.Date(imgrange.get('min')).getInfo()
end = ee.Date(imgrange.get('max')).getInfo()
print(f"Date range: {start, end}")
Date range: ({'type': 'Date', 'value': 694224000000}, {'type': 'Date', 'value': 1356998400000})

Huh? One thing the GEE editor does is unpacks the JSON output (like a Python dictionary).

We can do that here by accessing the “value” key from our .getInfo() object. Let’s try it again:

imgrange = dmsp.reduceColumns(ee.Reducer.minMax(), ["system:time_start"])
start = ee.Date(imgrange.get('min')).getInfo()['value']
end = ee.Date(imgrange.get('max')).getInfo()['value']
print(f"Date range: {start, end}")
Date range: (694224000000, 1356998400000)

OK…these numbers still don’t make any sense. What’s happening?

1.1.3. Date conversions…

Tracking time in a world with many timezones is incredibly complicated. And you can imagine the nightmare scenario if applications (like airplane navigation systems) relied on arbitrary or inconsistent methods of tracking time. Therefore, software applications (including those that preprocessed our DMSP-OLS files as well as GEE) use what’s called Unix time, which is what we see here.

We'll spare you the hairy details, but basically, Unix time a.k.a. Epoch time refers to the time that has elapsed (not counting leap seconds) since the begining of the "Unix epoch," which started at 00 hrs 00 minutes 00 seconds on 1 January, 1970 in Coordinated Universal Time (UTC). There are two important conventions you should assume unless stated otherwise: 1) this number is measured in milliseconds (a positive integer after Jan 1, 1970, negative before) and 2) it is in UTC time (equivalent to Greenwich Mean Time).

It’s an extra step to convert this to a date that we humans understand in Python, but understanding timestamps is critical when working with temporal data – especially data produced by satelliltes orbiting the planet, so hopefully this helps you appreciate working with data that spans time and space!

To convert these Unix time values to readable dates, we’ll:

  1. Divide by 1000 to convert to seconds and

  2. Use Python’s handy datetime library to convert to datetime object using the utcfromtimestamp method (recall the time is in UTC!)

  3. Convert this datetime object to a string that we can read with the pattern: Year (%Y) - Month (%m) - Day (%d) Hour (%H): Minute (%M): Second (%S)

imgrange = dmsp.reduceColumns(ee.Reducer.minMax(), ["system:time_start"])
start = ee.Date(imgrange.get('min')).getInfo()['value']
end = ee.Date(imgrange.get('max')).getInfo()['value']

# convert date
from datetime import datetime
start = datetime.utcfromtimestamp(start/1000).strftime('%Y-%m-%d %H:%M:%S')
end = datetime.utcfromtimestamp(end/1000).strftime('%Y-%m-%d %H:%M:%S')
print(f"Date range: {start, end}")
Date range: ('1992-01-01 00:00:00', '2013-01-01 00:00:00')

OK! That makes sense.

So to recap…

In the GEE editor (in JavaScript):

var imgrange = dmsp.reduceColumns(ee.Reducer.minMax(), ["system:time_start"]);
var start = ee.Date(imgrange.get('min'));
var end = ee.Date(imgrange.get('max'));
print('Date range: ', start, end);

In Python:

from datetime import datetime

imgrange = dmsp.reduceColumns(ee.Reducer.minMax(), ["system:time_start"])
start = ee.Date(imgrange.get('min')).getInfo()['value']
end = ee.Date(imgrange.get('max')).getInfo()['value']

start = datetime.utcfromtimestamp(start/1000).strftime('%Y-%m-%d %H:%M:%S')
end = datetime.utcfromtimestamp(end/1000).strftime('%Y-%m-%d %H:%M:%S')
print(f"Date range: {start, end}")

If this is a method you’ll want to run more than once, you should create a function from these lines of code that only needs the image collection as a parameter.

def get_date_range(img_collection):
    imgrange = img_collection.reduceColumns(ee.Reducer.minMax(), ["system:time_start"])
    start = ee.Date(imgrange.get('min')).getInfo()['value']
    end = ee.Date(imgrange.get('max')).getInfo()['value']

    start = datetime.utcfromtimestamp(start/1000).strftime('%Y-%m-%d %H:%M:%S')
    end = datetime.utcfromtimestamp(end/1000).strftime('%Y-%m-%d %H:%M:%S')
    print(f"Date range: {start, end}")

1.1.4. Get DMSP-OLS annual composite for 1996

Recall from Introduction to nighttime light data (20 min):
NOAA`s National Center for Environmental Information has the reference table of DMSP satellites:

DMSP-OLS satellites by year

F10

F12

F14

F15

F16

F18

1992

F101992

1993

F101993

1994

F101994

F121994

1995

F121995

1996

F121996

1997

F121997

F141997

1998

F121998

F141998

1999

F121999

F141999

2000

F142000

F152000

2001

F142001

F152001

2002

F142002

F152002

2003

F142003

F152003

2004

F152004

F162004

2005

F152005

F162005

2006

F152006

F162006

2007

F152007

F162007

2008

F162008

2009

F162009

2010

F182010

2011

F182011

2012

F182012

2013

F182013

In 1996, we have but one option to choose: F121996

1.2. Add this composite to our map

Lets add this image as a layer to our map.

Let’s also apply the mask to areas with no data and adjust the opacity to 75%.

Recall how to do this from Practical exercise: image visualization (10 min)

dmsp1996 = ee.Image("NOAA/DMSP-OLS/NIGHTTIME_LIGHTS/F101992")

# initialize our map
map2 = geemap.Map(center=[center_lat,center_lon], zoom=zoomlevel)
map2.add_basemap('SATELLITE')

map2.addLayer(dmsp1996.mask(dmsp1996), {}, "DMSP-OLS 1996", opacity=0.75)

map2.addLayerControl()
map2

1.3. Add a 2nd annual composite from another year (2010) and create a slider panel to view and compare both.

1.3.1. Warning, this is based on ipyleaflet a Python library that does not play well with Google Colab, so this code will not work in the Google Colab environment but should on your local machine.

# get the satellite name from the reference table
dmsp2010 = ee.Image("NOAA/DMSP-OLS/NIGHTTIME_LIGHTS/F182010")

map2.addLayer(dmsp2010.mask(dmsp2010), {}, "DMSP-OLS 2010", opacity=0.75)

Recall how to create a split panel slider from Practical exercise: image visualization (10 min):

Warning, this is based on `ipyleaflet` a Python library that does not play well with Google Colab, so the split panel display will not work in the Google Colab environment but should on your local machine.
# initialize our map
map3 = geemap.Map(center=[center_lat,center_lon], zoom=zoomlevel)
map3.add_basemap('SATELLITE')

# generate tile layers
dmsp1996_tile = geemap.ee_tile_layer(dmsp1996.mask(dmsp1996), {}, 'DMSP-OLS 1996', opacity=0.75)
dmsp2010_tile = geemap.ee_tile_layer(dmsp2010.mask(dmsp2010), {}, 'DMSP-OLS 2010', opacity=0.75)

# create split map
map3.split_map(left_layer=dmsp1996_tile, right_layer=dmsp2010_tile)
map3

Great! Now you’ve found and visualized an annual composite as well as learned a few helpful ways to access meta-data. This slider is cool as a visualization, but as noted in the first tutorial, the DMSP-OLS satellites do not have on-board calibration. As a result, the change from satellite to satellite (as is the case when comparing 1996 to 2010) might include sensor variations that do not represent actual observed changes in light.

This is why you should be very careful when comparing or analyzing change in DMSP-OLS series directly.

The way to address this is by adjusting our DMSP-OLS annual composites through a process we’ll call “intercalibration” and this is the subject of a later tutorial.

1.4. References:

EBK+97

Christopher D Elvidge, Kimberly E Baugh, Eric A Kihn, Herbert W Kroehl, and Ethan R Davis. Mapping city lights with nighttime data from the dmsp operational linescan system. Photogrammetric Engineering and Remote Sensing, 63(6):727–734, 1997.