Adding Cloud-Optimized GeoTIFFs to the MAAP Dashboard

The following notebook steps through how to add a dataset to the MAAP Dashboard.

Note, there are 2 scenarios:

  1. Adding a single Cloud-Optimized GeoTIFF (COG) and,
  2. Adding many distinct COGs as a “mosaic” with mosaicJSON.

High-level, the steps are:

  1. Inspect your COG(s) to understand the best rescale and colormap name parameters. Optionally create a mosaic.
  2. Define a colormap. Colormaps provide mappings of data values to RGB values.
  3. Create a PR to the datasets repo to add or update your dataset.

The MAAP dashboard has 3 environments:

  1. Developer-in-test (DIT): https://biomass.dit.maap-project.org
  2. Staging: https://uat.earthdata.nasa.gov/maap-biomass
  3. Production: https://earthdata.nasa.gov/maap-biomass

These instructions will guide you towards adding your dataset to biomass.dit.maap-project.org. The MAAP Dashboard team will “promote” changes to staging and production periodically (release schedule forthcoming).

Setup

[ ]:

%%capture
!pip install rasterio rio-cogeo supermercado

[1]:

# import the MAAP package
import ipycmc
from maap.maap import MAAP

import glob
import json
import os
import matplotlib
import urllib

# create MAAP class
maap = MAAP(maap_host='api.ops.maap-project.org')

[2]:

project_dir = "/projects/shared-buckets/<your_name>/<project_dir>"
# e.g.
project_dir = "/projects/shared-buckets/alexdevseed/landsat8/viz/"

[3]:

# Search for files to include, use recursive if nested folders (common in DPS output)
files = glob.glob(os.path.join(project_dir, "*.tif"), recursive=False)
files = [os.path.basename(f) for f in files]
#files

[4]:

my_tif = files[1]

[5]:

%%bash -s "$project_dir" "$my_tif"
rio cogeo validate $1/$2

/projects/shared-buckets/alexdevseed/landsat8/viz/Copernicus_30439_covars_cog_topo_stack.tif is a valid cloud optimized GeoTIFF

[6]:

%%bash -s "$project_dir" "$my_tif"
gdalinfo $1/$2 -stats

Driver: GTiff/GeoTIFF
Files: /projects/shared-buckets/alexdevseed/landsat8/viz//Copernicus_30439_covars_cog_topo_stack.tif
       /projects/shared-buckets/alexdevseed/landsat8/viz//Copernicus_30439_covars_cog_topo_stack.tif.aux.xml
