Reducción de la resolución de los datos NetCDF con GDAL

Question

Estoy descargando los datos 0.05 Deg NASA MODIS CMG en formato netCDF ( https://lpdaac.usgs.gov/products/vnp13c1v001/ ). Las coordenadas geográficas de los datos MODIS CMG son tales que la esquina superior izquierda del píxel superior izquierdo son -180,00 grados de longitud, 90,00 grados de latitud, y las coordenadas geográficas de la esquina inferior derecha del píxel inferior derecho son 180,00 grados de longitud, -90,00 grados de latitud. Hay 7200 x 3600 puntos de datos en el archivo de datos, y cada punto se sitúa en el centro de un píxel. Mi propósito es visualizar estos datos en Web Mercator utilizando OpenLayers, y me gustaría reducir en gran medida la resolución de los datos.

De la documentación de gdal_translate ( https://gdal.org/programs/gdal_translate.html ), parece que si utilizo el indicador -tr, puedo establecer la resolución de destino del archivo de salida, por ejemplo así: -tr 360 180 . Pero cuando corrí:

gdal_translate -of XYZ -tr 360 180 HDF5:"VNP13C1.A2020065.001.2020081184151.h5"://HDFEOS/GRIDS/NPP_Grid_16Day_VI_CMG/Data_Fields/CMG_0.05_Deg_16_days_EVI2 output.xyz

Recibí el mensaje de error:

The -tr option was used, but there's no geotransform or it is rotated.  This configuration is not supported.

¿Cómo podría reducir la resolución de los datos con gdal?

Sería estupendo que la reducción de resolución utilizara la interpolación para crear puntos de datos, en lugar de limitarse a filtrarlos.

Answer 1

Espero haber descargado datos similares a los tuyos.

Primero eche un vistazo a un subconjunto con gdalinfo. He aquí un informe reducido

gdalinfo HDF5:"VNP13C1.A2020065.001.2020081184151.h5"://HDFEOS/GRIDS/NPP_Grid_16Day_VI_CMG/Data_Fields/CMG_0.05_Deg_16_days_EVI2
Driver: HDF5Image/HDF5 Dataset
Files: VNP13C1.A2020065.001.2020081184151.h5
Size is 7200, 3600
Metadata:

  EastBoundingCoord=180
  LongName=VIIRS/NPP Vegetation Indices 16-Day L3 Global 0.05Deg CMG
  NorthBoundingCoord=90
  SouthBoundingCoord=-90
  WestBoundingCoord=-180
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 3600.0)
Upper Right ( 7200.0,    0.0)
Lower Right ( 7200.0, 3600.0)
Center      ( 3600.0, 1800.0)
Band 1 Block=7200x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-15000

Se puede ver que la extensión de los datos se almacena en los metadatos y los datos abarcan todo el mundo. Sin embargo, GDAL no ha podido encontrar el sistema de referencia de coordenadas e informa de las coordenadas de las esquinas en el espacio de píxeles. Pero como se conocen los límites, es posible crear una versión georreferenciada y submuestreada a partir de un subconjunto. Yo lo haría por pasos.

1) Convierta toda la imagen en un archivo TIFF temporal. Crear el archivo TIFF en mosaico y comprimido agiliza los pasos siguientes y ahorra espacio en disco. El sistema de referencia de coordenadas y los coordiantes de las esquinas superior izquierda e inferior derecha pueden adjuntarse al mismo:

gdal_translate -of GTiff -co tiled=yes -co compress=deflate -a_srs epsg:4326 -a_ullr -180 90 180 -90 HDF5:"VNP13C1.A2020065.001.2020081184151.h5"://HDFEOS/GRIDS/NPP_Grid_16Day_VI_CMG/Data_Fields/CMG_0.05_Deg_16_days_EVI2 viirs_temp.tif
Input file size is 7200, 3600
0...10...20...30...40...50...60...70...80...90...100 - done.

Compruebe el resultado visualmente

y con gdalinfo

gdalinfo viirs_temp.tif
Driver: GTiff/GeoTIFF
Files: viir_temp2.tif
Size is 7200, 3600
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.050000000000000,-0.050000000000000)
Metadata:

Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=256x256 Type=Int16, ColorInterp=Gray
  NoData Value=-15000

2) Crear un archivo XYZ submuestreado a partir del GeoTIFF. Si he entendido su descripción que desea obtener una salida que es más pequeño en tamaño, 360x180 píxeles. Entonces la opción correcta de gdal_translate es -outsize utilizado junto con -r para seleccionar el método de remuestreo. Más información en la documentación https://gdal.org/programs/gdal_translate.html .

gdal_translate -of XYZ -outsize 360 180 -r average viirs_temp.tif viirs.xyz
Input file size is 7200, 3600
0...10...20...30...40...50...60...70...80...90...100 - done.

gdalinfo viirs.xyz
Driver: XYZ/ASCII Gridded XYZ
Files: viirs.xyz
Size is 360, 180
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000)
Lower Left  (-180.0000000, -90.0000000)
Upper Right ( 180.0000000,  90.0000000)
Lower Right ( 180.0000000, -90.0000000)
Center      (   0.0000000,   0.0000000)
Band 1 Block=360x1 Type=Int16, ColorInterp=Undefined
  Min=-15000.000 Max=6707.000

Escribiste que apuntabas a Web Mercator. Puedes usar el programa gdalwarp para el warping y puedes usar el archivo tiff temporal como entrada. Lee la documentación sobre los métodos de remuestreo y otros parámetros. https://gdal.org/programs/gdalwarp.html .