/******************************************************************************
 *
 * Project:  netCDF read/write Driver
 * Purpose:  Sentinel 3 SRAL/MWR Level 2 products
 * Author:   Even Rouault <even.rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2019, Even Rouault <even.rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

// Example product:
// https://scihub.copernicus.eu/dhus/odata/v1/Products('65b615b0-0db9-4ced-8020-eb17818f0c26')/$value
// Specification:
// https://sentinel.esa.int/documents/247904/2753172/Sentinel-3-Product-Data-Format-Specification-Level-2-Land

#include "netcdfdataset.h"

#include <limits>

/************************************************************************/
/*                      Sentinel3_SRAL_MWR_Layer                        */
/************************************************************************/

class Sentinel3_SRAL_MWR_Layer final : public OGRLayer
{
    OGRFeatureDefn *m_poFDefn = nullptr;
    int m_cdfid;
    size_t m_nCurIdx = 0;
    size_t m_nFeatureCount = 0;
    CPLStringList m_aosMetadata{};

    struct VariableInfo
    {
        int varid;
        nc_type nctype;
        double scale;
        double offset;
        double nodata;
    };

    std::vector<VariableInfo> m_asVarInfo{};
    int m_iLongitude = -1;
    int m_iLatitude = -1;
    double m_dfLongScale = 1.0;
    double m_dfLongOffset = 0.0;
    double m_dfLatScale = 1.0;
    double m_dfLatOffset = 0.0;

    OGRFeature *TranslateFeature(size_t nIndex);
    OGRFeature *GetNextRawFeature();

  public:
    Sentinel3_SRAL_MWR_Layer(const std::string &name, int cdfid, int dimid);
    ~Sentinel3_SRAL_MWR_Layer();

    OGRFeatureDefn *GetLayerDefn() override
    {
        return m_poFDefn;
    }

    void ResetReading() override;
    OGRFeature *GetNextFeature() override;
    OGRFeature *GetFeature(GIntBig nFID) override;
    GIntBig GetFeatureCount(int bForce) override;
    int TestCapability(const char *pszCap) override;
    char **GetMetadata(const char *pszDomain) override;
    const char *GetMetadataItem(const char *pszKey,
                                const char *pszDomain) override;
};

/************************************************************************/
/*                      Sentinel3_SRAL_MWR_Layer()                      */
/************************************************************************/

