SEDAC Gridded World Population Density

Global, 1 km resolution human population density estimates based on national censuses and population registers, version 4.11.
Author

Siddharth Chaudhary, Vishal Gaur

Published

June 29, 2023

Access this Notebook

You can launch this notebook in the US GHG Center JupyterHub (requires access) by clicking the following link: SEDAC Gridded World Population Density. If you are a new user, you should first sign up for the hub by filling out this request form and providing the required information.

If you do not have a US GHG Center Jupyterhub account, you can access this notebook through MyBinder by clicking the button below.

Binder

Table of Contents

Data Summary and Application

  • Spatial coverage: Global
  • Spatial resolution: 30 arc-seconds (~1 km at equator)
  • Temporal extent: 2000 - 2020
  • Temporal resolution: Annual, every 5 years
  • Unit: Number of persons per square kilometer (persons/km²)
  • Utility: Climate Research

For more information, visit the SEDAC Gridded World Population Density data overview page.

Approach

  1. Identify available dates and temporal frequency of observations for the given collection using the US Greenhouse Gas Center (GHGC) Application Programming Interface (API) /stac endpoint. The collection processed in this notebook is SEDAC Gridded World Population Density.
  2. Pass the STAC item into raster API /collections/{collection_id}/items/{item_id}/{tile_matrix_set_id}/tilejson.json endpoint
  3. Using folium.plugins.DualMap, visualize two tiles (side-by-side), allowing time point comparison.
  4. After the visualization, we’ll perform zonal statistics for a given polygon.
  5. Create a time-series analysis.

About the Data

SEDAC Gridded World Population Density

The SEDAC Gridded Population of the World: Population Density, v4.11 dataset provides annual estimates of population density for the years 2000, 2005, 2010, 2015, and 2020 on a 30 arc-second (~1 km) grid. These data can be used to assess disaster impacts, risk mapping, and any other applications that include human dimension. This population density dataset is provided by NASA’s Socioeconomic Data and Applications Center (SEDAC), hosted by the Center for International Earth Science Information Network (CIESIN) at Columbia University. The population estimates are provided as a continuous raster for the entire globe.

For more information regarding this dataset, please visit the SEDAC Gridded World Population Density data overview page.

Terminology

Navigating data via the US Greenhouse Gas Center (GHGC) Application Programming Interface (API), you will encounter terminology that is different from browsing in a typical filesystem. We’ll define some terms here which are used throughout this notebook.

  • catalog: All datasets available at the /stac endpoint
  • collection: A specific dataset, e.g. SEDAC Gridded World Population Density
  • item: One data file (i.e. granule) in the dataset, e.g. one annual file of population density
  • asset: A variable available within the granule, e.g. population density
  • STAC API: SpatioTemporal Asset Catalogs - Endpoint for fetching metadata about available datasets
  • Raster API: Endpoint for fetching data itself, for imagery and statistics

Install the Required Libraries

Required libraries are pre-installed on the US GHG Center Hub. If you need to run this notebook elsewhere, please install them with this line in a code cell:

%pip install requests folium rasterstats pystac_client pandas matplotlib –quiet

# Import the following libraries
# For fetching from the Raster API
import requests
# For making maps
import folium
import folium.plugins
from folium import Map, TileLayer
# For talking to the STAC API
from pystac_client import Client
# For working with data
import pandas as pd
# For making time series
import matplotlib.pyplot as plt
# For formatting date/time data
import datetime
# Custom functions for working with GHGC data via the API
import ghgc_utils

Query the STAC API

STAC API Collection Names

Now, you must fetch the dataset from the STAC API by defining its associated STAC API collection ID as a variable. The collection ID, also known as the collection name, for the SEDAC Gridded World Population Density dataset is sedac-popdensity-yeargrid5yr-v4.11.

# Provide the STAC and RASTER API endpoints
# The endpoint is referring to a location within the API that executes a request on a data collection nesting on the server.

