2. Calculate rate of change (5 min)¶

Nighttime lights can be useful for change detection. A basic measure of change over time is to calculate the per-pixel rate of change or the slope (sometimes referred to as the “Slope of Change”), which measures the change in rise (measures of radiance) over the change in run (time).

In this short exercise, we’ll calculate the rate of change looking at examples for both DMSP-OLS and VIIRS-DNB.

Our tasks in this exercise:

Calculate the rate of change in DMSP-OLS from 2000 to 2007
Calculate the rate of change in VIIRS-DNB from January 2014 to January 2020

2.1. DMSP-OLS rate of change 2000 to 2007¶

2.1.1. Selecting the satellite¶

You’ll recall from previous tutorials DMSP-OLS annual composites in Google Earth Engine (5 min) that the DMSP-OLS data are available in GEE as annual composites.

We’ll get that ImageCollection and filter it on our dates, 2000 to 2007.

# 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()

# retrieve DMSP-OLS data, filtered to 2000 to 2007 (inclusive)
dmsp = ee.ImageCollection("NOAA/DMSP-OLS/NIGHTTIME_LIGHTS").filterDate('2000-01-01','2007-12-31')

print(f"For 8 years, we have a time series of {dmsp.size().getInfo()} annual composites.")

For 8 years, we have a time series of 16 annual composites.

Wait, if we’re using annual composites, we do we have 16 images for 8 years [2000 thru 2007]?

Recall the DMSP program has a series of satellites that overlapped. Introduction to nighttime light data (20 min).

Intercalibration: since the different satellite sensors are not calibrated on-board (also discussed in Module 1), it is very ill-advised to compare DMSP-OLS data across satellites without performing calibration. We did this in DMSP-OLS intercalibration (10 min), but for this tutorial we’ll just stick with a single satellite.

Satellite F15 was operational from 2000 to 2007 so we’ll focus on data it collected for now. Satellites F14 and F16 also overlapped those years, which is why we have more images than years.

We’ll filter our collection to only include images from F15.

We can do this by using the .filterMetadata() function on our ImageCollection to filter for images that contain only F15 in their name (you’ll recall the naming convention for DMSP-OLS includes the satellite name as prefix for the image.

# retrieve DMSP-OLS data, filtered to 2000 to 2007 (inclusive) AND filtered for F15 satellite images only
# we filter on the property name "system:index" - image name -- and get those that contain "F15"
dmsp = ee.ImageCollection("NOAA/DMSP-OLS/NIGHTTIME_LIGHTS").filterDate(
    '2000-01-01','2007-12-31').filterMetadata(
    'system:index','contains','F15')

print(f"For {2008-2000} years, we have a time series of {dmsp.size().getInfo()} annual composites.")

For 8 years, we have a time series of 8 annual composites.

Excellent! Now we have only F15 images for 2000-2007.

2.1.1.1. On your own¶

Since we know how to perform basic intercalibration, you may want to try calibrating a series of data across satellites…

2.1.2. Rise over run¶

The rate of change is the slope, i.e. rise/run. Our run (in years) is 8. Our rise will be the difference from in Digital Number (DN) values for the last image relative to the first.

To calculate the rise, we’ll get the first and last images in our collection. We can explicitly select these since we know the names (“F152000” and “F152007”), but we can also do so programmatically, which is good practice so this operation can be generalized to other series without needing to reference the exact image name.

As pulled from GEE, our ImageCollection should be sorted by date in ascending order already, but we’ll sort it anyway since it’s a good idea never to assume any order or structure to data blindly!

# sort by image "time_end"
first_img = dmsp.sort('system:time_end').first()

# reverse sort so that last=first
last_img = dmsp.sort('system:time_end',False).first()

Now we subtract our initial values (first image) from our last to get the gain (or loss) in 2007 relative to 2000, and divide by our run of 8 years to get the annualized rate of change.

We can use simple math operations as we learned in Cell statistics and basic band math (10 min).

Note A: we’ll select the stable_lights band for this calculation.

First we’ll get the difference and then we’ll divide by 8 to get the annualized rate.

Note that in this form, the operations are performed sequentially, not according to the order of ops. So the difference will be calculated first and then the quotient from that.

It doesnt hurt to set that first operation in parenthesis anyway (last minus first); however, so that the order of operations is explicitly clear to anywone following this.

dmsp_slope = (last_img.select('stable_lights').subtract(first_img.select('stable_lights'))).divide(8)

2.1.3. Let’s initialize a map centered on Shanghai, China.¶

We’ll add our rate of change layer.

We’ll also set some visualization parameters to help the view. We’ll set a color palette of red <-> blue (red=positive change, blue=negative change). We know the min / max DN value is -63 and 63 respectively, and since we’re dividing by 8, we’ll scale this by 8 and set our min/max to -8 and 8 respectively.

This will “stretch” our color gradient so we can see the change more clearly.

# center on Shanghai
lat = 31.18
lon = 121.49

dmspMap = geemap.Map(center=[lat,lon],zoom=8)

viz_params = {'min':-8,
             'max':8,
             'palette':['1d4877','1b8a5a','f68838','ee3e32']}

# make it opaque so we can see underlying basemap
dmspMap.addLayer(dmsp_slope, viz_params, '2000-2007 DMSP-OLS annual rate of change',opacity=.66)
dmspMap.addLayerControl()
dmspMap

The red/orange indicates increased lights from 2000 to 2007, whereas green/blue indicates relative decrease. In this color scheme, zero change is green.

2.1.4. DMSP-OLS saturation¶

Note that the city center is green – this is likely because we know that the DMSP-OLS data is limited and easily saturated. So bright areas are “maxed” out at a DN of 63 in both 2000 as well as 2007. If that’s the case, we wouldnt see any change. So it’s more likely that the urban center of Shanghai is not showing change because of saturation and not because there actually was no true change. This is a known limitation to DMSP-OLS data (even after calibration), and something to be wary of.

But you can certainly see the increase in lights in the areas around the city and suburbs!

2.2. VIIRS-DNB rate of change 2014 to 2020¶

Now let’s calculate the rate of change for the VIIRS-DNB data. These data are calibrated, so we can more confidently compare across years without additional adjustments. These data are also provided in GEE in monthly composites.

For a simple rate of change calculation, this doesnt change much and we can show the rate in months instead of years.

First, we’ll get the collection filtered on our dates using the stray-light collected monthly composites.

We only have one satellite series, so overlapping data is not an issue for the VIIRS collection.

viirs = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG").filterDate(
    '2014-01-01','2020-01-31')