Sentinel3_SRAL_MWR_Layer::Sentinel3_SRAL_MWR_Layer(const std::string &name,
                                                   int cdfid, int dimid)
    : m_cdfid(cdfid)
{
    m_poFDefn = new OGRFeatureDefn(name.c_str());
    m_poFDefn->SetGeomType(wkbPoint);
    m_poFDefn->Reference();
    SetDescription(name.c_str());

    nc_inq_dimlen(cdfid, dimid, &m_nFeatureCount);

    int nVars = 0;
    NCDF_ERR(nc_inq(cdfid, nullptr, &nVars, nullptr, nullptr));
    for (int iVar = 0; iVar < nVars; iVar++)
    {
        int nVarDims = 0;
        NCDF_ERR(nc_inq_varndims(cdfid, iVar, &nVarDims));
        if (nVarDims != 1)
            continue;

        int vardimid = -1;
        NCDF_ERR(nc_inq_vardimid(cdfid, iVar, &vardimid));
        if (vardimid != dimid)
            continue;

        char szVarName[NC_MAX_NAME + 1] = {};
        NCDF_ERR(nc_inq_varname(cdfid, iVar, szVarName));

        nc_type vartype = NC_NAT;
        nc_inq_vartype(cdfid, iVar, &vartype);

        int nbAttr = 0;
        NCDF_ERR(nc_inq_varnatts(cdfid, iVar, &nbAttr));
        std::string scaleFactor;
        std::string offset;
        std::string fillValue;
        CPLStringList aosMetadata;
        for (int iAttr = 0; iAttr < nbAttr; iAttr++)
        {
            char szAttrName[NC_MAX_NAME + 1];
            szAttrName[0] = 0;
            NCDF_ERR(nc_inq_attname(cdfid, iVar, iAttr, szAttrName));
            char *pszMetaTemp = nullptr;
            if (NCDFGetAttr(cdfid, iVar, szAttrName, &pszMetaTemp) == CE_None &&
                pszMetaTemp)
            {
                if (EQUAL(szAttrName, "scale_factor"))
                {
                    scaleFactor = pszMetaTemp;
                }
                else if (EQUAL(szAttrName, "add_offset"))
                {
                    offset = pszMetaTemp;
                }
                else if (EQUAL(szAttrName, "_FillValue"))
                {
                    fillValue = pszMetaTemp;
                }
                else if (!EQUAL(szAttrName, "coordinates"))
                {
                    aosMetadata.SetNameValue(szAttrName, pszMetaTemp);
                }
            }
            CPLFree(pszMetaTemp);
        }

        const char *pszStandardName =
            aosMetadata.FetchNameValue("standard_name");
        if (pszStandardName)
        {
            if (EQUAL(pszStandardName, "longitude") && vartype == NC_INT)
            {
                m_iLongitude = iVar;
                if (!scaleFactor.empty())
                    m_dfLongScale = CPLAtof(scaleFactor.c_str());
                if (!offset.empty())
                    m_dfLongOffset = CPLAtof(offset.c_str());
                continue;
            }
            if (EQUAL(pszStandardName, "latitude") && vartype == NC_INT)
            {
                m_iLatitude = iVar;
                if (!scaleFactor.empty())
                    m_dfLatScale = CPLAtof(scaleFactor.c_str());
                if (!offset.empty())
                    m_dfLatOffset = CPLAtof(offset.c_str());
                continue;
            }
        }

        for (int i = 0; i < aosMetadata.size(); i++)
            m_aosMetadata.AddString(
                (std::string(szVarName) + '_' + aosMetadata[i]).c_str());

        OGRFieldType eType = OFTReal;
        if (!scaleFactor.empty())
        {
            // do nothing
        }
        else if (!offset.empty())
        {
            // do nothing
        }
        else if (vartype == NC_BYTE || vartype == NC_SHORT ||
                 vartype == NC_INT || vartype == NC_USHORT ||
                 vartype == NC_UINT)
        {
            eType = OFTInteger;
        }
        OGRFieldDefn oField(szVarName, eType);
        m_poFDefn->AddFieldDefn(&oField);
        VariableInfo varInfo;
        varInfo.varid = iVar;
        varInfo.nctype = vartype;
        varInfo.scale =
            scaleFactor.empty() ? 1.0 : CPLAtof(scaleFactor.c_str());
        varInfo.offset = offset.empty() ? 0.0 : CPLAtof(offset.c_str());
        varInfo.nodata = fillValue.empty()
                             ? std::numeric_limits<double>::quiet_NaN()
                             : CPLAtof(fillValue.c_str());
        m_asVarInfo.emplace_back(varInfo);
    }
}

/************************************************************************/
/*                     ~Sentinel3_SRAL_MWR_Layer()                      */
/************************************************************************/

Sentinel3_SRAL_MWR_Layer::~Sentinel3_SRAL_MWR_Layer()
{
    m_poFDefn->Release();
}

/************************************************************************/
/*                           GetMetadata()                              */
/************************************************************************/

char **Sentinel3_SRAL_MWR_Layer::GetMetadata(const char *pszDomain)
{
    if (pszDomain == nullptr || EQUAL(pszDomain, ""))
        return m_aosMetadata.List();
    return OGRLayer::GetMetadata(pszDomain);
}

/************************************************************************/
/*                           GetMetadataItem()                          */
/************************************************************************/

const char *Sentinel3_SRAL_MWR_Layer::GetMetadataItem(const char *pszKey,
                                                      const char *pszDomain)
{
    if (pszDomain == nullptr || EQUAL(pszDomain, ""))
        return m_aosMetadata.FetchNameValue(pszKey);
    return OGRLayer::GetMetadataItem(pszKey, pszDomain);
}

/************************************************************************/
/*                           ResetReading()                             */
/************************************************************************/

void Sentinel3_SRAL_MWR_Layer::ResetReading()
{
    m_nCurIdx = 0;
}

/************************************************************************/
/*                          GetFeatureCount()                           */
/************************************************************************/

GIntBig Sentinel3_SRAL_MWR_Layer::GetFeatureCount(int bForce)
{
    if (m_poFilterGeom == nullptr && m_poAttrQuery == nullptr)
        return m_nFeatureCount;
    return OGRLayer::GetFeatureCount(bForce);
}

