GDAL with HDF5 and VRT

from osgeo import gdal
from pyproj.crs import CRS

gdal.UseExceptions()

def warp_resample():
    src = "vrt://earthaccess_data/20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc?a_ullr=-179.995,89.995,180.005,-89.995&sd_name=analysed_sst"
    output = ""
    output_format = "MEM"
    dstSRS = "EPSG:3857"

    width = height = 256
    te = [
        -20037508.342789244,
        -20037508.342789244,
        20037508.342789244,
        20037508.342789244,
    ]
    gt = [
        te[0],
        (te[2] - te[0]) / width,
        0,
        te[3],
        0,
        -(te[3] - te[1]) / height,
    ]
    output_crs = CRS(dstSRS).to_wkt()

    output_ds = gdal.GetDriverByName(output_format).Create(
        output, width, height, 1, gdal.GDT_Byte
    )
    output_ds.SetProjection(output_crs)
    output_ds.SetGeoTransform(gt)
    input_ds = gdal.Open(src)
    return gdal.Warp(output_ds, input_ds)

if __name__ == "__main__":
    warp_resample()