# The STAC API is a catalog of all the existing data collections that are stored in the GHG Center.
STAC_API_URL = "https://earth.gov/ghgcenter/api/stac"

# The RASTER API is used to fetch collections for visualization
RASTER_API_URL = "https://earth.gov/ghgcenter/api/raster"

# The collection name is used to fetch the dataset from the STAC API. First, we define the collection name as a variable
# Name of the collection for SEDAC population density dataset 
collection_name = "sedac-popdensity-yeargrid5yr-v4.11"
# Using PySTAC client
# Fetch the collection from the STAC API using the appropriate endpoint
# The 'pystac' library makes an HTTP request 
catalog = Client.open(STAC_API_URL)
collection = catalog.get_collection(collection_name)

# Print the properties of the collection to the console
collection
<CollectionClient id=sedac-popdensity-yeargrid5yr-v4.11>

Examining the contents of our collection under summaries we see that the data is available from January 2000 to December 2020. By looking at the dashboard:time density we observe that the data is available for the years 2000, 2005, 2010, 2015, 2020.

items = list(collection.get_items())  # Convert the iterator to a list
print(f"Found {len(items)} items")
Found 5 items
# Examine the first item in the collection
# Keep in mind that a list starts from 0, 1, 2... therefore items[0] is referring to the first item in the list/collection
items[0]
<Item id=sedac-popdensity-yeargrid5yr-v4.11-gpw_v4_population_density_rev11_2020_30_sec_2020> 
# Restructure our items into a dictionary where keys are the datetime items
# Then we can query more easily by date/time, e.g. "2020"
items_dict = {item.properties["start_datetime"][:7]: item for item in collection.get_items()}
# Before we go further, let's pick which asset to focus on for the remainder of the notebook.
# This dataset only has one available asset:
asset_name = "population-density"

Creating Maps Using Folium

You will explore changes in population density in urban regions and visualize the outputs on a map using folium.

Fetch Imagery Using Raster API

Here we get information from the Raster API, which we will add to our map in the next section.

# Specify two date/times that you would like to visualize, using the format of items_dict.keys()
dates = ["2000-01","2020-01"]

Below, we use some statistics of the raster data to set upper and lower limits for our color bar. These are saved as rescale_values, and will be passed to the Raster API in the following step(s).

# Extract collection name and item ID for the first date
first_date = items_dict[dates[0]]
collection_id_1 = first_date.collection_id
item_id_1 = first_date.id
# Select relevant asset
object = first_date.assets[asset_name]
raster_bands = object.extra_fields.get("raster:bands", [{}])
# Print raster bands' information
raster_bands
[{'scale': 1.0,
  'nodata': -9999.0,
  'offset': 0.0,
  'sampling': 'area',
  'data_type': 'float32',
  'histogram': {'max': 20757.5234375,
   'min': -800.1041259765625,
   'count': 11,
   'buckets': [129058, 321, 40, 21, 10, 3, 1, 0, 1, 2]},
  'statistics': {'mean': 42.60443622206601,
   'stddev': 234.09900050691866,
   'maximum': 20757.5234375,
   'minimum': -800.1041259765625,
   'valid_percent': 24.69196319580078}}]
rescale_values = {
    "max": raster_bands[0]["statistics"]["maximum"],
    "min": 0,
}

print(rescale_values)
{'max': 20757.5234375, 'min': 0}

Now, you will pass the item id, collection name, asset name, and the rescale values to the Raster API endpoint, along with a colormap. This step is done twice, one for each date/time you will visualize, and tells the Raster API which collection, item, and asset you want to view, specifying the colormap and colorbar ranges to use for visualization. The API returns a JSON with information about the requested image. Each image will be referred to as a tile.

# Choose a colormap for displaying the data
# Make sure to capitalize per Matplotlib standard colormap names
# For more information on Colormaps in Matplotlib, please visit https://matplotlib.org/stable/users/explain/colors/colormaps.html
color_map = "rainbow" 
# Make a GET request to retrieve information for the date mentioned below
tile_matrix_set_id = "WebMercatorQuad"

