[gdal-dev] SDTS elevation values

[Karthik]

Skipped content of type multipart/alternative-------------- next part --------------
#include <iostream>
#include "sdts_al.h"
#include <stdlib.h>
#include <math.h>
using namespace std;


int east, south, west, north;
int latitudePosts, longitudePosts;
SDTSRasterReader oSDTSRasterReader;
char _szFMT[32];
double xres, yres, tlEasting, tlNorthing, brEasting, brNorthing;
double adfTransform[6];

//Reads a block of raw data from the XXXXCEL0.DDF and interprets
//it as int/float depending on the type of the data.
void getBlock(int row, float elevations[]) {
        int nBytes = oSDTSRasterReader.GetRasterType()==1?2:4;
        char *pabyData = (char*)CPLMalloc(nBytes * longitudePosts);
        oSDTSRasterReader.GetBlock(0, row, pabyData);
        int j=0;
        for(int i=0; i<nBytes * longitudePosts; i+=nBytes) {
                char buffer[nBytes];
                memcpy(buffer, pabyData+i, nBytes);
                if (nBytes == 2) {
                        short *elevation = (short *)buffer;
                        elevations[j++] = (float)*elevation;
                } else if(EQUAL(_szFMT,"BI32") ) {
                        int *elevation = (int *)buffer;
                        elevations[j++] = (float)*elevation;
                } else {
                        float *elevation = (float *)buffer;
                        elevations[j++] = (float)*elevation;
                }
        }
        CPLFree(pabyData);
}

//Helper routine to calculate absolute elevation at a point
double getElevationAtPoint(int col, int row) {
        float elevations[longitudePosts];
        getBlock(row-1, elevations);
        return elevations[col];
}

//The following code's logic follows directly from - http://en.wikipedia.org/wiki/Bilinear_interpolation
float interpolate (int x1, int y2, double x, double y) {

        int x2 = x1 + 1;
        int y1 = y2 + 1;

        cout << "Four surrounding posts are - (" << x1 << "," << y1 << ") " <<
                                             "(" << x2 << "," << y1 << ") " <<
                                             "(" << x1 << "," << y2 << ") " <<
                                             "(" << x2 << "," << y2 << ") " << endl;


        double a00 = getElevationAtPoint(x1, y1);
        double a10 = getElevationAtPoint(x2, y1);
        double a01 = getElevationAtPoint(x1, y2);
        double a11 = getElevationAtPoint(x2, y2);

        cout << "The elevations at these points are : " << a00 << " " << a10 << " " << a01 << " " << a11 << endl;
        double x1_ = tlEasting + xres * x1;
        double y1_ = tlNorthing - yres * y1;
        double x2_ = tlEasting + xres * x2;
        double y2_ = tlNorthing - yres * y2;

        double den = (x2_ - x1_) * (y2_ - y1_);
        double term1 = a00 * (x2_ - x) * (y2_ - y)/den;
        double term2 = a10 * (x - x1_) * (y2_ - y)/den;
        double term3 = a01 * (x2_ - x) * (y - y1_)/den;
        double term4 = a11 * (x - x1_) * (y - y1_)/den;
        return term1 + term2 + term3 + term4;
}
//Get the interpolated elevation at the given UTM easting and northing
float getElevation(double easting, double northing) {
        int row = floor((tlNorthing-northing)/xres);
        int col = floor((easting - tlEasting)/yres);
        return interpolate(col, row, easting, northing);
}

int main() {
        DDFRecord *poRecord;
        //Get the CATD file
        char *pszCATD = (char*)malloc(100);
        cout << "Enter the *CATD file\n";
        cin >> pszCATD;

        //Create the SDTS_CATD and SDTS_IREF objects
        SDTS_CATD *poCATD  = new SDTS_CATD();
        SDTS_IREF *poIREF  = new SDTS_IREF();

        //Read the CATD File
        poCATD -> Read (pszCATD);

        //Read the IREF File
        poIREF -> Read (poCATD -> GetModuleFilePath("IREF"));

        //Read the RSDF File
        DDFModule oRSDF;
        oRSDF.Open(poCATD->GetModuleFilePath("RSDF"));
        poRecord = oRSDF.ReadRecord();
        double dfZ;
        //Get the topleft coordinates
        poIREF -> GetSADR(poRecord->FindField( "SADR" ), 1, adfTransform+0, adfTransform+3, &dfZ);
        tlEasting = adfTransform[0];
        tlNorthing = adfTransform[3];

        //Get the x and y resolutions
        xres = poIREF -> dfXRes;
        yres = poIREF -> dfYRes;

        cout << "X and Y Resolution are : " << xres << "," << yres << endl;

        //Get the format in which the elevations are stored.
        DDFModule oDDSH;
        oDDSH.Open(poCATD->GetModuleFilePath("DDSH"));
        poRecord = oDDSH.ReadRecord();
        strcpy( _szFMT, poRecord->GetStringSubfield("DDSH",0,"FMT",0) );


        DDFModule oLDEF;
        oLDEF.Open(poCATD->GetModuleFilePath("LDEF"));
        poRecord = oLDEF.ReadRecord();
        //Get the number of rows and columns (latitude and longitude posts)
		longitudePosts = poRecord->GetIntSubfield( "LDEF", 0, "NCOL", 0 );
        latitudePosts = poRecord->GetIntSubfield( "LDEF", 0, "NROW", 0 );

        char szINTR[5];
        strcpy(szINTR, poRecord -> GetStringSubfield(  "LDEF", 0, "INTR", 0 ) );
        //if the TL coordinates are actually the center pixel, shift them to left and
        //top by half pixel
        if (EQUAL(szINTR, "CE")) {
                tlEasting -= xres * 0.5;
                tlNorthing -= yres * 0.5;
        }
        cout << "Top Left UTM Coordinates are : " << tlEasting << "," << tlNorthing << endl;

        //Instantiate the Raster Reader for reading a block of elevation values from the
        //CEL0 Module.
        oSDTSRasterReader.Open(poCATD, poIREF, "CEL0");

        cout << "Enter the point at which elevation is needed" << endl;
        double _easting, _northing;
        //Read the points where the elevation is needed
        cin >> _easting >> _northing;

        //calculate the interpolated elevation at this point.
        cout << "Interpolated Elevation : " << getElevation(_easting, _northing) << endl;

}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: makefile
Type: application/octet-stream
Size: 176 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/gdal-dev/attachments/20080803/5275b4a6/makefile.obj