/******************************************************************************

 *

 * 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);

}