month1_tile = requests.get(
    f"{RASTER_API_URL}/collections/{collection_id_1}/items/{item_id_1}/{tile_matrix_set_id}/tilejson.json?"
    f"&assets={asset_name}"
    f"&color_formula=gamma+r+1.05&colormap_name={color_map}"
    f"&rescale={rescale_values['min']},{rescale_values['max']}"
).json()

# Print the properties of the retrieved granule to the console
month1_tile
{'tilejson': '2.2.0',
 'version': '1.0.0',
 'scheme': 'xyz',
 'tiles': ['https://earth.gov/ghgcenter/api/raster/collections/sedac-popdensity-yeargrid5yr-v4.11/items/sedac-popdensity-yeargrid5yr-v4.11-gpw_v4_population_density_rev11_2000_30_sec_2000/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?assets=population-density&color_formula=gamma+r+1.05&colormap_name=rainbow&rescale=0%2C20757.5234375'],
 'minzoom': 0,
 'maxzoom': 24,
 'bounds': [-180.0, -90.0, 179.99999999999983, 89.99999999999991],
 'center': [-8.526512829121202e-14, -4.263256414560601e-14, 0]}
# Repeat the above for your second date/time
# Note that we do not calculate new rescale_values for this tile
# We want date tiles 1 and 2 to have the same colorbar range for best visual comparison.
second_date = items_dict[dates[1]]
collection_id_2 = second_date.collection_id
item_id_2 = second_date.id

month2_tile = requests.get(
    f"{RASTER_API_URL}/collections/{collection_id_2}/items/{item_id_2}/{tile_matrix_set_id}/tilejson.json?"
    f"&assets={asset_name}"
    f"&color_formula=gamma+r+1.05&colormap_name={color_map}"
    f"&rescale={rescale_values['min']},{rescale_values['max']}"
).json()

# Print the properties of the retrieved granule to the console
month2_tile
{'tilejson': '2.2.0',
 'version': '1.0.0',
 'scheme': 'xyz',
 'tiles': ['https://earth.gov/ghgcenter/api/raster/collections/sedac-popdensity-yeargrid5yr-v4.11/items/sedac-popdensity-yeargrid5yr-v4.11-gpw_v4_population_density_rev11_2020_30_sec_2020/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?assets=population-density&color_formula=gamma+r+1.05&colormap_name=rainbow&rescale=0%2C20757.5234375'],
 'minzoom': 0,
 'maxzoom': 24,
 'bounds': [-180.0, -90.0, 179.99999999999983, 89.99999999999991],
 'center': [-8.526512829121202e-14, -4.263256414560601e-14, 0]}

Generate Map

# Initialize the map, specifying the center of the map and the starting zoom level.
# 'folium.plugins' allows mapping side-by-side via 'DualMap'
# Map is centered on the position specified by "location=(lat,lon)"
map_ = folium.plugins.DualMap(location=(38.97, -77), zoom_start=9) #Los Angeles

# Define the first map layer
map_layer_1 = TileLayer(
    tiles=month1_tile["tiles"][0], # Path to retrieve the tile
    attr="US GHG Center", # Set the attribution
    opacity=0.8, # Adjust the transparency of the layer
    overlay=True,
    name=f"{collection.title}, {dates[0]}"
)

# Add the first layer to the Dual Map
map_layer_1.add_to(map_.m1)

# Define the second map layer
map_layer_2 = TileLayer(
    tiles=month2_tile["tiles"][0], # Path to retrieve the tile
    attr="US GHG Center", # Set the attribution
    opacity=0.8, # Adjust the transparency of the layer
    overlay=True,
    name=f"{collection.title}, {dates[1]}"
)

# Add the second layer to the Dual Map
map_layer_2.add_to(map_.m2)

# Add a layer control to switch between map layers
folium.LayerControl(collapsed=False).add_to(map_)