/************************************************************************/
/*                           TestCapability()                           */
/************************************************************************/

int Sentinel3_SRAL_MWR_Layer::TestCapability(const char *pszCap)
{
    if (EQUAL(pszCap, OLCFastFeatureCount))
        return m_poFilterGeom == nullptr && m_poAttrQuery == nullptr;
    if (EQUAL(pszCap, OLCRandomRead))
        return true;
    return false;
}

/************************************************************************/
/*                        TranslateFeature()                            */
/************************************************************************/

OGRFeature *Sentinel3_SRAL_MWR_Layer::TranslateFeature(size_t nIndex)
{
    OGRFeature *poFeat = new OGRFeature(m_poFDefn);
    poFeat->SetFID(nIndex + 1);
    if (m_iLongitude >= 0 && m_iLatitude >= 0)
    {
        int nLong = 0;
        int status = nc_get_var1_int(m_cdfid, m_iLongitude, &nIndex, &nLong);
        if (status == NC_NOERR)
        {
            int nLat = 0;
            status = nc_get_var1_int(m_cdfid, m_iLatitude, &nIndex, &nLat);
            if (status == NC_NOERR)
            {
                const double dfLong = nLong * m_dfLongScale + m_dfLongOffset;
                const double dfLat = nLat * m_dfLatScale + m_dfLatOffset;
                auto poGeom = new OGRPoint(dfLong, dfLat);
                auto poGeomField = m_poFDefn->GetGeomFieldDefn(0);
                poGeom->assignSpatialReference(poGeomField->GetSpatialRef());
                poFeat->SetGeometryDirectly(poGeom);
            }
        }
    }

    for (int i = 0; i < static_cast<int>(m_asVarInfo.size()); i++)
    {
        if (m_asVarInfo[i].nctype == NC_BYTE)
        {
            signed char nVal = 0;
            int status = nc_get_var1_schar(m_cdfid, m_asVarInfo[i].varid,
                                           &nIndex, &nVal);
            if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
            {
                poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
                                        m_asVarInfo[i].offset);
            }
        }
        else if (m_asVarInfo[i].nctype == NC_SHORT)
        {
            short nVal = 0;
            int status = nc_get_var1_short(m_cdfid, m_asVarInfo[i].varid,
                                           &nIndex, &nVal);
            if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
            {
                poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
                                        m_asVarInfo[i].offset);
            }
        }
        else if (m_asVarInfo[i].nctype == NC_USHORT)
        {
            unsigned short nVal = 0;
            int status = nc_get_var1_ushort(m_cdfid, m_asVarInfo[i].varid,
                                            &nIndex, &nVal);
            if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
            {
                poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
                                        m_asVarInfo[i].offset);
            }
        }
        else if (m_asVarInfo[i].nctype == NC_INT)
        {
            int nVal = 0;
            int status =
                nc_get_var1_int(m_cdfid, m_asVarInfo[i].varid, &nIndex, &nVal);
            if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
            {
                poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
                                        m_asVarInfo[i].offset);
            }
        }
        else if (m_asVarInfo[i].nctype == NC_UINT)
        {
            unsigned int nVal = 0;
            int status =
                nc_get_var1_uint(m_cdfid, m_asVarInfo[i].varid, &nIndex, &nVal);
            if (status == NC_NOERR && nVal != m_asVarInfo[i].nodata)
            {
                poFeat->SetField(i, nVal * m_asVarInfo[i].scale +
                                        m_asVarInfo[i].offset);
            }
        }
        else if (m_asVarInfo[i].nctype == NC_DOUBLE)
        {
            double dfVal = 0.0;
            int status = nc_get_var1_double(m_cdfid, m_asVarInfo[i].varid,
                                            &nIndex, &dfVal);
            if (status == NC_NOERR && dfVal != m_asVarInfo[i].nodata)
            {
                poFeat->SetField(i, dfVal * m_asVarInfo[i].scale +
                                        m_asVarInfo[i].offset);
            }
        }
        else
        {
            CPLDebug("netCDF", "Unhandled data type %d for %s",
                     m_asVarInfo[i].nctype,
                     m_poFDefn->GetFieldDefn(i)->GetNameRef());
        }
    }

    return poFeat;
}

