/****************************************************************************** * levellerdataset.cpp,v 1.0 * * Project: Leveller TER Driver * Purpose: Reader for Leveller TER documents * Author: Ray Gardener, Daylon Graphics Ltd. * * Portions of this module derived from GDAL drivers by * Frank Warmerdam, see http://www.gdal.org * ****************************************************************************** * Copyright (c) 2005 Daylon Graphics Ltd. * * 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: levellerdataset.cpp,v $ * Revision 1.4 2006/02/02 22:30:49 fwarmerdam * byte swapping fix for mac * * Revision 1.3 2005/10/21 00:03:41 fwarmerdam * added coordinate system support from Ray * * Revision 1.2 2005/10/20 20:18:14 fwarmerdam * Applied some patches from Ray. * * Revision 1.1 2005/10/20 13:44:29 fwarmerdam * New * */ #include "gdal_pam.h" #include "ogr_spatialref.h" CPL_CVSID("$Id: levellerdataset.cpp,v 1.4 2006/02/02 22:30:49 fwarmerdam Exp $"); CPL_C_START void GDALRegister_Leveller(void); CPL_C_END #define str_equal(_s1, _s2) (0 == strcmp((_s1),(_s2))) #define array_size(_a) (sizeof(_a) / sizeof(_a[0])) /*GDALDataset *LevellerCreateCopy( const char *, GDALDataset *, int, char **, GDALProgressFunc pfnProgress, void * pProgressData ); */ /************************************************************************/ /* ==================================================================== */ /* LevellerDataset */ /* ==================================================================== */ /************************************************************************/ class LevellerRasterBand; class LevellerDataset : public GDALPamDataset { friend class LevellerRasterBand; int m_version; char* m_pszProjection; char m_szWorldscaleUnits[32]; double m_dWorldscale; double m_dElevScale; double m_dElevBase; double m_adfTransform[6]; FILE* m_fp; vsi_l_offset m_nDataOffset; int LoadFromFile(FILE*); int locate_data(vsi_l_offset&, size_t&, FILE*, const char*); int get(int&, FILE*, const char*); int get(double&, FILE*, const char*); int get(char*, size_t, FILE*, const char*); double convert_measure(double, const char* pszUnitsFrom); public: LevellerDataset(); ~LevellerDataset(); static GDALDataset* Open( GDALOpenInfo* ); static GDALDataset* Create( const char* pszFilename, int nXSize, int nYSize, int nBands, GDALDataType eType, char** papszOptions ); virtual CPLErr GetGeoTransform( double* ); virtual const char* GetProjectionRef(void); }; /************************************************************************/ /* ==================================================================== */ /* LevellerRasterBand */ /* ==================================================================== */ /************************************************************************/ class LevellerRasterBand : public GDALPamRasterBand { friend class LevellerDataset; public: LevellerRasterBand(LevellerDataset*); // Geomeasure support. virtual const char* GetUnitType(); virtual double GetScale(int* pbSuccess = NULL); virtual double GetOffset(int* pbSuccess = NULL); virtual CPLErr IReadBlock( int, int, void * ); }; /************************************************************************/ /* LevellerRasterBand() */ /************************************************************************/ LevellerRasterBand::LevellerRasterBand( LevellerDataset *poDS ) { this->poDS = poDS; this->nBand = 1; eDataType = GDT_Float32; nBlockXSize = poDS->GetRasterXSize(); nBlockYSize = 1;//poDS->GetRasterYSize(); } /************************************************************************/ /* IReadBlock() */ /************************************************************************/ CPLErr LevellerRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, void* pImage ) { CPLAssert( sizeof(float) == sizeof(GInt32) ); CPLAssert( nBlockXOff == 0 ); CPLAssert( pImage != NULL ); LevellerDataset *poGDS = (LevellerDataset *) poDS; /* -------------------------------------------------------------------- */ /* Seek to scanline. */ /* -------------------------------------------------------------------- */ const size_t rowbytes = nBlockXSize * sizeof(float); if(0 != VSIFSeekL( poGDS->m_fp, poGDS->m_nDataOffset + nBlockYOff * rowbytes, SEEK_SET)) { CPLError( CE_Failure, CPLE_FileIO, ".bt Seek failed:%s", VSIStrerror( errno ) ); return CE_Failure; } /* -------------------------------------------------------------------- */ /* Read the scanline into the image buffer. */ /* -------------------------------------------------------------------- */ if( VSIFReadL( pImage, rowbytes, 1, poGDS->m_fp ) != 1 ) { CPLError( CE_Failure, CPLE_FileIO, "Leveller read failed:%s", VSIStrerror( errno ) ); return CE_Failure; } /* -------------------------------------------------------------------- */ /* Swap on MSB platforms. */ /* -------------------------------------------------------------------- */ #ifdef CPL_MSB GDALSwapWords( pImage, 4, nRasterXSize, 4 ); #endif /* -------------------------------------------------------------------- */ /* Convert from legacy-format fixed-point if necessary. */ /* -------------------------------------------------------------------- */ float* pf = (float*)pImage; if(poGDS->m_version < 6) { GInt32* pi = (int*)pImage; for(int i = 0; i < nBlockXSize; i++) pf[i] = (float)pi[i] / 65536; } /* -------------------------------------------------------------------- */ /* Convert raw elevations to realworld elevs. */ /* -------------------------------------------------------------------- */ for(int i = 0; i < nBlockXSize; i++) pf[i] *= poGDS->m_dWorldscale; //this->GetScale(); return CE_None; } /************************************************************************/ /* GetUnitType() */ /************************************************************************/ const char *LevellerRasterBand::GetUnitType() { // Return elevation units. // For Leveller documents, it's the same as the ground units. LevellerDataset *poGDS = (LevellerDataset *) poDS; return poGDS->m_szWorldscaleUnits; } /************************************************************************/ /* GetScale() */ /************************************************************************/ double LevellerRasterBand::GetScale(int* pbSuccess) { LevellerDataset *poGDS = (LevellerDataset *) poDS; if(pbSuccess != NULL) *pbSuccess = TRUE; return poGDS->m_dElevScale; } /************************************************************************/ /* GetOffset() */ /************************************************************************/ double LevellerRasterBand::GetOffset(int* pbSuccess) { LevellerDataset *poGDS = (LevellerDataset *) poDS; if(pbSuccess != NULL) *pbSuccess = TRUE; return poGDS->m_dElevBase; } /************************************************************************/ /* ==================================================================== */ /* LevellerDataset */ /* ==================================================================== */ /************************************************************************/ /************************************************************************/ /* LevellerDataset() */ /************************************************************************/ LevellerDataset::LevellerDataset() { m_fp = NULL; m_pszProjection = NULL; } /************************************************************************/ /* ~LevellerDataset() */ /************************************************************************/ LevellerDataset::~LevellerDataset() { FlushCache(); CPLFree(m_pszProjection); if( m_fp != NULL ) VSIFCloseL( m_fp ); } /************************************************************************/ /* LoadFromFile() */ /* */ /* If the data from DEM is in meters, then values are stored as */ /* shorts. If DEM data is in feet, then height data will be */ /* stored in float, to preserve the precision of the original */ /* data. returns true if the file was successfully opened and */ /* read. */ /************************************************************************/ int LevellerDataset::locate_data(vsi_l_offset& offset, size_t& len, FILE* fp, const char* pszTag) { // Locate the file offset of the desired tag's data. // If it is not available, return false. // If the tag is found, leave the filemark at the // start of its data. if(0 != VSIFSeekL(fp, 5, SEEK_SET)) return 0; const int kMaxDescLen = 64; for(;;) { unsigned char c; if(1 != VSIFReadL(&c, sizeof(c), 1, fp)) return 0; const int descriptorLen = c; if(descriptorLen == 0 || descriptorLen > kMaxDescLen) return 0; char descriptor[kMaxDescLen+1]; if(1 != VSIFReadL(descriptor, descriptorLen, 1, fp)) return 0; GInt32 datalen; if(1 != VSIFReadL(&datalen, sizeof(datalen), 1, fp)) return 0; datalen = CPL_LSBWORD32(datalen); descriptor[descriptorLen] = 0; if(str_equal(descriptor, pszTag)) { len = (size_t)datalen; offset = VSIFTellL(fp); return 1; } else { // Seek to next tag. if(0 != VSIFSeekL(fp, (vsi_l_offset)datalen, SEEK_CUR)) return 0; } } return 0; } /************************************************************************/ /* get() */ /************************************************************************/ int LevellerDataset::get(int& n, FILE* fp, const char* psz) { vsi_l_offset offset; size_t len; if(this->locate_data(offset, len, fp, psz)) { GInt32 value; if(1 == VSIFReadL(&value, sizeof(value), 1, fp)) { CPL_LSBPTR32(&value); n = (int)value; return 1; } } return 0; } /************************************************************************/ /* get() */ /************************************************************************/ int LevellerDataset::get(double& d, FILE* fp, const char* pszTag) { vsi_l_offset offset; size_t len; if(this->locate_data(offset, len, fp, pszTag)) { if(1 == VSIFReadL(&d, sizeof(d), 1, fp)) { CPL_LSBPTR64(&d); return 1; } } return 0; } /************************************************************************/ /* get() */ /************************************************************************/ int LevellerDataset::get(char* pszValue, size_t maxchars, FILE* fp, const char* pszTag) { char szTag[65]; // We can assume 8-bit encoding, so just go straight // to the *_d tag. sprintf(szTag, "%sd", pszTag); vsi_l_offset offset; size_t len; if(this->locate_data(offset, len, fp, szTag)) { if(len > maxchars) return 0; if(1 == VSIFReadL(pszValue, len, 1, fp)) { pszValue[len] = 0; // terminate C-string return 1; } } return 0; } /************************************************************************/ /* convert_measure() */ /************************************************************************/ double LevellerDataset::convert_measure ( double d, const char* pszSpace ) { // Convert a measure to meters. const char* szLabels[] = { "fm", "pm", "A", "nm", "u", "um", "mm", "cm", "in", "dft", "dm", "sp", "ft", "sft", "yd", "m", "r", "dam", "dkm", "f", "km", "mi", "nmi", "ls", "lm", "AU", "lhr", "ld", "ly", "pc", "kly" }; const double dLYtoM = 9.4608953536e+15; const double dScales[] = { 1.0e-15, // fm (femtometer) 1.0e-12, // pm (picometer) 1.0e-10, // A (angstrom) 1.0e-9, // nm (nanometer) 1.0e-6, // u (micron) 1.0e-6, // um (micrometer) 1.0e-3, // mm (millimeter) 1.0e-2, // cm (centimeter) 0.0254, // in (inch) 0.03048, // dft (decifoot) 1.0e-1, // dm (decimeter) 0.2286, // sp (span) 0.3048, // ft (foot) 1200.0 / 3937.0, // sft (survey foot) 0.9144, // yd (yard) 1.0, // m (meter) 5.029, // r (rod) 10.0, // dam (decameter) 10.0, // dkm (decameter) 201.168, // f (furlong) 1.0e+3, // km 1609.344,// mi (mile) //1611.7874086021505376344086021505, 1853.0, // nmi (nautical mile) dLYtoM / 3.16E+07, // ls (light second) // 299,367,088.60759493670886075949367 meters. // 299,791,819.008 if we use 186,282 miles per ls. dLYtoM / 5.26E+05, // lm (light minute) 1.50e+11, // AU (astronomical units) dLYtoM / 5.26E+05 * 60 , // lhr (light hour) dLYtoM / 5.26E+05 * 60 * 24 , // ld (light day) dLYtoM, // ly (light year) dLYtoM * 3.26, // 3.08e+16, // pc (parsec) dLYtoM * 1000, // kly (kilo light year) }; CPLAssert(array_size(dScales) == array_size(szLabels)); for(size_t i = 0; i < array_size(dScales); i++) { if(str_equal(pszSpace, szLabels[i])) return d * dScales[i]; } CPLAssert(0); return d; } /************************************************************************/ /* LoadFromFile() */ /************************************************************************/ int LevellerDataset::LoadFromFile(FILE* file) { // get hf dimensions if(!this->get(nRasterXSize, file, "hf_w")) return 0; if(!this->get(nRasterYSize, file, "hf_b")) return 0; // record start of pixel data size_t datalen; if(!this->locate_data(m_nDataOffset, datalen, file, "hf_data")) return 0; // sanity check: do we have enough pixels? if(datalen != nRasterXSize * nRasterYSize * sizeof(float)) return 0; // Read any world scaling. m_dWorldscale = 1.0; strcpy(m_szWorldscaleUnits, "m"); m_dElevBase = 0.0; m_adfTransform[0] = 0.0; m_adfTransform[1] = 1.0; m_adfTransform[2] = 0.0; m_adfTransform[3] = 0.0; m_adfTransform[4] = 0.0; m_adfTransform[5] = 1.0; if(this->get(m_dWorldscale, file, "hf_worldspacing")) { //m_bHasWorldscale = true; if(this->get(m_szWorldscaleUnits, sizeof(m_szWorldscaleUnits)-1, file, "hf_worldspacinglabel")) { // Drop long name, if present. char* p = strchr(m_szWorldscaleUnits, ' '); if(p != NULL) *p = 0; } // If the units are something besides m/ft/sft, // then convert them to meters. if(!str_equal("m", m_szWorldscaleUnits) && !str_equal("ft", m_szWorldscaleUnits) && !str_equal("sft", m_szWorldscaleUnits)) { m_dWorldscale = this->convert_measure(m_dWorldscale, m_szWorldscaleUnits); strcpy(m_szWorldscaleUnits, "m"); } // Our extents are such that the origin is at the // center of the heightfield. m_adfTransform[0] = -0.5 * m_dWorldscale * nRasterXSize; m_adfTransform[3] = -0.5 * m_dWorldscale * nRasterYSize; m_adfTransform[1] = m_dWorldscale; m_adfTransform[5] = m_dWorldscale; } m_dElevScale = 1.0;//m_dWorldscale; /* -------------------------------------------------------------------- */ /* Set projection. */ /* -------------------------------------------------------------------- */ // Leveller files as of Oct 2005 are not currently georeferenced, // but we can't indicate ground units without a spatialref, // so make a default geocs (local). OGRSpatialReference sr; sr.SetLocalCS("Leveller world space"); if(OGRERR_NONE != sr.SetLinearUnits(m_szWorldscaleUnits, this->convert_measure(1.0, m_szWorldscaleUnits))) return 0; if(OGRERR_NONE != sr.exportToWkt(&m_pszProjection)) return 0; return TRUE; } /************************************************************************/ /* GetProjectionRef() */ /************************************************************************/ const char* LevellerDataset::GetProjectionRef(void) { if(m_pszProjection == NULL ) return ""; else return m_pszProjection; } /************************************************************************/ /* GetGeoTransform() */ /************************************************************************/ CPLErr LevellerDataset::GetGeoTransform(double* padfTransform) { // Return identity transform, since Leveller heightfields // do not currently have geodata. memcpy(padfTransform, m_adfTransform, sizeof(m_adfTransform)); return CE_None; } /************************************************************************/ /* Open() */ /************************************************************************/ GDALDataset *LevellerDataset::Open( GDALOpenInfo * poOpenInfo ) { // The file should have at least 5 header bytes // and hf_w, hf_b, and hf_data tags. if( poOpenInfo->nHeaderBytes < 5+13+13+16 ) return NULL; if( !EQUALN((const char *) poOpenInfo->pabyHeader, "trrn",4) ) return NULL; const int version = poOpenInfo->pabyHeader[4]; if(version < 4 || version > 6) return NULL; /* -------------------------------------------------------------------- */ /* Create a corresponding GDALDataset. */ /* -------------------------------------------------------------------- */ LevellerDataset *poDS; poDS = new LevellerDataset(); poDS->m_version = version; // Reopen for large file access. if( poOpenInfo->eAccess == GA_Update ) poDS->m_fp = VSIFOpenL( poOpenInfo->pszFilename, "rb+" ); else poDS->m_fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" ); if( poDS->m_fp == NULL ) { CPLError( CE_Failure, CPLE_OpenFailed, "Failed to re-open %s within Leveller driver.\n", poOpenInfo->pszFilename ); return NULL; } poDS->eAccess = poOpenInfo->eAccess; /* -------------------------------------------------------------------- */ /* Read the file. */ /* -------------------------------------------------------------------- */ if( !poDS->LoadFromFile( poDS->m_fp ) ) { delete poDS; return NULL; } /* -------------------------------------------------------------------- */ /* Create band information objects. */ /* -------------------------------------------------------------------- */ poDS->SetBand( 1, new LevellerRasterBand( poDS )); poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT ); /* -------------------------------------------------------------------- */ /* Initialize any PAM information. */ /* -------------------------------------------------------------------- */ poDS->SetDescription( poOpenInfo->pszFilename ); poDS->TryLoadXML(); return( poDS ); } /************************************************************************/ /* GDALRegister_Leveller() */ /************************************************************************/ void GDALRegister_Leveller() { GDALDriver *poDriver; if( GDALGetDriverByName( "Leveller" ) == NULL ) { poDriver = new GDALDriver(); poDriver->SetDescription( "Leveller" ); poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "ter" ); poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "Leveller heightfield" ); poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_leveller.html" ); poDriver->pfnOpen = LevellerDataset::Open; GetGDALDriverManager()->RegisterDriver( poDriver ); } }