Rasterio with NetCDF, VSIS3, and earthaccess

import rasterio as rio
import rioxarray as rx
from rasterio.session import AWSSession
from rasterio.warp import calculate_default_transform

def configure_auth():
    import boto3
    import earthaccess

    auth = earthaccess.login()
    s3_credentials = auth.get_s3_credentials("PODAAC")
    session = boto3.Session(
        aws_access_key_id=s3_credentials["accessKeyId"],
        aws_secret_access_key=s3_credentials["secretAccessKey"],
        aws_session_token=s3_credentials["sessionToken"],
        region_name="us-west-2",
    )
    rio_env = rio.Env(
        AWSSession(session),
    )
    rio_env.__enter__()
    return rio_env

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"
    dstSRS = "EPSG:3857"
    srcSRS = "EPSG:4326"
    width = height = 256
    te = [
        -20037508.342789244,
        -20037508.342789244,
        20037508.342789244,
        20037508.342789244,
    ]

    da = rx.open_rasterio(src)  # , mask_and_scale=True)
    da = da.rio.write_crs(srcSRS)
    da = da.rio.clip_box(
        *te,
        crs=dstSRS,
    )
    dst_transform, w, h = calculate_default_transform(
        srcSRS,
        dstSRS,
        da.rio.width,
        da.rio.height,
        *da.rio.bounds(),
        dst_width=width,
        dst_height=height,
    )
    return da.rio.reproject(dstSRS, shape=(h, w), transform=dst_transform)

if __name__ == "__main__":
    rio_env = configure_auth()
    warp_resample()
    rio_env.__exit__()