La proyección polar provoca la pérdida de datos visualizados, el recorte impar cerca de los polos para MODIS (pero no DEM) en GEE

Question

Quiero crear visualizaciones (estáticas y timelapse) de MODIS, etc. con una proyección polar (como WGS 84 / EPSG Alaska Polar Stereographic, como la utilizada por Explorador ArcticDEM ) en el motor de Google Earth.

Cuando creo la miniatura animada de una colección de imágenes MODIS y elijo la CRS polar ( 'EPSG:5936' ), obtengo un gran agujero negro en medio de mi escena, al igual que este ejemplo de GDAL . Intenté recortar el conjunto de datos MODIS a mi ROI pero GEE parece comportarse de forma extraña si le das polígonos demasiado cercanos al polo norte. Por ejemplo, el geometry se representa tanto en la consola como en un clíper real como una fina línea que rodea el globo. Sin embargo, si defino una geometría como un rectángulo delgado en latitud, puedo obtener datos cercanos al Polo Norte para renderizar (por ejemplo, si sólo recorto Groenlandia puedo obtener toda Groenlandia). Se puede ver esto en acción si se aumenta la latitud del punto para que esté cada vez más cerca de 90; paradójicamente, se renderizarán las latitudes más altas cuanto más baja sea la latitud elegida para el centro del punto.

He observado que el conjunto de datos ArcticDEM, que viene con la proyección estereográfica, no tiene este problema. Así que intenté reproyectar las escenas MODIS antes de renderizarlas por si ese era el problema (sabiendo que se desaconseja la reproyección). Pero esto da como resultado una miniatura de la imagen rota cuando se imprime en la consola.

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);

¿Qué ocurre? ¿Cómo puedo definir un ROI que cubra el Círculo Polar Ártico y que realmente rinda (¿reproyecte correctamente?)?

Answer 1

Sospecho que esto tiene que ver con el hecho de que los datos MODIS no se ingieren con un CRS polar, lo que da lugar a transformaciones erróneas para los niveles bajos. capas de la pirámide . Puede ver en este script del Editor de Código que después de cambiar la proyección, las capas piramidales de alto nivel de zoom están completas, pero en los niveles de zoom bajos, las latitudes altas están enmascaradas. Al generar las miniaturas en CRS polar, GEE probablemente está dibujando a partir de capas de pirámide de bajo nivel de zoom.

También puede utilizar la API de GEE Python con cartopy para reproyectar las miniaturas de GEE Mercator a un variedad de proyecciones y luego animar los cuadros con ffmpeg. Aquí hay una lista para ejecutar Cuaderno Colab que lo demuestra; el siguiente es el mismo código, menos la configuración del entorno.

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