Ticket #4160: csk_hdf5.patch

File csk_hdf5.patch, 26.5 KB (added by alexmantaut, 5 years ago)

Patch to get georreference and geotransform into (1st version)

  • frmts/hdf5/hdf5dataset.cpp

    diff -rupN /c/Documents and Settings/alexmantaut/Mis documentos/codigo suelto/CSKHDF5/hdf5_clean/hdf5/hdf5dataset.cpp /c/Documents and Settings/alexmantaut/Mis documentos/codigo suelto/CSKHDF5/hdf5/hdf5dataset.cpp
     
    3838#include "gdal_priv.h"
    3939#include "cpl_string.h"
    4040#include "hdf5dataset.h"
     41#include <ogr_spatialref.h>
    4142
    4243CPL_CVSID("$Id: hdf5dataset.cpp 20437 2010-08-25 17:20:02Z chaitanya $");
    4344
    herr_t HDF5AttrIterate( hid_t hH5ObjID,  
    611612
    612613    poDS = (HDF5Dataset *) pDS;
    613614
    614     // Convert "/" into "_" fro the path component
    615     const char* pszPath = poDS->poH5CurrentObject->pszUnderscorePath;
    616     if(pszPath != NULL && strlen(pszPath) > 0)
    617     {
    618         papszTokens = CSLTokenizeString2( pszPath, "/", CSLT_HONOURSTRINGS );
     615    osKey = HDF5ReplaceGroupSeparators(poDS->poH5CurrentObject->pszUnderscorePath);
    619616
    620         for( int i = 0; papszTokens != NULL && papszTokens[i] != NULL; ++i )
    621         {
    622             if( i != 0)
    623                 osKey += '_';
    624             osKey += papszTokens[i];
    625         }
    626         CSLDestroy( papszTokens );
    627     }
    628    
    629617    // Convert whitespaces into "_" for the attribute name component
    630618    papszTokens = CSLTokenizeString2( pszAttrName, " ",
    631619                            CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES );
    CPLErr HDF5Dataset::ReadGlobalAttributes 
    10181006    HDF5ListGroupObjects( poRootGroup, bSUBDATASET );
    10191007    return CE_None;
    10201008}
     1009
     1010char *HDF5ReplaceGroupSeparators(const char * pszPath)
     1011{
     1012   char          **papszTokens;
     1013   CPLString       osKey;
     1014
     1015   // Convert "/" into "_" fro the path component
     1016   //const char* pszPath = poDS->poH5CurrentObject->pszUnderscorePath;
     1017
     1018        if(pszPath != NULL && strlen(pszPath) > 0)
     1019   {
     1020        papszTokens = CSLTokenizeString2( pszPath, "/", CSLT_HONOURSTRINGS );
     1021
     1022        for( int i = 0; papszTokens != NULL && papszTokens[i] != NULL; ++i )
     1023        {
     1024            if( i != 0)
     1025                osKey += HDF5_INTERAL_GROUP_SEPARATOR;
     1026            osKey += papszTokens[i];
     1027        }
     1028        CSLDestroy( papszTokens );
     1029   }
     1030
     1031   return CPLStrdup(osKey);
     1032}
     1033
  • frmts/hdf5/hdf5dataset.h

    diff -rupN /c/Documents and Settings/alexmantaut/Mis documentos/codigo suelto/CSKHDF5/hdf5_clean/hdf5/hdf5dataset.h /c/Documents and Settings/alexmantaut/Mis documentos/codigo suelto/CSKHDF5/hdf5/hdf5dataset.h
     
    3232
    3333#include "gdal_pam.h"
    3434#include "cpl_list.h"
     35#define HDF5_INTERAL_GROUP_SEPARATOR '_'
    3536
    3637typedef struct HDF5GroupObjects {
    3738  char         *pszName;
    class HDF5Dataset : public GDALPamDatase 
    8081
    8182  HDF5GroupObjects* HDF5FindDatasetObjects( HDF5GroupObjects *, const char * );
    8283  HDF5GroupObjects* HDF5FindDatasetObjectsbyPath( HDF5GroupObjects *, const char * );
    83   char* CreatePath(HDF5GroupObjects *);
     84
    8485  void DestroyH5Objects(HDF5GroupObjects *);
    8586 
    8687  GDALDataType GetDataType(hid_t);
    8788  const char * GetDataTypeName(hid_t);
     89//  char *FormatAttributePath(const char *);
    8890
    8991  public:
    9092
    class HDF5Dataset : public GDALPamDatase 
    9698   
    9799  static GDALDataset *Open(GDALOpenInfo *);
    98100  static int Identify(GDALOpenInfo *);
    99 };
     101 
     102//  const char *GetProjectionRef();
    100103
     104};
    101105
     106char *HDF5ReplaceGroupSeparators(const char *);
    102107
    103108#endif /* _HDF5DATASET_H_INCLUDED_ */
    104109
  • frmts/hdf5/hdf5imagedataset.cpp

    diff -rupN /c/Documents and Settings/alexmantaut/Mis documentos/codigo suelto/CSKHDF5/hdf5_clean/hdf5/hdf5imagedataset.cpp /c/Documents and Settings/alexmantaut/Mis documentos/codigo suelto/CSKHDF5/hdf5/hdf5imagedataset.cpp
     
    2727 * DEALINGS IN THE SOFTWARE.
    2828 ****************************************************************************/
    2929
    30 #define H5_USE_16_API
    31 
    32 #include "hdf5.h"
    33 
    34 #include "gdal_pam.h"
    35 #include "gdal_priv.h"
    36 #include "cpl_string.h"
    37 #include "hdf5dataset.h"
    38 #include "ogr_spatialref.h"
    39 
    40 CPL_CVSID("$Id: hdf5imagedataset.cpp 21169 2010-11-24 15:57:11Z warmerdam $");
    41 
    42 CPL_C_START
    43 void    GDALRegister_HDF5Image(void);
    44 CPL_C_END
    45 
    46 /* release 1.6.3 or 1.6.4 changed the type of count in some api functions */
    47 
    48 #if H5_VERS_MAJOR == 1 && H5_VERS_MINOR <= 6 \
    49        && (H5_VERS_MINOR < 6 || H5_VERS_RELEASE < 3)
    50 #  define H5OFFSET_TYPE hssize_t
    51 #else
    52 #  define H5OFFSET_TYPE  hsize_t
    53 #endif
    54 
    55 class HDF5ImageDataset : public HDF5Dataset
    56 {
    57 
    58     friend class HDF5ImageRasterBand;
    59 
    60     char        *pszProjection;
    61     char        *pszGCPProjection;
    62     GDAL_GCP    *pasGCPList;
    63     int         nGCPCount;
    64     OGRSpatialReference oSRS;
    65 
    66     hsize_t      *dims,*maxdims;
    67     HDF5GroupObjects *poH5Objects;
    68     int          ndims,dimensions;
    69     hid_t        dataset_id;
    70     hid_t        dataspace_id;
    71     hsize_t      size;
    72     int          address;
    73     hid_t        datatype;
    74     hid_t        native;
    75     H5T_class_t  clas;
    76 
    77 public:
    78     HDF5ImageDataset();
    79     ~HDF5ImageDataset();
    80    
    81     CPLErr CreateProjections( );
    82     static GDALDataset  *Open( GDALOpenInfo * );
    83     static int           Identify( GDALOpenInfo * );
    84 
    85     const char          *GetProjectionRef();
    86     virtual int         GetGCPCount( );
    87     virtual const char  *GetGCPProjection();
    88     virtual const GDAL_GCP *GetGCPs( );
    89    
    90 };
     30#include "hdf5imagedataset.h"
    9131
    9232/************************************************************************/
    9333/* ==================================================================== */
    HDF5ImageDataset::HDF5ImageDataset() 
    10949    dims            = NULL;
    11050    maxdims         = NULL;
    11151    papszMetadata   = NULL;
     52    adfGeoTransform[0] = 0.0;
     53    adfGeoTransform[1] = 1.0;
     54    adfGeoTransform[2] = 0.0;
     55    adfGeoTransform[3] = 0.0;
     56    adfGeoTransform[4] = 0.0;
     57    adfGeoTransform[5] = 1.0;
     58        iSubdatasetType     = UNKNOWN_PRODUCT;
     59    bHasGeoTransform =  false;
    11260}
    11361
    11462/************************************************************************/
    GDALDataset *HDF5ImageDataset::Open( GDA 
    440388    poDS->bIsHDFEOS=TRUE;
    441389    poDS->ReadGlobalAttributes( FALSE );
    442390 
     391    //Fecth the open subdataset
     392    HDF5GroupObjects *poH5Aux = poDS->HDF5FindDatasetObjectsbyPath( poDS->poH5RootGroup,
     393                                                    (char *)poDS->GetSubdatasetName() );
     394
     395    //Get the metadata asociated to the subdataset
     396    if(poH5Aux!=NULL)
     397         poDS->CreateMetadata( poH5Aux, H5G_DATASET );
     398
    443399/* -------------------------------------------------------------------- */
    444400/*      Create HDF5 Data Hierarchy in a link list                       */
    445401/* -------------------------------------------------------------------- */
    GDALDataset *HDF5ImageDataset::Open( GDA 
    487443                poBand->SetNoDataValue( 255 );
    488444    }
    489445
    490     poDS->CreateProjections( );
    491 
     446        //Get all the metadata
    492447    poDS->SetMetadata( poDS->papszMetadata );
    493448
     449    //Check if the hdf5 is a well known product type
     450    poDS->IdentifyProductType();
     451
     452    poDS->CreateProjections( );
     453
    494454/* -------------------------------------------------------------------- */
    495455/*      Setup/check for pam .aux.xml.                                   */
    496456/* -------------------------------------------------------------------- */
    void GDALRegister_HDF5Image( ) 
    537497/************************************************************************/
    538498CPLErr HDF5ImageDataset::CreateProjections()
    539499{
    540 #define NBGCPLAT 100
    541 #define NBGCPLON 30
    542 
    543     hid_t LatitudeDatasetID  = -1;
    544     hid_t LongitudeDatasetID = -1;
    545     hid_t LatitudeDataspaceID;
    546     hid_t LongitudeDataspaceID;
    547     float* Latitude;
    548     float* Longitude;
    549     int    i,j;
    550     int   nDeltaLat;
    551     int   nDeltaLon;
    552 
    553     nDeltaLat = nRasterYSize / NBGCPLAT;
    554     nDeltaLon = nRasterXSize / NBGCPLON;
    555 
    556     if( nDeltaLat == 0 || nDeltaLon == 0 )
    557         return CE_None;
    558 
    559 /* -------------------------------------------------------------------- */
    560 /*      Create HDF5 Data Hierarchy in a link list                       */
    561 /* -------------------------------------------------------------------- */
    562     poH5Objects=HDF5FindDatasetObjects( poH5RootGroup,  "Latitude" );
    563     if( !poH5Objects ) {
    564         return CE_None;
    565     }
    566 /* -------------------------------------------------------------------- */
    567 /*      The Lattitude and Longitude arrays must have a rank of 2 to     */
    568 /*      retrieve GCPs.                                                  */
    569 /* -------------------------------------------------------------------- */
    570     if( poH5Objects->nRank != 2 ) {
    571         return CE_None;
    572     }
    573    
    574 /* -------------------------------------------------------------------- */
    575 /*      Retrieve HDF5 data information                                  */
    576 /* -------------------------------------------------------------------- */
    577     LatitudeDatasetID   = H5Dopen( hHDF5,poH5Objects->pszPath );
    578     LatitudeDataspaceID = H5Dget_space( dataset_id );                       
    579 
    580     poH5Objects=HDF5FindDatasetObjects( poH5RootGroup, "Longitude" );
    581     LongitudeDatasetID   = H5Dopen( hHDF5,poH5Objects->pszPath );
    582     LongitudeDataspaceID = H5Dget_space( dataset_id );                       
    583 
    584     if( ( LatitudeDatasetID > 0 ) && ( LongitudeDatasetID > 0) ) {
    585        
    586         Latitude         = ( float * ) CPLCalloc(  nRasterYSize*nRasterXSize,
    587                                                 sizeof( float ) );
    588         Longitude         = ( float * ) CPLCalloc( nRasterYSize*nRasterXSize,
    589                                                  sizeof( float ) );
    590         memset( Latitude, 0, nRasterXSize*nRasterYSize*sizeof(  float ) );
    591         memset( Longitude, 0, nRasterXSize*nRasterYSize*sizeof( float ) );
     500        switch(iSubdatasetType)
     501        {
     502        case CSK_PRODUCT:
     503        {
     504                const char *osMissionLevel;
     505                int productType;
     506
     507                if(GetMetadataItem("Product_Type")!=NULL)
     508                {
     509                        //Get the format's level
     510                        osMissionLevel = HDF5Dataset::GetMetadataItem("Product_Type");
     511
     512                        if(EQUALN(osMissionLevel,"RAW",3))
     513                                productType  = PROD_CSK_L0;
     514
     515                        if(EQUALN(osMissionLevel,"SSC",3))
     516                                productType  = PROD_CSK_L1A;
     517
     518                        if(EQUALN(osMissionLevel,"DGM",3))
     519                                productType  = PROD_CSK_L1B;
     520
     521                        if(EQUALN(osMissionLevel,"GEC",3))
     522                                productType  = PROD_CSK_L1C;
     523
     524                        if(EQUALN(osMissionLevel,"GTC",3))
     525                                productType  = PROD_CSK_L1D;
     526                }
     527                else
     528                        productType =PROD_UNKNOWN;
     529                       
     530                CaptureCSKGeoTransform(productType);
     531                CaptureCSKGeolocation(productType);
     532               
     533                break;
     534        }
     535        case UNKNOWN_PRODUCT:
     536                {
     537                #define NBGCPLAT 100
     538                #define NBGCPLON 30
     539
     540                 hid_t LatitudeDatasetID  = -1;
     541                 hid_t LongitudeDatasetID = -1;
     542                 hid_t LatitudeDataspaceID;
     543                 hid_t LongitudeDataspaceID;
     544                 float* Latitude;
     545                 float* Longitude;
     546                 int    i,j;
     547                 int   nDeltaLat;
     548                 int   nDeltaLon;
    592549       
    593         H5Dread ( LatitudeDatasetID,
    594                   H5T_NATIVE_FLOAT,
    595                   H5S_ALL,
    596                   H5S_ALL,
    597                   H5P_DEFAULT,
    598                   Latitude );
     550                 nDeltaLat = nRasterYSize / NBGCPLAT;
     551                 nDeltaLon = nRasterXSize / NBGCPLON;
    599552       
    600         H5Dread ( LongitudeDatasetID,
    601                   H5T_NATIVE_FLOAT,
    602                   H5S_ALL,
    603                   H5S_ALL,
    604                   H5P_DEFAULT,
    605                   Longitude );
     553                 if( nDeltaLat == 0 || nDeltaLon == 0 )
     554                          return CE_None;
    606555       
    607         oSRS.SetWellKnownGeogCS( "WGS84" );
    608         CPLFree(pszProjection);
    609         CPLFree(pszGCPProjection);
    610         oSRS.exportToWkt( &pszProjection );
    611         oSRS.exportToWkt( &pszGCPProjection );
     556        /* -------------------------------------------------------------------- */
     557        /*      Create HDF5 Data Hierarchy in a link list                       */
     558        /* -------------------------------------------------------------------- */
     559                 poH5Objects=HDF5FindDatasetObjects( poH5RootGroup,  "Latitude" );
     560                 if( !poH5Objects ) {
     561                return CE_None;
     562                 }
     563        /* -------------------------------------------------------------------- */
     564        /*      The Lattitude and Longitude arrays must have a rank of 2 to     */
     565        /*      retrieve GCPs.                                                  */
     566        /* -------------------------------------------------------------------- */
     567                 if( poH5Objects->nRank != 2 ) {
     568                return CE_None;
     569                 }
     570
     571        /* -------------------------------------------------------------------- */
     572        /*      Retrieve HDF5 data information                                  */
     573        /* -------------------------------------------------------------------- */
     574                 LatitudeDatasetID   = H5Dopen( hHDF5,poH5Objects->pszPath );
     575                 LatitudeDataspaceID = H5Dget_space( dataset_id );
    612576       
    613 /* -------------------------------------------------------------------- */
    614 /*  Fill the GCPs list.                                                 */
    615 /* -------------------------------------------------------------------- */
    616         nGCPCount = nRasterYSize/nDeltaLat * nRasterXSize/nDeltaLon;
    617 
    618         pasGCPList = ( GDAL_GCP * )
    619             CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
     577                 poH5Objects=HDF5FindDatasetObjects( poH5RootGroup, "Longitude" );
     578                 LongitudeDatasetID   = H5Dopen( hHDF5,poH5Objects->pszPath );
     579                 LongitudeDataspaceID = H5Dget_space( dataset_id );
    620580       
    621         GDALInitGCPs( nGCPCount, pasGCPList );
    622         int k=0;
     581                 if( ( LatitudeDatasetID > 0 ) && ( LongitudeDatasetID > 0) ) {
    623582
    624         int nYLimit = ((int)nRasterYSize/nDeltaLat) * nDeltaLat;
    625         int nXLimit = ((int)nRasterXSize/nDeltaLon) * nDeltaLon;
    626         for( j = 0; j < nYLimit; j+=nDeltaLat ) {
    627             for( i = 0; i < nXLimit; i+=nDeltaLon ) {
    628                 int iGCP =  j * nRasterXSize + i;
    629                 pasGCPList[k].dfGCPX = ( double ) Longitude[iGCP]+180.0;
    630                 pasGCPList[k].dfGCPY = ( double ) Latitude[iGCP];
     583                Latitude         = ( float * ) CPLCalloc(  nRasterYSize*nRasterXSize,
     584                                                        sizeof( float ) );
     585                Longitude         = ( float * ) CPLCalloc( nRasterYSize*nRasterXSize,
     586                                                         sizeof( float ) );
     587                memset( Latitude, 0, nRasterXSize*nRasterYSize*sizeof(  float ) );
     588                memset( Longitude, 0, nRasterXSize*nRasterYSize*sizeof( float ) );
     589
     590                H5Dread ( LatitudeDatasetID,
     591                          H5T_NATIVE_FLOAT,
     592                          H5S_ALL,
     593                          H5S_ALL,
     594                          H5P_DEFAULT,
     595                          Latitude );
     596
     597                H5Dread ( LongitudeDatasetID,
     598                          H5T_NATIVE_FLOAT,
     599                          H5S_ALL,
     600                          H5S_ALL,
     601                          H5P_DEFAULT,
     602                          Longitude );
     603
     604                oSRS.SetWellKnownGeogCS( "WGS84" );
     605                          CPLFree(pszProjection);
     606                          CPLFree(pszGCPProjection);
     607                oSRS.exportToWkt( &pszProjection );
     608                oSRS.exportToWkt( &pszGCPProjection );
    631609               
    632                 pasGCPList[k].dfGCPPixel = i + 0.5;
    633                 pasGCPList[k++].dfGCPLine =  j + 0.5;
     610        /* -------------------------------------------------------------------- */
     611        /*  Fill the GCPs list.                                                 */
     612        /* -------------------------------------------------------------------- */
     613                nGCPCount = nRasterYSize/nDeltaLat * nRasterXSize/nDeltaLon;
     614
     615                pasGCPList = ( GDAL_GCP * )
     616                         CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
    634617               
    635             }
    636         }
     618                GDALInitGCPs( nGCPCount, pasGCPList );
     619                int k=0;
    637620       
    638         CPLFree( Latitude );
    639         CPLFree( Longitude );
    640     }
    641     return CE_None;
    642 
     621                int nYLimit = ((int)nRasterYSize/nDeltaLat) * nDeltaLat;
     622                int nXLimit = ((int)nRasterXSize/nDeltaLon) * nDeltaLon;
     623                for( j = 0; j < nYLimit; j+=nDeltaLat ) {
     624                         for( i = 0; i < nXLimit; i+=nDeltaLon ) {
     625                        int iGCP =  j * nRasterXSize + i;
     626                        pasGCPList[k].dfGCPX = ( double ) Longitude[iGCP]+180.0;
     627                        pasGCPList[k].dfGCPY = ( double ) Latitude[iGCP];
     628
     629                        pasGCPList[k].dfGCPPixel = i + 0.5;
     630                        pasGCPList[k++].dfGCPLine =  j + 0.5;
     631
     632                         }
     633                }
     634
     635                CPLFree( Latitude );
     636                CPLFree( Longitude );
     637                 }
     638                 
     639                break;
     640                }
     641        }
     642        return CE_None;
    643643}
    644644/************************************************************************/
    645645/*                          GetProjectionRef()                          */
    const GDAL_GCP *HDF5ImageDataset::GetGCP 
    694694
    695695
    696696
     697CPLErr HDF5ImageDataset::GetGeoTransform( double * padfTransform )
     698{
     699        CPLErr retVal;
     700       
     701        if ( bHasGeoTransform )
     702        {
     703                memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
     704                retVal = CE_None;
     705        }
     706        else
     707                retVal = CE_Failure;
     708
     709        return retVal;
     710}
     711
     712/**
     713 * Identify if the subdataset has a known product format
     714 * It stores a product identifier in iSubdatasetType,
     715 * UNKNOWN_PRODUCT, if it isn't a recognizable format.
     716 */
     717void HDF5ImageDataset::IdentifyProductType()
     718{
     719        iSubdatasetType = UNKNOWN_PRODUCT;
     720
     721/************************************************************************/
     722/*                               COSMO-SKYMED                           */
     723/************************************************************************/
     724   const char *pszMissionId;
     725
     726        //Get the Mission Id as a char *, because the
     727        //field may not exist
     728        pszMissionId = HDF5Dataset::GetMetadataItem("Mission_ID");
     729
     730        //If there is a Mission_ID field
     731         if(pszMissionId != NULL)
     732                 //Check if the mission type is CSK
     733                 if(EQUAL(pszMissionId,"CSK"))
     734                         iSubdatasetType = CSK_PRODUCT;
     735}
     736
     737/**
     738 * Captures Geolocation information from a COSMO-SKYMED
     739 * file.
     740 * The geoid will allways be WGS84
     741 * The projection type may be UTM or UPS, depending on the
     742 * latitude from the center of the image.
     743 * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
     744 */
     745void HDF5ImageDataset::CaptureCSKGeolocation(int iProductType)
     746{
     747        //Check if all the metadata attributes are present
     748        if(GetMetadataItem("Projection_ID")== NULL ||
     749                GetMetadataItem("Map_Projection_False_East-North")== NULL ||
     750                GetMetadataItem("Map_Projection_Centre") ==NULL ||
     751                GetMetadataItem("Map_Projection_Scale_Factor") == NULL)
     752        {
     753
     754                pszProjection = CPLStrdup("");
     755                pszGCPProjection = CPLStrdup("");
     756                CPLError( CE_Failure, CPLE_OpenFailed,
     757                                         "The CSK hdf5 file geolocation information is malformed\n" );
     758        }
     759        else
     760        {
     761
     762                //Fetch projection Type
     763                CPLString osProjectionID = GetMetadataItem("Projection_ID");
     764
     765                //Fetch projection properties
     766                CPLString osProjectionFalseEastNorth = GetMetadataItem("Map_Projection_False_East-North");
     767                CPLString osProjectionScale = GetMetadataItem("Map_Projection_Scale_Factor");
     768                double dfProjFalseEast, dfProjFalseNorth, dfProjScaleFactor;
     769                CPLString osProjCenter = GetMetadataItem("Map_Projection_Centre");
     770                char ** papszCenterCoord = CSLTokenizeString2(osProjCenter," ",
     771                                                                                                CSLT_STRIPLEADSPACES);
     772
     773                double dfCenterLong, dfCenterLat;
     774
     775                //Convert projection parameters to double
     776                char ** papszProjFalseEastNorth = CSLTokenizeString2(osProjectionFalseEastNorth," ",
     777                                                                                                                CSLT_STRIPLEADSPACES);
     778
     779                dfProjFalseEast = CPLAtof(papszProjFalseEastNorth[0]);
     780                dfProjFalseNorth = CPLAtof(papszProjFalseEastNorth[1]);
     781                dfProjScaleFactor =  CPLAtof(osProjectionScale);
     782                dfCenterLong = CPLAtof(papszCenterCoord[1]);
     783                dfCenterLat = CPLAtof(papszCenterCoord[0]);
     784
     785                CPLFree(papszCenterCoord);
     786                CPLFree(papszProjFalseEastNorth);
     787
     788                //Set the ellipsoid to WGS84
     789                oSRS.SetWellKnownGeogCS( "WGS84" );
     790
     791                //If the projection is UTM
     792                if(EQUAL(osProjectionID,"UTM"))
     793                {
     794                        oSRS.SetTM(dfCenterLat,dfCenterLong,dfProjScaleFactor,
     795                                                  dfProjFalseEast,dfProjFalseNorth);
     796                }
     797                else
     798                {
     799                        //TODO Test! I didn't had any UPS projected files to test!
     800                        //If the projection is UPS
     801                        if(EQUAL(osProjectionID,"UPS"))
     802                        {
     803                                oSRS.SetPS(dfCenterLat,dfCenterLong,dfProjScaleFactor,
     804                                                        dfProjFalseEast,dfProjFalseNorth);
     805                        }
     806                }
     807
     808                //Export to Wkt.
     809                //In case of error then clean the projection
     810                if (oSRS.exportToWkt(&pszProjection) != OGRERR_NONE
     811                        || oSRS.exportToWkt(&pszGCPProjection) != OGRERR_NONE)
     812                {
     813                        pszProjection = CPLStrdup("");
     814                        pszGCPProjection = CPLStrdup("");
     815                }
     816        }
     817
     818
     819}
     820
     821/**
     822* Get Geotransform information for COSMO-SKYMED files
     823* In case of sucess it stores the transformation
     824* in adfGeoTransform. In case of failure it doesn't
     825* modify adfGeoTransform
     826* @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
     827*/
     828void HDF5ImageDataset::CaptureCSKGeoTransform(int iProductType)
     829{
     830        double      dfInULX, dfInULY, dfInLRX, dfInLRY;
     831        double          dfOutULX, dfOutULY, dfOutLRX, dfOutLRY;
     832
     833        const char *pszAux;
     834        const char *pszAux2;
     835        CPLString osULCoord, osLRCoord;
     836    CPLString osAttributeName;
     837        const char *pszSubdatasetName = GetSubdatasetName();
     838
     839        bHasGeoTransform = FALSE;
     840        //If the product level is not L1C or L1D then
     841        //it doesn't have a valid projection
     842        if(iProductType == PROD_CSK_L1C||iProductType == PROD_CSK_L1D)
     843        {
     844                //If there is a subdataset
     845                if(pszSubdatasetName != NULL)
     846                {
     847                  CPLString osFormatedSubdatasetName =
     848                                        HDF5ReplaceGroupSeparators(pszSubdatasetName);
     849
     850                  osFormatedSubdatasetName += HDF5_INTERAL_GROUP_SEPARATOR;
     851
     852                  osAttributeName = osFormatedSubdatasetName +
     853                                                                                  "Top_Left_East-North";
     854
     855                  //Fetch the coordinates
     856                  pszAux = HDF5Dataset::GetMetadataItem(osAttributeName.c_str());
     857
     858                  osAttributeName = osFormatedSubdatasetName +
     859                                                                                  "Bottom_Right_East-North";
     860
     861                  pszAux2 = HDF5Dataset::GetMetadataItem(osAttributeName.c_str());
     862
     863                        //If it could find the attributes on the metadata
     864                        if( pszAux == NULL || pszAux2 == NULL)
     865                        {
     866                                bHasGeoTransform = FALSE;
     867                        }
     868                        else
     869                        {
     870                                osULCoord = pszAux;
     871                                osLRCoord = pszAux2;
     872
     873                                //Split the array of coordinates
     874                                char ** papszULCoord = CSLTokenizeString2(osULCoord," ",
     875                                                                                                                  CSLT_STRIPLEADSPACES);
     876                                char ** papszLRCoord = CSLTokenizeString2(osLRCoord," ",
     877                                                                                                                  CSLT_STRIPLEADSPACES);
     878
     879                                //Get the pixel-line dimensions (IN)
     880                                dfInULX = 0;
     881                                dfInULY = 0;
     882                                dfInLRX = GetRasterXSize();
     883                                dfInLRY = GetRasterYSize();
     884
     885                                //Convert geographic coordinates to double (OUT)
     886                                dfOutULX = CPLAtof(papszULCoord[0]);
     887                                dfOutULY = CPLAtof(papszULCoord[1]);
     888                                dfOutLRX = CPLAtof(papszLRCoord[0]);
     889                                dfOutLRY = CPLAtof(papszLRCoord[1]);
     890
     891                                //Build the matrix
     892                                adfGeoTransform[1] = dfOutLRX/(dfInLRX - dfInULX) -
     893                                                                                                dfOutULX/(dfInLRX - dfInULX);
     894                                adfGeoTransform[0] = dfOutULX - adfGeoTransform[1] * dfInULX;
     895                                adfGeoTransform[2] = 0;
     896                                adfGeoTransform[5] = dfOutLRY/(dfInLRY - dfInULY) -
     897                                                                                                dfOutULY/(dfInLRY - dfInULY);
     898                                adfGeoTransform[3] = dfOutULY - adfGeoTransform[5] * dfInULY;
     899                                adfGeoTransform[4] = 0;
     900
     901                                CPLFree(papszULCoord);
     902                                CPLFree(papszLRCoord);
     903
     904                                bHasGeoTransform = TRUE;
     905                        }
     906                }
     907        }
     908}
    697909
    698910
  • frmts/hdf5/hdf5imagedataset.h

    diff -rupN /c/Documents and Settings/alexmantaut/Mis documentos/codigo suelto/CSKHDF5/hdf5_clean/hdf5/hdf5imagedataset.h /c/Documents and Settings/alexmantaut/Mis documentos/codigo suelto/CSKHDF5/hdf5/hdf5imagedataset.h
     
     1/******************************************************************************
     2 * $Id: hdf5imagedataset.cpp 21169 2010-11-24 15:57:11Z warmerdam $
     3 *
     4 * Project:  Hierarchical Data Format Release 5 (HDF5)
     5 * Purpose:  Read subdatasets of HDF5 file.
     6 * Author:   Denis Nadeau <denis.nadeau@gmail.com>
     7 *
     8 ******************************************************************************
     9 * Copyright (c) 2005, Frank Warmerdam <warmerdam@pobox.com>
     10 *
     11 * Permission is hereby granted, free of charge, to any person obtaining a
     12 * copy of this software and associated documentation files (the "Software"),
     13 * to deal in the Software without restriction, including without limitation
     14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
     15 * and/or sell copies of the Software, and to permit persons to whom the
     16 * Software is furnished to do so, subject to the following conditions:
     17 *
     18 * The above copyright notice and this permission notice shall be included
     19 * in all copies or substantial portions of the Software.
     20 *
     21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
     22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
     23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
     24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
     25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
     26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
     27 * DEALINGS IN THE SOFTWARE.
     28 ****************************************************************************/
     29
     30#define H5_USE_16_API
     31
     32#ifndef _HDF5IMAGEDATASET_H_INCLUDED_
     33#define _HDF5IMAGEDATASET_H_INCLUDED_
     34
     35#include "hdf5.h"
     36
     37#include "gdal_pam.h"
     38#include "gdal_priv.h"
     39#include "cpl_string.h"
     40#include "hdf5dataset.h"
     41#include "ogr_spatialref.h"
     42
     43CPL_CVSID("$Id: hdf5imagedataset.cpp 21169 2010-11-24 15:57:11Z warmerdam $");
     44
     45CPL_C_START
     46void    GDALRegister_HDF5Image(void);
     47CPL_C_END
     48
     49/* release 1.6.3 or 1.6.4 changed the type of count in some api functions */
     50
     51#if H5_VERS_MAJOR == 1 && H5_VERS_MINOR <= 6 \
     52       && (H5_VERS_MINOR < 6 || H5_VERS_RELEASE < 3)
     53#  define H5OFFSET_TYPE hssize_t
     54#else
     55#  define H5OFFSET_TYPE  hsize_t
     56#endif
     57
     58class HDF5ImageDataset : public HDF5Dataset
     59{
     60        typedef enum
     61        {
     62                UNKNOWN_PRODUCT=0,
     63                CSK_PRODUCT
     64        }Hdf5ProductType;
     65
     66        typedef enum
     67        {
     68            PROD_UNKNOWN=0,
     69            PROD_CSK_L0,
     70            PROD_CSK_L1A,
     71            PROD_CSK_L1B,
     72            PROD_CSK_L1C,
     73            PROD_CSK_L1D
     74        }HDF5CSKProductEnum;
     75
     76    friend class HDF5ImageRasterBand;
     77
     78    char        *pszProjection;
     79    char        *pszGCPProjection;
     80    GDAL_GCP    *pasGCPList;
     81    int         nGCPCount;
     82    OGRSpatialReference oSRS;
     83
     84    hsize_t      *dims,*maxdims;
     85    HDF5GroupObjects *poH5Objects;
     86    int          ndims,dimensions;
     87    hid_t        dataset_id;
     88    hid_t        dataspace_id;
     89    hsize_t      size;
     90    int          address;
     91    hid_t        datatype;
     92    hid_t        native;
     93    H5T_class_t  clas;
     94    int          iSubdatasetType;
     95    double                adfGeoTransform[6];
     96    bool                  bHasGeoTransform;
     97
     98public:
     99    HDF5ImageDataset();
     100    ~HDF5ImageDataset();
     101   
     102    CPLErr CreateProjections( );
     103    static GDALDataset  *Open( GDALOpenInfo * );
     104    static int           Identify( GDALOpenInfo * );
     105
     106    const char          *GetProjectionRef();
     107    virtual int         GetGCPCount( );
     108    virtual const char  *GetGCPProjection();
     109    virtual const GDAL_GCP *GetGCPs( );
     110    virtual CPLErr GetGeoTransform( double * padfTransform );
     111
     112    /**
     113     * Identify if the subdataset has a known product format
     114     * It stores a product identifier in iSubdatasetType,
     115     * UNKNOWN_PRODUCT, if it isn't a recognizable format.
     116     */
     117    void IdentifyProductType();
     118
     119    /**
     120     * Captures Geolocation information from a COSMO-SKYMED
     121     * file.
     122     * The geoid will allways be WGS84
     123     * The projection type may be UTM or UPS, depending on the
     124     * latitude from the center of the image.
     125     * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
     126     */
     127    void CaptureCSKGeolocation(int iProductType);
     128
     129    /**
     130    * Get Geotransform information for COSMO-SKYMED files
     131    * In case of sucess it stores the transformation
     132    * in adfGeoTransform. In case of failure it doesn't
     133    * modify adfGeoTransform
     134    * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
     135    */
     136    void CaptureCSKGeoTransform(int iProductType);
     137   
     138};
     139
     140
     141
     142
     143#endif //_HDF5IMAGEDATASET_H_INCLUDED_