# Add colorbar
# First we'll rescale our numbers to make the colorbar labels nicer
re_rescale_values = {
    "min":rescale_values["min"]/1e3,
    "max":rescale_values["max"]/1e3
}
# We can use 'generate_html_colorbar' from the 'ghgc_utils' module 
# to create an HTML colorbar representation.
legend_html = ghgc_utils.generate_html_colorbar(color_map,re_rescale_values,label=f'{items[0].assets[asset_name].title} (10^3 persons/km2)',dark=True)

# Add colorbar to the map
map_.get_root().html.add_child(folium.Element(legend_html))

# Visualize the Dual Map
map_
Make this Notebook Trusted to load map: File -> Trust Notebook

This comparison indicates a slight increase in population density in Washington, D.C. from 2000 to 2020.

Calculate Zonal Statistics

To perform zonal statistics, first we need to create a polygon. In this use case we are creating a polygon around Washington D.C.

# Give the area of interest (AOI) a name, to be used in plotting later
aoi_name = "Washington, D.C."
# Define the AOI as a GeoJSON
aoi = {
    "type": "Feature", # Create a feature object
    "properties": {},
    "geometry": { # Set the bounding coordinates for the polygon
        "coordinates": [
            [
                [-77.25, 39],    # Northwest bounding coordinate
                [-77.25, 38.8],  # Southwest bounding coordinate
                [-76.8, 38.8],   # Southeast bounding coordinate
                [-76.8, 39],     # Northeast bounding coordinate
                [-77.25, 39]     # Northwest bounding coordinate (closing the polygon)
            ]
        ],
        "type": "Polygon",
    },
}
# Quick Folium map to visualize this AOI
map_ = folium.Map(location=(38.9, -77.1), zoom_start=10)
# Add AOI to map
folium.GeoJson(aoi, name=aoi_name,style_function=lambda feature: {"fillColor":"none","color":"#E467E8"}).add_to(map_)
# Add data layer (created in previous section) to visualize number of grid cells within AOI
map_layer_1.add_to(map_)
# Add a layer control to switch between map layers
folium.LayerControl(collapsed=False).add_to(map_)
# Add a quick colorbar
legend_html = ghgc_utils.generate_html_colorbar(color_map,re_rescale_values,label='10^3 persons/km2',dark=True)
map_.get_root().html.add_child(folium.Element(legend_html))
map_
Make this Notebook Trusted to load map: File -> Trust Notebook

We will perform zonal statistics and create a dataframe using our AOI for Washington, D.C., the items from our initial catalog search, and the “population-density” asset.

We can generate the statistics for the AOI using a function from the ghgc_utils module, which fetches the data and its statistics from the Raster API.

%%time
# %%time = Wall time (execution time) for running the code below

