Polar Projection Data Loss – How to Fix Data Loss and Odd Clipping Near Poles for MODIS in GEE

Question

I want to create visualizations (static and timelapse) of MODIS etc. with a polar projection (such as WGS 84 / EPSG Alaska Polar Stereographic, as used by ArcticDEM Explorer) in Google Earth Engine.

When I create the animated thumbnail of a MODIS image collection and choose the polar CRS ('EPSG:5936'), I get a big black hole in the middle of my scene, just like this GDAL example. I tried to clip the MODIS dataset to my ROI but GEE seems to behave oddly if you give it polygons too near to the north pole. For example, the geometry object below renders both in the console and as an actual clipper as a thin line circling the globe. However, if I define a geometry as a latitudinally thin rectangle, I can get data close the North Pole to render (e.g. if I just clip Greenland I can get all of Greenland). You can see this in action if you increase the latitude of the point to be closer and closer to 90; paradoxically, higher latitudes will be rendered the lower the latitude that is chosen for the point center.

I noted that the ArcticDEM dataset, which comes with the stereographic projection, does not have this issue. So, I tried to reproject MODIS scenes before rendering them in case that was the issue (knowing that reprojection is advised against). But this results in a broken image thumbnail when printed to the console.

var geometry = ee.Geometry.Point([180, 90]).buffer(2500000);

// Display a clipped version of the mosaic.
Map.addLayer(geometry ,{}, "circle");

var arcticDEM = ee.Image('UMN/PGC/ArcticDEM/V3/2m_mosaic');

var collection = ee.ImageCollection("MODIS/006/MOD13Q1")
  .filterDate('2018-01-01', '2019-01-01')
  .select('NDVI')
  //fails if reprojected
  //.map(function(image) { return image.reproject('EPSG:5936', null, 250);    })
  .map(function(image) { return image.clip(geometry); }); 

// Visualization parameters.
var args = {
  crs: 'EPSG:5936',  // WGS 84 / EPSG Alaska Polar Stereographic
  //crs: 'EPSG:3857',  // Maps Mercator
  dimensions: '300',
  region: geometry,
  min: -2000,
  max: 8000,
  palette: 'black, blanchedalmond, #8FBC8F, #006400',
  framesPerSecond: 6,
};

var thumb = ui.Thumbnail({
  image: collection,//try arcticDEM, it works!
  params: args,
  style: {
    width: '300px'
  }
  });
print(thumb);

What's going on? How can I define an ROI that covers the Arctic Circle and actually renders (correctly re-projects?)?

Answer 1

I suspect this has something to do with MODIS data not being ingested with a polar CRS resulting in wonky transformations for low level pyramid layers. You can see in this Code Editor script that after changing the projection, high zoom level pyramid layers are complete, but at low zoom levels, high latitudes are masked. When generating the thumbnails in polar CRS, GEE is probably drawing from low zoom level pyramid layers.

You can alternatively use the GEE Python API with cartopy to reproject GEE Mercator thumbnails to a variety of projections, and then animate the frames with ffmpeg. Here is a ready-to-run Colab notebook that demonstrates it; the following is the same code, minus environment setup.

# Import packages
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import ee
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import numpy as np
import os
import requests

# Initialize EE
ee.Initialize()

# Make visualization collection in Earth Engine
col = (ee.ImageCollection('MODIS/006/MOD13Q1')
           .filterDate('2018-01-01', '2019-01-01'))

vis_params = {
    'bands': 'NDVI',
    'min': -2000,
    'max': 8000,
    'palette': ['black', 'blanchedalmond', '8FBC8F', '006400'],
    'forceRgbOutput': True
}

aoi = {
    'xmin': -180,
    'ymin': 60,
    'xmax': 180,
    'ymax': 90
}

def vis_rgb(img):
    return (img.visualize(**vis_params)
                 .unmask(0)
                 .set({'date': img.date().format('YYYY-MM-dd')}))

col_vis = (col.map(vis_rgb))

#Get thumbnails of the visualization collection downloaded locally
thumb_params = {
  'dimensions': 5000,
  'region': ee.Geometry.BBox(aoi['xmin'], aoi['ymin'], aoi['xmax'], aoi['ymax']),
  'crs': 'EPSG:4326',
  'format': 'PNG'
}

dates = col_vis.aggregate_array('date').getInfo()

frames = []
for count, date in enumerate(dates):
    img_vis = col_vis.filter(ee.Filter.eq('date', date)).first()
    img_url = img_vis.getThumbURL(thumb_params)
    img_name = str(count).zfill(3) + '.png'
    frames.append(img_name)
    print(img_name)
    response = requests.get(img_url)
    with open(img_name, 'wb') as fd:
        fd.write(response.content)

# Make a directory to hold altered images
alt_dir = 'alt'
if not os.path.exists(alt_dir):
    os.mkdir(alt_dir)

# Use cartopy to change the projection of the downloaded images to North Polar Stereo, set the extent of the output plot, and clip by a circle.
def make_plot(frame):
    fig = plt.figure(figsize=[8, 8])
    ax = plt.axes(projection=ccrs.NorthPolarStereo())
    ax.spines['geo'].set_edgecolor('none')
    ax.set_extent([aoi['xmin'], aoi['xmax'], aoi['ymin'], aoi['ymax']], ccrs.PlateCarree())
    ax.add_feature(cfeature.COASTLINE, edgecolor='grey')
    ax.gridlines()

    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpath.Path(verts * radius + center)

    ax.set_boundary(circle, transform=ax.transAxes)
    img = plt.imread(frame)
    ax.imshow(img, extent=(aoi['xmin'], aoi['xmax'], aoi['ymin'], aoi['ymax']), origin='upper', transform=ccrs.PlateCarree())
    fig.savefig(os.path.join(alt_dir, frame), bbox_inches='tight', pad_inches=0.5, facecolor='grey', edgecolor='none')
    plt.close()

for frame in frames:
    print(frame)
    make_plot(frame)

# Make an animation from the frames
outfile = 'animation'
cmd = 'ffmpeg -framerate 10 -i alt/%03d.png -movflags faststart -pix_fmt yuv420p -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2"'

cmd_mp4 = f'{cmd} {outfile}.mp4' 
os.system(cmd_mp4)

cmd_gif = f'{cmd} {outfile}.gif' 
os.system(cmd_gif)

# Download the MP4 and GIF animations
from google.colab import files
files.download(f'{outfile}.mp4')
files.download(f'{outfile}.gif')