ntf_raster.cpp

GIS系统支持库Geospatial Data Abstraction Library代码.GDAL is a translator library for raster geospatial dat

/******************************************************************************
 * $Id: ntf_raster.cpp,v 1.13 2003/07/08 15:23:58 warmerda Exp $
 *
 * Project:  NTF Translator
 * Purpose:  Handle UK Ordnance Survey Raster DTM products.  Includes some
 *           raster related methods from NTFFileReader and the implementation
 *           of OGRNTFRasterLayer.
 * Author:   Frank Warmerdam, [e-mail]
 *
 ******************************************************************************
 * Copyright (c) 1999, Frank Warmerdam
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ******************************************************************************
 *
 * $Log: ntf_raster.cpp,v $
 * Revision 1.13  2003/07/08 15:23:58  warmerda
 * avoid casting warning
 *
 * Revision 1.12  2003/04/14 14:14:40  warmerda
 * fixed problem with Y origin for Landranger DTM
 *
 * Revision 1.11  2003/02/27 21:09:47  warmerda
 * Apply ZMult to LANDFORM_PROFILE_DTM data to scale to meters.
 * adfGeoTransform[5] is left positive to correct coordinates ... the first
 * entry in the first profile column is the *south* west corner, not the
 * north west corner.  See also: FME bug RT3887.
 * Note: adfGeoTransform[] change applied to LANDRANGER and LANDFORM products
 * but only check on LANDFORM products.
 *
 * Revision 1.10  2001/07/18 04:55:16  warmerda
 * added CPL_CSVID
 *
 * Revision 1.9  2001/01/19 20:31:12  warmerda
 * expand tabs
 *
 * Revision 1.8  1999/11/03 19:27:21  warmerda
 * Pass column in SetFPPos() as the feature id to ensure that the seek
 * isn't arbitrarily ignored.
 *
 * Revision 1.7  1999/11/03 19:07:58  warmerda
 * open raster file if not open when trying to read
 *
 * Revision 1.6  1999/10/04 18:19:22  warmerda
 * Avoid type conversion warnings.
 *
 * Revision 1.5  1999/10/04 13:38:39  warmerda
 * Fixed handling for filling holes in the column offset array.
 *
 * Revision 1.4  1999/10/04 13:28:43  warmerda
 * added DEM_SAMPLE support
 *
 * Revision 1.3  1999/10/04 12:52:22  warmerda
 * Added DTM_ as a prefix to raster layer names.
 *
 * Revision 1.2  1999/10/04 11:36:53  warmerda
 * Added HEIGHT attribute.
 *
 * Revision 1.1  1999/10/04 03:07:49  warmerda
 * New
 *
 */

#include "ntf.h"

CPL_CVSID("$Id: ntf_raster.cpp,v 1.13 2003/07/08 15:23:58 warmerda Exp $");

/************************************************************************/
/* ==================================================================== */
/*                     NTFFileReader Raster Methods                     */
/* ==================================================================== */
/************************************************************************/

/************************************************************************/
/*                          IsRasterProduct()                           */
/************************************************************************/

int NTFFileReader::IsRasterProduct()

{
    return GetProductId() == NPC_LANDRANGER_DTM
        || GetProductId() == NPC_LANDFORM_PROFILE_DTM;
}

/************************************************************************/
/*                       EstablishRasterAccess()                        */
/************************************************************************/

void NTFFileReader::EstablishRasterAccess()

{
/* -------------------------------------------------------------------- */
/*      Read the type 50 record.                                        */
/* -------------------------------------------------------------------- */
    NTFRecord   *poRecord;

    while( (poRecord = ReadRecord()) != NULL
           && poRecord->GetType() != NRT_GRIDHREC
           && poRecord->GetType() != NRT_VTR )
    {
        delete poRecord;
    }

    if( poRecord->GetType() != NRT_GRIDHREC )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Unable to find GRIDHREC (type 50) record in what appears\n"
                  "to be an NTF Raster DTM product." );
        return;
    }

