GDAL with NetCDF, VSIS3, and earthaccess

from osgeo import gdal
from pyproj.crs import CRS

gdal.UseExceptions()

def configure_auth():
    import earthaccess

    auth = earthaccess.login()
    s3_credentials = auth.get_s3_credentials("PODAAC")
    gdal.SetConfigOption("AWS_REGION", "us-west-2")
    gdal.SetConfigOption("AWS_SECRET_ACCESS_KEY", s3_credentials["secretAccessKey"])
    gdal.SetConfigOption("AWS_ACCESS_KEY_ID", s3_credentials["accessKeyId"])
    gdal.SetConfigOption("AWS_SESSION_TOKEN", s3_credentials["sessionToken"])

def warp_resample():
    bucket = "podaac-ops-cumulus-protected"
    input_uri = "MUR-JPL-L4-GLOB-v4.1/20020601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc"
    src = f"NETCDF:/vsis3/{bucket}/{input_uri}: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__":
    configure_auth()
    warp_resample()