Size is 1000, 1000
Coordinate System is:
PROJCRS["unknown",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101004,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["Albers Equal Area",
        METHOD["Albers Equal Area",
            ID["EPSG",9822]],
        PARAMETER["Latitude of false origin",40,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",180,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",50,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",70,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (698522.000000000000000,3153304.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  OVR_RESAMPLING_ALG=NEAREST
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  LAYOUT=COG
Corner Coordinates:
Upper Left  (  698522.000, 3153304.000) (163d20'28.79"W, 67d29'44.89"N)
Lower Left  (  698522.000, 3123304.000) (163d30'48.19"W, 67d14'14.17"N)
Upper Right (  728522.000, 3153304.000) (162d39'23.92"W, 67d25'44.51"N)
Lower Right (  728522.000, 3123304.000) (162d50' 6.57"W, 67d10'16.36"N)
Center      (  713522.000, 3138304.000) (163d 5'11.87"W, 67d20' 1.17"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = elevation
  Min=-1.830 Max=557.863
  Minimum=-1.830, Maximum=557.863, Mean=74.140, StdDev=84.128
  Overviews: 500x500, 250x250
  Metadata:
    STATISTICS_MAXIMUM=557.86291503906
    STATISTICS_MEAN=74.140305286651
    STATISTICS_MINIMUM=-1.8297636508942
    STATISTICS_STDDEV=84.128259469301
    STATISTICS_VALID_PERCENT=100
Band 2 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = slope
  Min=0.000 Max=40.566
  Minimum=0.000, Maximum=40.566, Mean=2.972, StdDev=4.108
  Overviews: 500x500, 250x250
  Metadata:
    STATISTICS_MAXIMUM=40.566223144531
    STATISTICS_MEAN=2.9720626801483
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=4.1079777350157
    STATISTICS_VALID_PERCENT=100
Band 3 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = tsri
  Min=0.000 Max=1.000
  Minimum=0.000, Maximum=1.000, Mean=0.466, StdDev=0.333
  Overviews: 500x500, 250x250
  Metadata:
    STATISTICS_MAXIMUM=0.99999368190765
    STATISTICS_MEAN=0.46556238474922
    STATISTICS_MINIMUM=4.3117461245856e-06
    STATISTICS_STDDEV=0.33263731194436
    STATISTICS_VALID_PERCENT=100
Band 4 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = tpi
  Min=-4.774 Max=7.687
  Minimum=-4.774, Maximum=7.687, Mean=-0.000, StdDev=0.349
  Overviews: 500x500, 250x250
  Metadata:
    STATISTICS_MAXIMUM=7.6866292953491
    STATISTICS_MEAN=-4.5027400747145e-05
    STATISTICS_MINIMUM=-4.7744207382202
    STATISTICS_STDDEV=0.34907900787421
    STATISTICS_VALID_PERCENT=100
Band 5 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = slopemask
  Min=0.000 Max=1.000
  Minimum=0.000, Maximum=1.000, Mean=0.982, StdDev=0.128
  Overviews: 500x500, 250x250
  Metadata:
    STATISTICS_MAXIMUM=1
    STATISTICS_MEAN=0.98238805740667
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0.12799686533701
    STATISTICS_VALID_PERCENT=100

Create some parameters for the dynamic tiler.

These parameters will be pased to titiler_url for visualization.

Note the values below: band_min, band_max and colormap_name are set as the current defaults for biomass in Mg/hectare for the dashboard. For other datasets, users should inspect the output of the gdalinfo for band_min and band_max and modify the colormap_name as makes sense for the current dataset. This notebook includes section on what colormaps are available and how to configure different types of colormaps and legends.

[7]:

titiler_url = f"https://titiler.maap-project.org"
band_min = 0
band_max = 400
rescale = f"{band_min},{band_max}"
bidx = 1
colormap_name = 'gist_earth_r'

Helper functions

Step 1, Scenario 1: Adding a single COG

1. Upload the file

Only use the following steps if you only have one COG to share to the dashboard. If you want to create a mosaic from multiple COGs, skip to Step 1 Option 2.

Upload the file to S3 and make note of the location.

[8]:

def local_to_s3(url):
    ''' A Function to convert local paths to s3 urls'''
    return url.replace("/projects/shared-buckets", "s3://maap-ops-workspace/shared")

location = local_to_s3(f"{project_dir}{my_tif}")
location

[8]:

's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30439_covars_cog_topo_stack.tif'

2. Test the COG with the titiler and ipyCMC

[9]:

## Build a WMTS call
wmts_url = f"{titiler_url}/cog/WMTSCapabilities.xml"
params = {
    "tile_format": "png",
    "tile_scale": "1",
    "pixel_selection": "first",
    "TileMatrixSetId": "WebMercatorQuad",
    "url": location,
    "bidx": bidx, # Select which band to use
    "resampling_method":"nearest",
    "rescale": rescale,
    "return_mask": "true",
    "colormap_name": colormap_name # Any colormap from matplotlib will work
}

wmts_call = "?".join([wmts_url, urllib.parse.urlencode(params)])

# Note Jupyter bug add amp; incorrectly when printing the url
wmts_call

[9]:

'https://titiler.maap-project.org/cog/WMTSCapabilities.xml?tile_format=png&tile_scale=1&pixel_selection=first&TileMatrixSetId=WebMercatorQuad&url=s3%3A%2F%2Fmaap-ops-workspace%2Fshared%2Falexdevseed%2Flandsat8%2Fviz%2FCopernicus_30439_covars_cog_topo_stack.tif&bidx=1&resampling_method=nearest&rescale=0%2C400&return_mask=true&colormap_name=gist_earth_r'

[10]:

w = ipycmc.MapCMC()
w

[11]:

w.load_layer_config(wmts_call, "wmts/xml")

Step 1, Scenario 2: Adding data from multiple COGs by creating a mosaic

Many datasets are comprised of many tiles distributed spatially over the globe. In order to visualize them all together, we can use mosaicJSON to create a mosaic for the dynamic tiler API. The dynamic tiler API knows how to read this mosaicJSON and select which tiles to render based on the current zoom, x and y coordinates across spatially distinct COGs.

[12]:

%%capture
!pip install cogeo-mosaic

[13]:

import glob
import os
import urllib

from cogeo_mosaic.mosaic import MosaicJSON
from cogeo_mosaic.backends import MosaicBackend

[14]:

# Skipping this step since we don't want to upload these large files more than once!
full_path_files = [f'{project_dir}{file}' for file in files]

[15]:

def local_to_s3(url):
    ''' A Function to convert local paths to s3 urls'''
    return url.replace("/projects/shared-buckets", "s3://maap-ops-workspace/shared")

tiles = [local_to_s3(file) for file in full_path_files]
print(tiles)

['s3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30438_covars_cog_topo_stack.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30439_covars_cog_topo_stack.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30440_covars_cog_topo_stack.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30542_covars_cog_topo_stack.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30543_covars_cog_topo_stack.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30821_covars_cog_topo_stack.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30822_covars_cog_topo_stack.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Copernicus_30823_covars_cog_topo_stack.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Landsat8_30542_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Landsat8_30543_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Landsat8_30822_comp_cog_2015-2020_dps.tif', 's3://maap-ops-workspace/shared/alexdevseed/landsat8/viz/Landsat8_30823_comp_cog_2015-2020_dps.tif']

[16]:

mosaicdata = MosaicJSON.from_urls(tiles, minzoom=1, maxzoom=16)

Using mosaicJSON with titiler

There are 2 options for using mosaicJSON with titiler:

  1. (Preferred) Post mosaicJSON to titiler mosaics endpoint and use the mosaicjson/mosaics endpoint for dynamic tiling.
  2. Upload mosaicJSON to S3 and pass the S3 url to the titiler mosaicjson/tiles endpoint.

Post mosaicJSON to titiler

[17]:

import requests
from pprint import pprint

r = requests.post(
    url=f"{titiler_url}/mosaics",
    headers={
        "Content-Type": "application/vnd.titiler.mosaicjson+json",
    },
    json=mosaicdata.dict(exclude_none=True)).json()

pprint(r)

{'id': 'edf63f99-84dc-447c-afec-90b7109feab4',
 'links': [{'href': 'https://titiler.maap-project.org/mosaics/edf63f99-84dc-447c-afec-90b7109feab4',
            'rel': 'self',
            'title': 'Self',
            'type': 'application/json'},
           {'href': 'https://titiler.maap-project.org/mosaics/edf63f99-84dc-447c-afec-90b7109feab4/mosaicjson',
            'rel': 'mosaicjson',
            'title': 'MosaicJSON',
            'type': 'application/json'},
           {'href': 'https://titiler.maap-project.org/mosaics/edf63f99-84dc-447c-afec-90b7109feab4/tilejson.json',
            'rel': 'tilejson',
            'title': 'TileJSON',
            'type': 'application/json'},
           {'href': 'https://titiler.maap-project.org/mosaics/edf63f99-84dc-447c-afec-90b7109feab4/tiles/{z}/{x}/{y}',
            'rel': 'tiles',
            'title': 'Tiles',
            'type': 'application/json'},
           {'href': 'https://titiler.maap-project.org/mosaics/edf63f99-84dc-447c-afec-90b7109feab4/WMTSCapabilities.xml',
            'rel': 'wmts',
            'title': 'WMTS',
            'type': 'application/json'}]}

[18]:

mosaic_link = list(filter(lambda x: x['rel'] == 'tiles', r['links']))[0]['href']
wmts_mosaic_link = list(filter(lambda x: x['rel'] == 'wmts', r['links']))[0]['href']

Test your mosaic

[19]:

params = {
    "tile_format": "png",
    "bidx": bidx, # Select which band to use
    "resampling_method":"nearest",
    "rescale": rescale,
    "return_mask": "true",
    "colormap_name": colormap_name # Any colormap from matplotlib will work
}

wmts_call = "?".join([wmts_mosaic_link, urllib.parse.urlencode(params)])

w.load_layer_config(wmts_call, "wmts/xml")

By default, the image will be displayed in greyscale if no colormap_name parameter is passed to the titiler API. Guidance below is provided to help determine what a valid colormap_name might be and how to create a legend for the dashboard.

Dashboard ColorRamps & Legends

When using the dashboard, there 2 components for implementing a color scheme for your map. There is the map render and there is the legend.

Titiler used for Cloud Optimized Geotiff (COG) rendering accepts any color scheme from the python matplotlib library, and custom color formulas.

Available colormap_name values for titiler: above, accent, accent_r, afmhot, afmhot_r, autumn, autumn_r, binary, binary_r, blues, blues_r, bone, bone_r, brbg, brbg_r, brg, brg_r, bugn, bugn_r, bupu, bupu_r, bwr, bwr_r, cfastie, cividis, cividis_r, cmrmap, cmrmap_r, cool, cool_r, coolwarm, coolwarm_r, copper, copper_r, cubehelix, cubehelix_r, dark2, dark2_r, flag, flag_r, gist_earth, gist_earth_r, gist_gray, gist_gray_r, gist_heat, gist_heat_r, gist_ncar, gist_ncar_r, gist_rainbow, gist_rainbow_r, gist_stern, gist_stern_r, gist_yarg, gist_yarg_r, gnbu, gnbu_r, gnuplot, gnuplot2, gnuplot2_r, gnuplot_r, gray, gray_r, greens, greens_r, greys, greys_r, hot, hot_r, hsv, hsv_r, inferno, inferno_r, jet, jet_r, magma, magma_r, nipy_spectral, nipy_spectral_r, ocean, ocean_r, oranges, oranges_r, orrd, orrd_r, paired, paired_r, pastel1, pastel1_r, pastel2, pastel2_r, pink, pink_r, piyg, piyg_r, plasma, plasma_r, prgn, prgn_r, prism, prism_r, pubu, pubu_r, pubugn, pubugn_r, puor, puor_r, purd, purd_r, purples, purples_r, rainbow, rainbow_r, rdbu, rdbu_r, rdgy, rdgy_r, rdpu, rdpu_r, rdylbu, rdylbu_r, rdylgn, rdylgn_r, reds, reds_r, rplumbo, schwarzwald, seismic, seismic_r, set1, set1_r, set2, set2_r, set3, set3_r, spectral, spectral_r, spring, spring_r, summer, summer_r, tab10, tab10_r, tab20, tab20_r, tab20b, tab20b_r, tab20c, tab20c_r, terrain, terrain_r, twilight, twilight_r, twilight_shifted, twilight_shifted_r, viridis, viridis_r, winter, winter_r, wistia, wistia_r, ylgn, ylgn_r, ylgnbu, ylgnbu_r, ylorbr, ylorbr_r, ylorrd, ylorrd_r

Example 1: Class based known colors

In this example, the raster represents classes of forest with 11 possible values. There are specific colors selected to correspond to each class. We combine the list of colors and the list of classes and format them for the legend parameter the dashboard needs.

https://github.com/MAAP-Project/dashboard-datasets-maap/blob/main/datasets/taiga-forest-classification.json

[ ]:

colors = [
    '#5255A3','#1796A3','#FDBF6F','#FF7F00', '#FFFFBF','#D9EF8B','#91CF60','#1A9850', '#C4C4C4','#FF0000', '#0000FF'
]
labels = [
    'Sparse & Uniform',
    'Sparse & Diffuse-gradual',
    'Sparse & Diffuse-rapid',
    'Sparse & Abrupt ',
    'Open & Uniform ',
    'Open & Diffuse-gradual',
    'Open & Diffuse-rapid',
    'Open & Abrupt',
    'Intermediate & Closed',
    'Non-forest edge (dry)',
    'Non-forest edge (wet)'
]

legend = [dict(color=colors[i],label=labels[i]) for i in range(0, len(colors))]
print(json.dumps(legend, indent=2))

# Copy and Paste the output below to your dashboard config.

Example 2: Discrete ColorRamp

In this example, the range of values is known, but the color scale has many non-sequential colors. Starting with the premade color list, we create a continuous color ramp that uses the known colors as stops points. Arbitrarly 12 breaks looked decent in the dashboard legend so we split it into 12 discrete colors. Then combine the list of values and colors into the correct json syntax.

https://github.com/MAAP-Project/dashboard-datasets-maap/blob/main/datasets/ATL08.json

[ ]:

forest_ht = matplotlib.colors.LinearSegmentedColormap.from_list('forest_ht', ['#636363','#FC8D59','#FEE08B','#FFFFBF','#D9EF8B','#91CF60','#1A9850','#005A32'], 12)
cols = [matplotlib.colors.to_hex(forest_ht(i)) for i in range(forest_ht.N)]

cats = range(0,25, (25//len(cols)))
legend = [[cats[i],cols[i]] for i in range(0, len(cols))]
text = (json.dumps(legend, separators=(',', ': ') ))

print(text.replace('],[','],\n['))

# Copy and Paste the output below to your dashboard config.

Example 3: Continuous ColorRamp

In this example, we are using a built in ColorRamp from matplotlib. So we just need to extract enough colors to fill the legend adequately, and convert the colors to hex codes.

https://github.com/MAAP-Project/dashboard-datasets-maap/blob/main/datasets/topo.json

[ ]:

cmap_name = 'gist_earth_r'
cmap = matplotlib.cm.get_cmap(cmap_name, 12)
cols = [matplotlib.colors.to_hex(cmap(i)) for i in range(cmap.N)]
print(cols)

# Copy and Paste the output below to your dashboard config.

Step 3: Create and submit your dashboard dataset json

[ ]:

# This example is for a continuous color ramps
dataset_id = "paraguay-estimated-biomass"
dataset_name = "Estimated Biomass in Paraguay"
dataset_type = "raster"
legend_type = "gradient-adjustable"
info = "Estimated biomass within 6km grids."
stops = cols
parameters = f"colormap_name={cmap_name}&rescale={band_min},{band_max}&bidx={bidx}"

[ ]:

# Single COG
tiles_link = f"{{titiler_server_url}}/cog/tiles/{{z}}/{{x}}/{{y}}.png?url={location}&{parameters}"

[ ]:

tiles_link

[ ]:

dataset_dict = {
    "id": dataset_id,
    "name": dataset_name,
    "type": dataset_type,
    "swatch": {
      "color": "#6976d7",
      "name": "Moody Blue"
    },
    "source": {
        "type": dataset_type,
        "tiles": [ tiles_link ]
    },
    "legend": {
      "type": legend_type,
       "min": band_min,
       "max": band_max,
      "stops": stops
    },
    "info": info
}
print(json.dumps(dataset_dict, indent=4))

Create a PR to the datasets repo

git clone git@github.com:MAAP-Project/biomass-dashboard-datasets.git
cd biomass-dashboard-datasets
git checkout -b ab/add-dataset
# select and copy json above
echo <copied_json> >> datasets/paraguay-estimated-biomass.json

Add to the datasets list in config.yml

In config.yml:

DATASETS:
- paraguay-estimated-biomass.json

Add to the product or country pilot

In country_pilots/paraguay/country_pilot.json:

{
    "id": "paraguay",
    "label": "Paraguay",
    //...
    "datasets": [
        {
            "id": "paraguay-forest-mask"
        },
        {
            "id": "paraguay-tree-cover"
        },
        {
            "id": "paraguay-estimated-biomass"
        }
    ]
}

Add content to the summary

There should be a summary.html file corresponding to the product or country pilot you are working on, for example: country_pilots/paraguay/summary.html. Add or modify content in that file as appropriate.

Submit a PR

Once you have added the dataset json file and summary content, submit a PR to https://github.com/MAAP-Project/biomass-dashboard-datasets. A member of the data team will review the PR and when it is merged your content will appear in biomass.dit.maap-project.org.