/* -------------------------------------------------------------------- */
/*      Parse if LANDRANGER_DTM                                         */
/* -------------------------------------------------------------------- */
    if( GetProductId() == NPC_LANDRANGER_DTM )
    {
        nRasterXSize = atoi(poRecord->GetField(13,16));
        nRasterYSize = atoi(poRecord->GetField(17,20));

        // NOTE: unusual use of GeoTransform - the pixel origin is the
        // bottom left corner!
        adfGeoTransform[0] = atoi(poRecord->GetField(25,34));
        adfGeoTransform[1] = 50;
        adfGeoTransform[2] = 0;
        adfGeoTransform[3] = atoi(poRecord->GetField(35,44));
        adfGeoTransform[4] = 0;
        adfGeoTransform[5] = 50;
        
        nRasterDataType = 3; /* GDT_Int16 */
    }

/* -------------------------------------------------------------------- */
/*      Parse if LANDFORM_PROFILE_DTM                                   */
/* -------------------------------------------------------------------- */
    else if( GetProductId() == NPC_LANDFORM_PROFILE_DTM )
    {
        nRasterXSize = atoi(poRecord->GetField(23,30));
        nRasterYSize = atoi(poRecord->GetField(31,38));

        // NOTE: unusual use of GeoTransform - the pixel origin is the
        // bottom left corner!
        adfGeoTransform[0] = atoi(poRecord->GetField(13,17))
                           + GetXOrigin();
        adfGeoTransform[1] = atoi(poRecord->GetField(39,42));
        adfGeoTransform[2] = 0;
        adfGeoTransform[3] = atoi(poRecord->GetField(18,22))
                           + GetYOrigin();
        adfGeoTransform[4] = 0;
        adfGeoTransform[5] = atoi(poRecord->GetField(43,46));
        
        nRasterDataType = 3; /* GDT_Int16 */
    }

/* -------------------------------------------------------------------- */
/*      Initialize column offsets table.                                */
/* -------------------------------------------------------------------- */
    panColumnOffset = (long *) CPLCalloc(sizeof(long),nRasterXSize);

    GetFPPos( panColumnOffset+0, NULL );

/* -------------------------------------------------------------------- */
/*      Create an OGRSFLayer for this file readers raster points.       */
/* -------------------------------------------------------------------- */
    if( poDS != NULL )
    {
        poRasterLayer = new OGRNTFRasterLayer( poDS, this );
        poDS->AddLayer( poRasterLayer );
    }
}

/************************************************************************/
/*                          ReadRasterColumn()                          */
/************************************************************************/

CPLErr NTFFileReader::ReadRasterColumn( int iColumn, float *pafElev )

{
/* -------------------------------------------------------------------- */
/*      If we don't already have the scanline offset of the previous    */
/*      line, force reading of previous records to establish it.        */
/* -------------------------------------------------------------------- */
    if( panColumnOffset[iColumn] == 0 )
    {
        int     iPrev;
        
        for( iPrev = 0; iPrev < iColumn-1; iPrev++ )
        {
            if( panColumnOffset[iPrev+1] == 0 )
            {
                CPLErr  eErr;
                
                eErr = ReadRasterColumn( iPrev, NULL );
                if( eErr != CE_None )
                    return eErr;
            }
        }
    }

/* -------------------------------------------------------------------- */
/*      If the dataset isn't open, open it now.                         */
/* -------------------------------------------------------------------- */
    if( GetFP() == NULL )
        Open();
    
/* -------------------------------------------------------------------- */
/*      Read requested record.                                          */
/* -------------------------------------------------------------------- */
    NTFRecord   *poRecord;
    
    SetFPPos( panColumnOffset[iColumn], iColumn );
    poRecord = ReadRecord();

    if( iColumn < nRasterXSize-1 )
    {
        GetFPPos( panColumnOffset+iColumn+1, NULL );
    }
    
/* -------------------------------------------------------------------- */
/*      Handle LANDRANGER DTM columns.                                  */
/* -------------------------------------------------------------------- */
    if( pafElev != NULL && GetProductId() == NPC_LANDRANGER_DTM )
    {
        double  dfVScale, dfVOffset;

        dfVOffset = atoi(poRecord->GetField(56,65));
        dfVScale = atoi(poRecord->GetField(66,75)) * 0.001;

        for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )