/******************************************************************************
 *
 * Project:     TerraSAR-X XML Product Support
 * Purpose:     Support for TerraSAR-X XML Metadata files
 * Author:      Philippe Vachon <philippe@cowpig.ca>
 * Description: This driver adds support for reading metadata and georef data
 *              associated with TerraSAR-X products.
 *
 ******************************************************************************
 * Copyright (c) 2007, Philippe Vachon <philippe@cowpig.ca>
 * Copyright (c) 2009-2012, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_minixml.h"
#include "gdal_frmts.h"
#include "gdal_pam.h"
#include "gdal_driver.h"
#include "gdal_drivermanager.h"
#include "gdal_openinfo.h"
#include "gdal_cpp_functions.h"
#include "ogr_spatialref.h"

#define MAX_GCPS 5000  // this should be more than enough ground control points

namespace gdal::TSX
{
enum ePolarization
{
    HH = 0,
    HV,
    VH,
    VV
};

enum eProductType
{
    eSSC = 0,
    eMGD,
    eEEC,
    eGEC,
    eUnknown
};
}  // namespace gdal::TSX

using namespace gdal::TSX;

/************************************************************************/
/* Helper Functions                                                     */
/************************************************************************/

/* GetFilePath: return a relative path to a file within an XML node.
 * Returns Null on failure
 */
static CPLString GetFilePath(CPLXMLNode *psXMLNode, const char **pszNodeType)
{
    const char *pszDirectory =
        CPLGetXMLValue(psXMLNode, "file.location.path", "");
    const char *pszFilename =
        CPLGetXMLValue(psXMLNode, "file.location.filename", "");
    *pszNodeType = CPLGetXMLValue(psXMLNode, "type", " ");

    if (pszDirectory == nullptr || pszFilename == nullptr)
    {
        return "";
    }

    return CPLString(pszDirectory) + '/' + pszFilename;
}

/************************************************************************/
/* ==================================================================== */
/*                                TSXDataset                                 */
/* ==================================================================== */
/************************************************************************/

class TSXDataset final : public GDALPamDataset
{
    int nGCPCount;
    GDAL_GCP *pasGCPList;

    OGRSpatialReference m_oGCPSRS{};

    OGRSpatialReference m_oSRS{};
    GDALGeoTransform m_gt{};
    bool bHaveGeoTransform;

    eProductType nProduct;

  public:
    TSXDataset();
    ~TSXDataset() override;

    int GetGCPCount() override;
    const OGRSpatialReference *GetGCPSpatialRef() const override;
    const GDAL_GCP *GetGCPs() override;

    CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;
    const OGRSpatialReference *GetSpatialRef() const override;

    static GDALDataset *Open(GDALOpenInfo *poOpenInfo);
    static int Identify(GDALOpenInfo *poOpenInfo);

  private:
    bool getGCPsFromGEOREF_XML(const char *pszGeorefFilename);
};

/************************************************************************/
/* ==================================================================== */
/*                                TSXRasterBand                           */
/* ==================================================================== */
/************************************************************************/

class TSXRasterBand final : public GDALPamRasterBand
{
    GDALDataset *poBand;
    ePolarization ePol;

  public:
    TSXRasterBand(TSXDataset *poDSIn, GDALDataType eDataType,
                  ePolarization ePol, GDALDataset *poBand);
    ~TSXRasterBand() override;

    CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage) override;

    static GDALDataset *Open(GDALOpenInfo *poOpenInfo);
};

/************************************************************************/
/*                            TSXRasterBand                             */
/************************************************************************/

TSXRasterBand::TSXRasterBand(TSXDataset *poDSIn, GDALDataType eDataTypeIn,
                             ePolarization ePolIn, GDALDataset *poBandIn)
    : poBand(poBandIn), ePol(ePolIn)
{
    poDS = poDSIn;
    eDataType = eDataTypeIn;

    switch (ePol)
    {
        case HH:
            SetMetadataItem("POLARIMETRIC_INTERP", "HH");
            break;
        case HV:
            SetMetadataItem("POLARIMETRIC_INTERP", "HV");
            break;
        case VH:
            SetMetadataItem("POLARIMETRIC_INTERP", "VH");
            break;
        case VV:
            SetMetadataItem("POLARIMETRIC_INTERP", "VV");
            break;
    }

    /* now setup the actual raster reader */
    GDALRasterBand *poSrcBand = poBandIn->GetRasterBand(1);
    poSrcBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
}

/************************************************************************/
/*                            TSXRasterBand()                           */
/************************************************************************/

TSXRasterBand::~TSXRasterBand()
{
    if (poBand != nullptr)
        GDALClose(reinterpret_cast<GDALRasterBandH>(poBand));
}

/************************************************************************/
/*                             IReadBlock()                             */
/************************************************************************/

CPLErr TSXRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
{
    int nRequestYSize;

    /* Check if the last strip is partial so we can avoid over-requesting */
    if ((nBlockYOff + 1) * nBlockYSize > nRasterYSize)
    {
        nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
        memset(pImage, 0,
               static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
                   nBlockXSize * nBlockYSize);
    }
    else
    {
        nRequestYSize = nBlockYSize;
    }

    /* Read Complex Data */
    if (eDataType == GDT_CInt16)
    {
        return poBand->RasterIO(
            GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
            nBlockXSize, nRequestYSize, pImage, nBlockXSize, nRequestYSize,
            GDT_CInt16, 1, nullptr, 4, nBlockXSize * 4, 0, nullptr);
    }

    // Detected Product
    return poBand->RasterIO(
        GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
        nBlockXSize, nRequestYSize, pImage, nBlockXSize, nRequestYSize,
        GDT_UInt16, 1, nullptr, 2, nBlockXSize * 2, 0, nullptr);
}

/************************************************************************/
/* ==================================================================== */
/*                                TSXDataset                                */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                             TSXDataset()                             */
/************************************************************************/

TSXDataset::TSXDataset()
    : nGCPCount(0), pasGCPList(nullptr), bHaveGeoTransform(false),
      nProduct(eUnknown)
{
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
    m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
}

/************************************************************************/
/*                            ~TSXDataset()                             */
/************************************************************************/

TSXDataset::~TSXDataset()
{
    FlushCache(true);

    if (nGCPCount > 0)
    {
        GDALDeinitGCPs(nGCPCount, pasGCPList);
        CPLFree(pasGCPList);
    }
}

/************************************************************************/
/*                              Identify()                              */
/************************************************************************/

int TSXDataset::Identify(GDALOpenInfo *poOpenInfo)
{
    if (poOpenInfo->fpL == nullptr || poOpenInfo->nHeaderBytes < 260)
    {
        if (poOpenInfo->bIsDirectory)
        {
            const CPLString osFilename = CPLFormCIFilenameSafe(
                poOpenInfo->pszFilename,
                CPLGetFilename(poOpenInfo->pszFilename), "xml");

            /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR
             * (TanDEM-X) or PAZ1_SAR (PAZ) */
            if (!(STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
                                 "TSX1_SAR") ||
                  STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
                                 "TDX1_SAR") ||
                  STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
                                 "PAZ1_SAR")))
                return 0;

            VSIStatBufL sStat;
            if (VSIStatL(osFilename, &sStat) == 0)
                return 1;
        }

        return 0;
    }

    /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR
     * (TanDEM-X) or PAZ1_SAR (PAZ) */
    if (!(STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
                         "TSX1_SAR") ||
          STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
                         "TDX1_SAR") ||
          STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
                         "PAZ1_SAR")))
        return 0;

    /* finally look for the <level1Product tag */
    if (!STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
                        "<level1Product"))
        return 0;

    return 1;
}

/************************************************************************/
/*                                getGCPsFromGEOREF_XML()               */
/*Reads georeferencing information from the TerraSAR-X GEOREF.XML file    */
/*and writes the information to the dataset's gcp list and projection     */
/*string.                                                                */
/*Returns true on success.                                                */
/************************************************************************/
bool TSXDataset::getGCPsFromGEOREF_XML(const char *pszGeorefFilename)
{
    // open GEOREF.xml
    CPLXMLNode *psGeorefData = CPLParseXMLFile(pszGeorefFilename);
    if (psGeorefData == nullptr)
        return false;

    // get the ellipsoid and semi-major, semi-minor axes
    OGRSpatialReference osr;
    CPLXMLNode *psSphere =
        CPLGetXMLNode(psGeorefData, "=geoReference.referenceFrames.sphere");
    if (psSphere != nullptr)
    {
        const char *pszEllipsoidName =
            CPLGetXMLValue(psSphere, "ellipsoidID", "");
        const double minor_axis =
            CPLAtof(CPLGetXMLValue(psSphere, "semiMinorAxis", "0.0"));
        const double major_axis =
            CPLAtof(CPLGetXMLValue(psSphere, "semiMajorAxis", "0.0"));
        // save datum parameters to the spatial reference
        if (EQUAL(pszEllipsoidName, "") || minor_axis == 0.0 ||
            major_axis == 0.0)
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "Warning- incomplete"
                     " ellipsoid information.  Using wgs-84 parameters.\n");
            osr.SetWellKnownGeogCS("WGS84");
        }
        else if (EQUAL(pszEllipsoidName, "WGS84"))
        {
            osr.SetWellKnownGeogCS("WGS84");
        }
        else
        {
            const double inv_flattening =
                major_axis / (major_axis - minor_axis);
            osr.SetGeogCS("", "", pszEllipsoidName, major_axis, inv_flattening);
        }
    }

    // get gcps
    CPLXMLNode *psGeolocationGrid =
        CPLGetXMLNode(psGeorefData, "=geoReference.geolocationGrid");
    if (psGeolocationGrid == nullptr)
    {
        CPLDestroyXMLNode(psGeorefData);
        return false;
    }
    nGCPCount = atoi(
        CPLGetXMLValue(psGeolocationGrid, "numberOfGridPoints.total", "0"));
    // count the gcps if the given count value is invalid
    CPLXMLNode *psNode = nullptr;
    if (nGCPCount <= 0)
    {
        for (psNode = psGeolocationGrid->psChild; psNode != nullptr;
             psNode = psNode->psNext)
            if (EQUAL(psNode->pszValue, "gridPoint"))
                nGCPCount++;
    }
    // if there are no gcps, fail
    if (nGCPCount <= 0)
    {
        CPLDestroyXMLNode(psGeorefData);
        return false;
    }

    // put some reasonable limits of the number of gcps
    if (nGCPCount > MAX_GCPS)
        nGCPCount = MAX_GCPS;
    // allocate memory for the gcps
    pasGCPList =
        static_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), nGCPCount));

    // loop through all gcps and set info

    // save the number allocated to ensure it does not run off the end of the
    // array
    const int gcps_allocated = nGCPCount;
    nGCPCount = 0;  // reset to zero and count
    // do a check on the grid point to make sure it has lat,long row, and column
    // it seems that only SSC products contain row, col - how to map lat long
    // otherwise?? for now fail if row and col are not present - just check the
    // first and assume the rest are the same
    for (psNode = psGeolocationGrid->psChild; psNode != nullptr;
         psNode = psNode->psNext)
    {
        if (!EQUAL(psNode->pszValue, "gridPoint"))
            continue;

        if (!strcmp(CPLGetXMLValue(psNode, "col", "error"), "error") ||
            !strcmp(CPLGetXMLValue(psNode, "row", "error"), "error") ||
            !strcmp(CPLGetXMLValue(psNode, "lon", "error"), "error") ||
            !strcmp(CPLGetXMLValue(psNode, "lat", "error"), "error"))
        {
            CPLDestroyXMLNode(psGeorefData);
            return false;
        }
    }
    for (psNode = psGeolocationGrid->psChild; psNode != nullptr;
         psNode = psNode->psNext)
    {
        // break out if the end of the array has been reached
        if (nGCPCount >= gcps_allocated)
        {
            CPLError(CE_Warning, CPLE_AppDefined,
                     "GDAL TSX driver: Truncating the number of GCPs.");
            break;
        }

        GDAL_GCP *psGCP = pasGCPList + nGCPCount;

        if (!EQUAL(psNode->pszValue, "gridPoint"))
            continue;

        nGCPCount++;

        char szID[32];
        snprintf(szID, sizeof(szID), "%d", nGCPCount);
        psGCP->pszId = CPLStrdup(szID);
        psGCP->pszInfo = CPLStrdup("");
        psGCP->dfGCPPixel = CPLAtof(CPLGetXMLValue(psNode, "col", "0"));
        psGCP->dfGCPLine = CPLAtof(CPLGetXMLValue(psNode, "row", "0"));
        psGCP->dfGCPX = CPLAtof(CPLGetXMLValue(psNode, "lon", ""));
        psGCP->dfGCPY = CPLAtof(CPLGetXMLValue(psNode, "lat", ""));
        // looks like height is in meters - should it be converted so xyz are
        // all on the same scale??
        psGCP->dfGCPZ = 0;
        // CPLAtof(CPLGetXMLValue(psNode,"height",""));
    }

    m_oGCPSRS = std::move(osr);

    CPLDestroyXMLNode(psGeorefData);

    return true;
}

/************************************************************************/
/*                                Open()                                */
/************************************************************************/

GDALDataset *TSXDataset::Open(GDALOpenInfo *poOpenInfo)
{
    /* -------------------------------------------------------------------- */
    /*      Is this a TerraSAR-X product file?                              */
    /* -------------------------------------------------------------------- */
    if (!TSXDataset::Identify(poOpenInfo))
    {
        return nullptr; /* nope */
    }

    /* -------------------------------------------------------------------- */
    /*      Confirm the requested access is supported.                      */
    /* -------------------------------------------------------------------- */
    if (poOpenInfo->eAccess == GA_Update)
    {
        ReportUpdateNotSupportedByDriver("TSX");
        return nullptr;
    }

    CPLString osFilename;

    if (poOpenInfo->bIsDirectory)
    {
        osFilename = CPLFormCIFilenameSafe(
            poOpenInfo->pszFilename, CPLGetFilename(poOpenInfo->pszFilename),
            "xml");
    }
    else
        osFilename = poOpenInfo->pszFilename;

    /* Ingest the XML */
    CPLXMLTreeCloser psData(CPLParseXMLFile(osFilename));
    if (psData == nullptr)
        return nullptr;

    /* find the product components */
    const CPLXMLNode *psComponents =
        CPLGetXMLNode(psData.get(), "=level1Product.productComponents");
    if (psComponents == nullptr)
    {
        CPLError(CE_Failure, CPLE_OpenFailed,
                 "Unable to find <productComponents> tag in file.\n");
        return nullptr;
    }

    /* find the product info tag */
    const CPLXMLNode *psProductInfo =
        CPLGetXMLNode(psData.get(), "=level1Product.productInfo");
    if (psProductInfo == nullptr)
    {
        CPLError(CE_Failure, CPLE_OpenFailed,
                 "Unable to find <productInfo> tag in file.\n");
        return nullptr;
    }

    /* -------------------------------------------------------------------- */
    /*      Create the dataset.                                             */
    /* -------------------------------------------------------------------- */

    auto poDS = std::make_unique<TSXDataset>();

    /* -------------------------------------------------------------------- */
    /*      Read in product info.                                           */
    /* -------------------------------------------------------------------- */

    poDS->SetMetadataItem(
        "SCENE_CENTRE_TIME",
        CPLGetXMLValue(psProductInfo,
                       "sceneInfo.sceneCenterCoord.azimuthTimeUTC", "unknown"));
    poDS->SetMetadataItem("OPERATIONAL_MODE",
                          CPLGetXMLValue(psProductInfo,
                                         "generationInfo.groundOperationsType",
                                         "unknown"));
    poDS->SetMetadataItem(
        "ORBIT_CYCLE",
        CPLGetXMLValue(psProductInfo, "missionInfo.orbitCycle", "unknown"));
    poDS->SetMetadataItem(
        "ABSOLUTE_ORBIT",
        CPLGetXMLValue(psProductInfo, "missionInfo.absOrbit", "unknown"));
    poDS->SetMetadataItem(
        "ORBIT_DIRECTION",
        CPLGetXMLValue(psProductInfo, "missionInfo.orbitDirection", "unknown"));
    poDS->SetMetadataItem("IMAGING_MODE",
                          CPLGetXMLValue(psProductInfo,
                                         "acquisitionInfo.imagingMode",
                                         "unknown"));
    poDS->SetMetadataItem("PRODUCT_VARIANT",
                          CPLGetXMLValue(psProductInfo,
                                         "productVariantInfo.productVariant",
                                         "unknown"));
    std::string osDataType =
        CPLGetXMLValue(psProductInfo, "imageDataInfo.imageDataType", "unknown");
    poDS->SetMetadataItem("IMAGE_TYPE", osDataType.c_str());

    /* Get raster information */
    int nRows = atoi(CPLGetXMLValue(
        psProductInfo, "imageDataInfo.imageRaster.numberOfRows", ""));
    int nCols = atoi(CPLGetXMLValue(
        psProductInfo, "imageDataInfo.imageRaster.numberOfColumns", ""));

    poDS->nRasterXSize = nCols;
    poDS->nRasterYSize = nRows;

    poDS->SetMetadataItem("ROW_SPACING",
                          CPLGetXMLValue(psProductInfo,
                                         "imageDataInfo.imageRaster.rowSpacing",
                                         "unknown"));
    poDS->SetMetadataItem(
        "COL_SPACING",
        CPLGetXMLValue(psProductInfo, "imageDataInfo.imageRaster.columnSpacing",
                       "unknown"));
    poDS->SetMetadataItem(
        "COL_SPACING_UNITS",
        CPLGetXMLValue(psProductInfo,
                       "imageDataInfo.imageRaster.columnSpacing.units",
                       "unknown"));

    /* Get equivalent number of looks */
    poDS->SetMetadataItem(
        "AZIMUTH_LOOKS",
        CPLGetXMLValue(psProductInfo, "imageDataInfo.imageRaster.azimuthLooks",
                       "unknown"));
    poDS->SetMetadataItem("RANGE_LOOKS",
                          CPLGetXMLValue(psProductInfo,
                                         "imageDataInfo.imageRaster.rangeLooks",
                                         "unknown"));

    const char *pszProductVariant = CPLGetXMLValue(
        psProductInfo, "productVariantInfo.productVariant", "unknown");

    poDS->SetMetadataItem("PRODUCT_VARIANT", pszProductVariant);

    /* Determine what product variant this is */
    if (STARTS_WITH_CI(pszProductVariant, "SSC"))
        poDS->nProduct = eSSC;
    else if (STARTS_WITH_CI(pszProductVariant, "MGD"))
        poDS->nProduct = eMGD;
    else if (STARTS_WITH_CI(pszProductVariant, "EEC"))
        poDS->nProduct = eEEC;
    else if (STARTS_WITH_CI(pszProductVariant, "GEC"))
        poDS->nProduct = eGEC;
    else
        poDS->nProduct = eUnknown;

    /* Start reading in the product components */
    std::string osGeorefFile;
    CPLErr geoTransformErr = CE_Failure;
    for (CPLXMLNode *psComponent = psComponents->psChild;
         psComponent != nullptr; psComponent = psComponent->psNext)
    {
        const char *pszType = nullptr;
        const std::string osFilePath = GetFilePath(psComponent, &pszType);
        if (CPLHasPathTraversal(osFilePath.c_str()))
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Path traversal detected in %s", osFilePath.c_str());
            return nullptr;
        }
        std::string osPath = CPLFormFilenameSafe(
            CPLGetDirnameSafe(osFilename).c_str(), osFilePath.c_str(), "");
        const char *pszPolLayer = CPLGetXMLValue(psComponent, "polLayer", " ");

        if (!STARTS_WITH_CI(pszType, " "))
        {
            if (STARTS_WITH_CI(pszType, "MAPPING_GRID"))
            {
                /* the mapping grid... save as a metadata item this path */
                poDS->SetMetadataItem("MAPPING_GRID", osPath.c_str());
            }
            else if (STARTS_WITH_CI(pszType, "GEOREF"))
            {
                /* save the path to the georef data for later use */
                osGeorefFile = std::move(osPath);
            }
        }
        else if (!STARTS_WITH_CI(pszPolLayer, " ") &&
                 STARTS_WITH_CI(psComponent->pszValue, "imageData"))
        {
            /* determine the polarization of this band */
            ePolarization ePol;
            if (STARTS_WITH_CI(pszPolLayer, "HH"))
            {
                ePol = HH;
            }
            else if (STARTS_WITH_CI(pszPolLayer, "HV"))
            {
                ePol = HV;
            }
            else if (STARTS_WITH_CI(pszPolLayer, "VH"))
            {
                ePol = VH;
            }
            else
            {
                ePol = VV;
            }

            GDALDataType eDataType =
                STARTS_WITH_CI(osDataType.c_str(), "COMPLEX") ? GDT_CInt16
                                                              : GDT_UInt16;

            /* try opening the file that represents that band */
            GDALDataset *poBandData =
                GDALDataset::FromHandle(GDALOpen(osPath.c_str(), GA_ReadOnly));
            if (poBandData != nullptr)
            {
                TSXRasterBand *poBand =
                    new TSXRasterBand(poDS.get(), eDataType, ePol, poBandData);
                poDS->SetBand(poDS->GetRasterCount() + 1, poBand);

                // copy georeferencing info from the band
                // need error checking??
                // it will just save the info from the last band
                const auto poSrcSRS = poBandData->GetSpatialRef();
                if (poSrcSRS)
                    poDS->m_oSRS = *poSrcSRS;

                geoTransformErr = poBandData->GetGeoTransform(poDS->m_gt);
            }
        }
    }

    // now check if there is a geotransform
    if (!poDS->m_oSRS.IsEmpty() && geoTransformErr == CE_None)
    {
        poDS->bHaveGeoTransform = TRUE;
    }
    else
    {
        poDS->bHaveGeoTransform = FALSE;
        poDS->m_oSRS.Clear();
        poDS->m_gt = GDALGeoTransform();
    }

    /* -------------------------------------------------------------------- */
    /*      Check and set matrix representation.                            */
    /* -------------------------------------------------------------------- */

    if (poDS->GetRasterCount() == 4)
    {
        poDS->SetMetadataItem("MATRIX_REPRESENTATION", "SCATTERING");
    }

    /* -------------------------------------------------------------------- */
    /*      Read the four corners and centre GCPs in                        */
    /* -------------------------------------------------------------------- */

    const CPLXMLNode *psSceneInfo =
        CPLGetXMLNode(psData.get(), "=level1Product.productInfo.sceneInfo");
    if (psSceneInfo != nullptr)
    {
        /* extract the GCPs from the provided file */
        bool success = false;
        if (!osGeorefFile.empty())
            success = poDS->getGCPsFromGEOREF_XML(osGeorefFile.c_str());

        // if the gcp's cannot be extracted from the georef file, try to get the
        // corner coordinates for now just SSC because the others don't have
        // refColumn and refRow
        if (!success && poDS->nProduct == eSSC)
        {
            int nGCP = 0;
            double dfAvgHeight = CPLAtof(
                CPLGetXMLValue(psSceneInfo, "sceneAverageHeight", "0.0"));

            // count and allocate gcps - there should be five - 4 corners and a
            // centre
            poDS->nGCPCount = 0;
            const CPLXMLNode *psNode = psSceneInfo->psChild;
            for (; psNode != nullptr; psNode = psNode->psNext)
            {
                if (!EQUAL(psNode->pszValue, "sceneCenterCoord") &&
                    !EQUAL(psNode->pszValue, "sceneCornerCoord"))
                    continue;

                poDS->nGCPCount++;
            }
            if (poDS->nGCPCount > 0)
            {
                poDS->pasGCPList =
                    (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount);

                /* iterate over GCPs */
                for (psNode = psSceneInfo->psChild; psNode != nullptr;
                     psNode = psNode->psNext)
                {
                    GDAL_GCP *psGCP = poDS->pasGCPList + nGCP;

                    if (!EQUAL(psNode->pszValue, "sceneCenterCoord") &&
                        !EQUAL(psNode->pszValue, "sceneCornerCoord"))
                        continue;

                    psGCP->dfGCPPixel =
                        CPLAtof(CPLGetXMLValue(psNode, "refColumn", "0.0"));
                    psGCP->dfGCPLine =
                        CPLAtof(CPLGetXMLValue(psNode, "refRow", "0.0"));
                    psGCP->dfGCPX =
                        CPLAtof(CPLGetXMLValue(psNode, "lon", "0.0"));
                    psGCP->dfGCPY =
                        CPLAtof(CPLGetXMLValue(psNode, "lat", "0.0"));
                    psGCP->dfGCPZ = dfAvgHeight;
                    psGCP->pszId = CPLStrdup(CPLSPrintf("%d", nGCP));
                    psGCP->pszInfo = CPLStrdup("");

                    nGCP++;
                }

                // set the projection string - the fields are lat/long - seems
                // to be WGS84 datum
                poDS->m_oGCPSRS.SetWellKnownGeogCS("WGS84");
            }
        }

        // gcps override geotransform - does it make sense to have both??
        if (poDS->nGCPCount > 0)
        {
            poDS->bHaveGeoTransform = FALSE;
            poDS->m_oSRS.Clear();
            poDS->m_gt = GDALGeoTransform();
        }
    }
    else
    {
        CPLError(CE_Warning, CPLE_AppDefined,
                 "Unable to find sceneInfo tag in XML document. "
                 "Proceeding with caution.");
    }

    /* -------------------------------------------------------------------- */
    /*      Initialize any PAM information.                                 */
    /* -------------------------------------------------------------------- */
    poDS->SetDescription(poOpenInfo->pszFilename);
    poDS->TryLoadXML();

    /* -------------------------------------------------------------------- */
    /*      Check for overviews.                                            */
    /* -------------------------------------------------------------------- */
    poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);

    return poDS.release();
}