# Statistics will be returned as a Pandas DataFrame
df = ghgc_utils.generate_stats(items,aoi,url=RASTER_API_URL,asset=asset_name)
# Print first five rows of our statistics DataFrame
df.head(5)
Generating stats...
Done!
CPU times: user 279 ms, sys: 16.3 ms, total: 296 ms
Wall time: 3.97 s
datetime min max mean count sum std median majority minority unique histogram valid_percent masked_pixels valid_pixels percentile_2 percentile_98 date
0 2020-01-01T00:00:00+00:00 0.00000000000000000000 18091.92382812500000000000 2342.92065429687500000000 1296.00000000000000000000 3036425.25000000000000000000 2331.32473070569903939031 1712.98815917968750000000 0.00000000000000000000 0.07085865736007690430 1261.00000000000000000000 [[682, 351, 143, 57, 38, 13, 5, 2, 2, 3], [0.0... 100.00000000000000000000 79.00000000000000000000 1296.00000000000000000000 0.00000000000000000000 9013.31640625000000000000 2020-01-01 00:00:00+00:00
1 2015-01-01T00:00:00+00:00 0.00000000000000000000 17640.75195312500000000000 2253.18554687500000000000 1296.00000000000000000000 2920128.50000000000000000000 2254.50648701661521045025 1643.98083496093750000000 0.00000000000000000000 0.06709918379783630371 1261.00000000000000000000 [[690, 351, 138, 55, 37, 13, 5, 2, 2, 3], [0.0... 100.00000000000000000000 79.00000000000000000000 1296.00000000000000000000 0.00000000000000000000 8628.78125000000000000000 2015-01-01 00:00:00+00:00
2 2010-01-01T00:00:00+00:00 0.00000000000000000000 17200.83203125000000000000 2167.17700195312500000000 1296.00000000000000000000 2808661.50000000000000000000 2180.98211822105486135115 1570.85021972656250000000 0.00000000000000000000 0.06353917717933654785 1261.00000000000000000000 [[701, 344, 138, 53, 36, 12, 5, 2, 2, 3], [0.0... 100.00000000000000000000 79.00000000000000000000 1296.00000000000000000000 0.00000000000000000000 8385.78710937500000000000 2010-01-01 00:00:00+00:00
3 2005-01-01T00:00:00+00:00 0.00000000000000000000 16771.88281250000000000000 2084.72875976562500000000 1296.00000000000000000000 2701808.50000000000000000000 2110.57326335761217706022 1506.44372558593750000000 0.00000000000000000000 0.06016804650425910950 1261.00000000000000000000 [[714, 341, 129, 53, 35, 13, 4, 2, 2, 3], [0.0... 100.00000000000000000000 79.00000000000000000000 1296.00000000000000000000 0.00000000000000000000 8039.70947265625000000000 2005-01-01 00:00:00+00:00
4 2000-01-01T00:00:00+00:00 0.00000000000000000000 16353.62988281250000000000 2005.68249511718750000000 1296.00000000000000000000 2599364.50000000000000000000 2043.11263761937516392209 1428.24743652343750000000 0.00000000000000000000 0.05697577446699142456 1261.00000000000000000000 [[722, 337, 131, 49, 35, 11, 4, 2, 2, 3], [0.0... 100.00000000000000000000 79.00000000000000000000 1296.00000000000000000000 0.00000000000000000000 7813.97998046875000000000 2000-01-01 00:00:00+00:00

Time-Series Analysis

Let’s explore the population density values in this data collection from 2000 to 2020. We can plot the dataset using the code below:

# Figure size: 10 is width, 5 is height
fig = plt.figure(figsize=(10,5))

# Use which_stat to choose any statistic from our DataFrame to visualize.
which_stat = "mean"

plt.plot(
    df["date"], # X-axis: sorted datetime
    df[which_stat], # Y-axis: maximum pop density
    color="red", # Line color
    linestyle="-", # Line style
    linewidth=2, # Line width
    label=f"{which_stat.capitalize()} {items[0].assets[asset_name].title}", # Legend label,
    marker='o' # Add circles to mark each data point
)

# Display legend
plt.legend()

# Insert label for the X-axis
plt.xlabel("Year")

# Insert label for the Y-axis
plt.ylabel("persons/km$^2$")

# Insert title for the plot
plt.title(f"{which_stat.capitalize()} {items[0].assets[asset_name].title} over {aoi_name}")


# Add data citation
plt.text(
    df["date"].iloc[0],           # X-coordinate of the text
    df[which_stat].min(),              # Y-coordinate of the text
    # Text to be displayed
    f"Source: {collection.title}",                  
    fontsize=9,                             # Font size
    horizontalalignment="right",             # Horizontal alignment
    verticalalignment="bottom",              # Vertical alignment
    color="blue",                            # Text color
)


# Plot the time series
plt.show()

Summary

In this notebook, we have successfully completed the following steps for the SEDAC Gridded World Population Density dataset:

  1. Install and import the necessary libraries
  2. Fetch the collection from STAC using the appropriate endpoints
  3. Count the number of existing granules within the collection
  4. Map and compare population density for two distinctive months/years
  5. Generate zonal statistics for the area of interest (AOI)
  6. Plot a time series of mean population density values for a specified region

If you have any questions regarding this user notebook, please contact us using the feedback form.

Back to top