/************************************************************************/
/*                        GetNextRawFeature()                           */
/************************************************************************/

OGRFeature *Sentinel3_SRAL_MWR_Layer::GetNextRawFeature()
{
    if (m_nCurIdx == m_nFeatureCount)
        return nullptr;
    OGRFeature *poFeat = TranslateFeature(m_nCurIdx);
    m_nCurIdx++;
    return poFeat;
}

/************************************************************************/
/*                           GetFeature()                               */
/************************************************************************/

OGRFeature *Sentinel3_SRAL_MWR_Layer::GetFeature(GIntBig nFID)
{
    if (nFID <= 0 || static_cast<size_t>(nFID) > m_nFeatureCount)
        return nullptr;
    return TranslateFeature(static_cast<size_t>(nFID - 1));
}

/************************************************************************/
/*                           GetNextFeature()                           */
/************************************************************************/

OGRFeature *Sentinel3_SRAL_MWR_Layer::GetNextFeature()
{
    while (true)
    {
        OGRFeature *poFeature = GetNextRawFeature();
        if (poFeature == nullptr)
            return nullptr;

        if ((m_poFilterGeom == nullptr ||
             FilterGeometry(poFeature->GetGeomFieldRef(m_iGeomFieldFilter))) &&
            (m_poAttrQuery == nullptr || m_poAttrQuery->Evaluate(poFeature)))
            return poFeature;

        delete poFeature;
    }
}

/************************************************************************/
/*                       ProcessSentinel3_SRAL_MWR()                    */
/************************************************************************/

void netCDFDataset::ProcessSentinel3_SRAL_MWR()
{
    int nDimCount = -1;
    int status = nc_inq_ndims(cdfid, &nDimCount);
    NCDF_ERR(status);
    if (status != NC_NOERR || nDimCount == 0 || nDimCount > 1000)
        return;
    std::vector<int> dimIds(nDimCount);
    int nDimCount2 = -1;
    status = nc_inq_dimids(cdfid, &nDimCount2, &dimIds[0], FALSE);
    NCDF_ERR(status);
    if (status != NC_NOERR)
        return;
    CPLAssert(nDimCount == nDimCount2);

    OGRSpatialReference *poSRS = nullptr;
    const char *pszSemiMajor =
        CSLFetchNameValue(papszMetadata, "NC_GLOBAL#semi_major_ellipsoid_axis");
    const char *pszFlattening =
        CSLFetchNameValue(papszMetadata, "NC_GLOBAL#ellipsoid_flattening");
    if (pszSemiMajor && EQUAL(pszSemiMajor, "6378137") && pszFlattening &&
        fabs(CPLAtof(pszFlattening) - 0.00335281066474748) < 1e-16)
    {
        int iAttr =
            CSLFindName(papszMetadata, "NC_GLOBAL#semi_major_ellipsoid_axis");
        if (iAttr >= 0)
            papszMetadata = CSLRemoveStrings(papszMetadata, iAttr, 1, nullptr);
        iAttr = CSLFindName(papszMetadata, "NC_GLOBAL#ellipsoid_flattening");
        if (iAttr >= 0)
            papszMetadata = CSLRemoveStrings(papszMetadata, iAttr, 1, nullptr);
        poSRS = new OGRSpatialReference();
        poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
        poSRS->importFromEPSG(4326);
    }

    for (int i = 0; i < nDimCount; i++)
    {
        char szDimName[NC_MAX_NAME + 1] = {};
        status = nc_inq_dimname(cdfid, dimIds[i], szDimName);
        NCDF_ERR(status);
        if (status != NC_NOERR)
            break;
        std::string name(CPLGetBasenameSafe(GetDescription()));
        name += '_';
        name += szDimName;
        std::shared_ptr<OGRLayer> poLayer(
            new Sentinel3_SRAL_MWR_Layer(name.c_str(), cdfid, dimIds[i]));
        auto poGeomField = poLayer->GetLayerDefn()->GetGeomFieldDefn(0);
        if (poGeomField)
            poGeomField->SetSpatialRef(poSRS);
        papoLayers.emplace_back(poLayer);
    }

    if (poSRS)
        poSRS->Release();
}