Reducing resolution of NetCDF data with GDAL

Question

I am downloading 0.05 Deg NASA MODIS CMG data in netCDF form (https://lpdaac.usgs.gov/products/vnp13c1v001/). The geographic coordinates of MODIS CMG data are such that the upper-left corner of the upper-left pixel are -180.00 degrees longitude, 90.00 degrees latitude, and the geographic coordinates of the lower-right corner of the lower right pixel are 180.00 degrees longitude, -90.00 degrees latitude. There are 7200 x 3600 data points in the data file, with each point landing in the center of a pixel. My purpose is to visualize this data in Web Mercator using OpenLayers, and I would like to greatly reduce the resolution of the data.

From the documentation of gdal_translate (https://gdal.org/programs/gdal_translate.html), it seems that if I use the flag -tr, I can set the target resolution of the output file, for example like so: -tr 360 180. But when I ran:

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

I got the error message:

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

How might I reduce the resolution of the data with gdal?

It would be great if the resolution reduction used interpolation to create data points as opposed to simply filtering out data points.

Answer 1

Hopefully I downloaded similar data than what you have.

First have a look at one subset with gdalinfo. Here is a reduced report

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

You can see that the extent of data is stored into the metadata and data covers the whole world. However, GDAL was not able to find the coordinate reference system and it is reporting the corner coordinates in the pixel space. But because the bounds are known it is possible to create a georeferenced and subsampled version from a subset. I would do it in steps.

1) Convert the whole image into a temporary TIFF file. Creating the TIFF file as tiled and compressed makes next steps faster and saves disk space. The coordinate reference system and the coordiantes of the upper left and lower right corners can be attached by the same:

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.

Check the result visually

and with 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) Create a subsampled XYZ file from the GeoTIFF. If I understood your description you want to get an output that is smaller in size, 360x180 pixels. Then the right gdal_translate option is -outsize used together with -r for selecting the resampling method. Read details from the documentation 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

You wrote that you were aiming at Web Mercator. You can use gdalwarp program for warping and you can use the temporaty tiff file as an input. Read the documentation about the resampling methods and other parameters https://gdal.org/programs/gdalwarp.html.