print(f"For {(2020-2014)*12+1} months, we have a time series of {viirs.size().getInfo()} monthly composites.")

For 73 months, we have a time series of 73 monthly composites.

Now we’ll get the first image, the last image and calculate our monthly rate of change.

We’ll use the avg_rad band.

# sort by image "time_end"
first_img = viirs.sort('system:time_end').first()

# reverse sort so that last=first
last_img = viirs.sort('system:time_end',False).first()

# get rate of change (diff over # months: 73)
viirs_slope = (last_img.select('avg_rad').subtract(first_img.select('avg_rad'))).divide(73)

And we’ll visualize this with the same color palette as DMSP-OLS, but scale the min/max values based on our VIIRS units (radiance vis-a-vis Watts/cm2/sr) and the divisor (73 months) in mind: min=-1,max=1

viirsMap = geemap.Map(center=[lat,lon],zoom=8)

viz_params = {'min':-1,
             'max':1,
             'palette':['1d4877','1b8a5a','f68838','ee3e32']}

# make it opaque so we can see underlying basemap
viirsMap.addLayer(viirs_slope, viz_params, '2014-2020 VIIRS-DNB monthly rate of change',opacity=.75)
viirsMap.addLayerControl()
viirsMap

Cool!

We already see a couple noticeable things:

The spatial resolution of VIIRS is much better than DMSP-OLS. We can actually see change at the neighborhood level, including road infrastructure!
Saturation is not an issue for VIIRS-DNB data (if you see anything “maxed” our here, that is just because of our choice in min/max values for visualization purposes). We’re able to see dynamics in nighttime lights almost everywhere.

Also note: The areas outside of our Area of Interest (AOI) are a bit distracting — and may be noise to filter for analysis, so if you find a geometry for your AOI (perhaps a boundary for Greater Shanghai in this case), you may want to clip your nighttime lights layer.

2.3. Bonus, on your own: quantifying the rate of change¶

We’ve calculated the rate of change and visualized this. But you can imagine that quantifying into numbers you can analyze can provide a lot of useful information, such as comparing the growth of two regions over the same period of time – or the growth of a single region compared across two separate periods of time.

To do this, you’ll want to calculate some aggregate measure or zonal statistic, such as the Sum of Lights (SOL) for a region. Can you do this on your own for a specific region and get the rate of change in SOL for a region?

Everything you need to do this you’ve already learned. HINT: check out Cell statistics and basic band math (10 min)

Open Nighttime Lights