/************************************************************************/
/*                            GetGCPCount()                             */
/************************************************************************/

int TSXDataset::GetGCPCount()
{
    return nGCPCount;
}

/************************************************************************/
/*                          GetGCPSpatialRef()                          */
/************************************************************************/

const OGRSpatialReference *TSXDataset::GetGCPSpatialRef() const
{
    return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
}

/************************************************************************/
/*                               GetGCPs()                              */
/************************************************************************/

const GDAL_GCP *TSXDataset::GetGCPs()
{
    return pasGCPList;
}

/************************************************************************/
/*                          GetSpatialRef()                             */
/************************************************************************/

const OGRSpatialReference *TSXDataset::GetSpatialRef() const

{
    return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
}

/************************************************************************/
/*                               GetGeotransform()                      */
/************************************************************************/
CPLErr TSXDataset::GetGeoTransform(GDALGeoTransform &gt) const
{
    gt = m_gt;

    if (bHaveGeoTransform)
        return CE_None;

    return CE_Failure;
}

/************************************************************************/
/*                         GDALRegister_TSX()                           */
/************************************************************************/

void GDALRegister_TSX()
{
    if (GDALGetDriverByName("TSX") != nullptr)
        return;

    GDALDriver *poDriver = new GDALDriver();

    poDriver->SetDescription("TSX");
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "TerraSAR-X Product");
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/tsx.html");
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");

    poDriver->pfnOpen = TSXDataset::Open;
    poDriver->pfnIdentify = TSXDataset::Identify;

    GetGDALDriverManager()->RegisterDriver(poDriver);
}
