[Gdal-dev] Parse DEM file problem
Chapman, Martin
MChapman at sanz.com
Wed May 5 15:21:57 EDT 2004
Bruce,
I use GDAL for elevation files and have not seen the problems you are
experiencing. I save a pointer to the dataset when I initialize my
class and call it's methods later in different methods with no problem.
However, when I use the dataset to read info from the file I always copy
the info to my own structures before I operate on the data. My code is
listed below:
In OnOpenDocument:
BOOL CImagePortDoc::OnOpenDocument(LPCTSTR lpszPathName)
{
try
{
if (!CDocument::OnOpenDocument(lpszPathName))
return FALSE;
m_sFileName = lpszPathName;
if (m_sFileName.size() == 0)
return false;
bool bUseOverView = false;
string ovr_file;
string ovr_name;
string file = m_sFileName;
char *pdest = strstr( file.c_str(), ".ovr" );
if (pdest != NULL)
bUseOverView = true;
string::size_type pos = file.find('.');
if (pos == string::npos)
{
ovr_file = file + ".ovr";
ovr_name = file;
}
else
{
ovr_file.insert(0, file.begin(), &file[pos]);
ovr_name = ovr_file;
ovr_file += ".ovr";
}
GDALAllRegister();
if (!bUseOverView)
{
CMainFrame* pMainFrame = (CMainFrame*)
AfxGetMainWnd();
if (pMainFrame->FileExists(ovr_file))
{
CString sMessage = _T("Viewport has
detected an associated\noverview file for ");
sMessage += m_sFileName.c_str();
sMessage += _T(" .\nDo you want to use
the overviews?");
if (AfxMessageBox(sMessage,
MB_YESNO|MB_ICONQUESTION) == IDYES)
{
bUseOverView = true;
}
}
}
m_pDataset = (GDALDataset*)
GDALOpen(m_sFileName.c_str(), GA_ReadOnly);
if( m_pDataset == NULL )
{
const char* pszError = CPLGetLastErrorMsg();
AfxMessageBox(pszError);
return false;
}
m_pGDALDefaultOverviews = new GDALDefaultOverviews();
int nOverviewCount = 0;
if (bUseOverView)
{
if (!m_pGDALDefaultOverviews->IsInitialized())
m_pGDALDefaultOverviews->Initialize(m_pDataset, ovr_file.c_str(), 1);
nOverviewCount =
m_pGDALDefaultOverviews->GetOverviewCount(1);
}
for (int k = 0; k < nOverviewCount; k++)
m_vResLevels.push_back(k);
int nNumberOfBands = m_pDataset->GetRasterCount();
for (int j = 0; j < nNumberOfBands; j++)
m_vBandList.push_back(j);
CString sFileExtension;
CString sFileName = m_sFileName.c_str();
int n = 0;
int nPos = sFileName.ReverseFind('.');
if (nPos != 1)
{
n = sFileName.GetLength() - nPos;
sFileExtension = sFileName.Right(n - 1);
}
else
{
AfxMessageBox(_T("Could not determine file
extension.\nPlease make sure the file has a valid extension."));
return NULL;
}
sFileExtension.MakeLower();
m_nFileType = GetFileType((string) sFileExtension);
m_DataType =
m_pDataset->GetRasterBand(1)->GetRasterDataType();
return TRUE;
}
catch(...)
{
//throw CString("An unspecified error occurred");
}
}
In another method I get the bytes for each band and then bit interleave
them.
NOTE that m_nFileType == 3 is elevation data:
CDIBSection* CImagePortDoc::GetAreaOfInterest(CRect rect, int nResLevel)
{
try
{
int nBandCount = m_vBandList.size();
int nWidth = rect.Width();
int nHeight = rect.Height();
int nXOffset = rect.left;
int nYOffset = rect.top;
void * ptr_dest = NULL;
void * ptr_src = NULL;
int xsize = 0;
int ysize = 0;
m_nResLevel = nResLevel;
int nSize = nBandCount;
void** pBandBuffer = NULL;
if (m_nFileType == 2)
pBandBuffer = (void**) new unsigned
char*[nSize];
else
{
if (m_DataType == GDT_Int16 || m_DataType ==
GDT_Byte)
pBandBuffer = (void**) new
short*[nSize];
else
pBandBuffer = (void**) new
float*[nSize];
}
GDALColorTable* pCT =
m_pDataset->GetRasterBand(1)->GetColorTable();
GDALRasterBand* pBand = NULL;
void* pBandBytes = NULL;
for (int n = 0; n < nSize; n++)
{
if (nBandCount > 1)
{
if
(m_pGDALDefaultOverviews->GetOverviewCount(nSize) > 0 && m_nResLevel >
0)
pBand =
m_pGDALDefaultOverviews->GetOverview(m_vBandList[n] + 1, m_nResLevel -
1);
else
pBand =
m_pDataset->GetRasterBand(m_vBandList[n] + 1);
}
else
{
if
(m_pGDALDefaultOverviews->GetOverviewCount(1) > 0 && m_nResLevel > 0)
pBand =
m_pGDALDefaultOverviews->GetOverview(m_vBandList[0] + 1, m_nResLevel -
1);
else
pBand =
m_pDataset->GetRasterBand(1);
}
int nRasterX = pBand->GetXSize();
int nRasterY = pBand->GetYSize();
if (nWidth > nRasterX)
{
nXOffset = 0;
nWidth = nRasterX;
}
if (nHeight > nRasterY)
{
nYOffset = 0;
nHeight = nRasterY;
}
if (m_nFileType == 2)
{
pBandBytes = (unsigned char*)
CPLMalloc(sizeof(unsigned char) * nWidth * nHeight);
pBand->RasterIO(GF_Read, nXOffset,
nYOffset, nWidth, nHeight, pBandBytes, nWidth , nHeight, GDT_Byte, 0,
0);
}
else
{
if (m_DataType == GDT_Int16 ||
m_DataType == GDT_Byte)
{
pBandBytes = (short*)
CPLMalloc(sizeof(short) * nWidth * nHeight);
pBand->RasterIO(GF_Read,
nXOffset, nYOffset, nWidth, nHeight, pBandBytes, nWidth , nHeight,
GDT_Int16, 0, 0);
}
else
{
pBandBytes = (float*)
CPLMalloc(sizeof(float) * nWidth * nHeight);
pBand->RasterIO(GF_Read,
nXOffset, nYOffset, nWidth, nHeight, pBandBytes, nWidth , nHeight,
GDT_Float32, 0, 0);
}
}
pBandBuffer[n] = pBandBytes;
}
unsigned char* pCombinedBytes = (unsigned char*)
CPLMalloc(sizeof(unsigned char) * nWidth * nHeight * nSize);
for (long i = 0, long j = 0; i < (nWidth * nHeight *
nSize); i += nSize, j++)
{
int k = i;
if (m_nFileType == 2)
{
for (int m = nSize - 1; m >= 0; m--,
k++)
{
unsigned char* pBytes =
(unsigned char*) pBandBuffer[m];
ptr_dest = pCombinedBytes + k;
unsigned char b = pBytes[j];
if (pCT != NULL)
{
GDALColorEntry ce;
unsigned int c =
(unsigned int) b;
c =
pCT->GetColorEntryAsRGB(c, &ce);
if (m == 0) c = ce.c1;
if (m == 1) c = ce.c2;
if (m == 2) c = ce.c3;
b = (unsigned char) c;
}
ptr_src = &b;
memcpy(ptr_dest, ptr_src, 1);
}
}
else
{
if (m_DataType == GDT_Int16 ||
m_DataType == GDT_Byte)
{
for (int m = nSize - 1; m >= 0;
m--, k++)
{
short* pBytes = (short*)
pBandBuffer[m];
ptr_dest =
pCombinedBytes + k;
short nVal = pBytes[j];
unsigned char uc = (nVal
& 0xffff) >> 4;
ptr_src = &uc;
memcpy(ptr_dest,
ptr_src, 1);
}
}
else
{
for (int m = nSize - 1; m >= 0;
m--, k++)
{
float* pBytes = (float*)
pBandBuffer[m];
ptr_dest =
pCombinedBytes + k;
float nVal = pBytes[j];
unsigned char uc =
(unsigned char) nVal >> 8;
ptr_src = &uc;
memcpy(ptr_dest,
ptr_src, 1);
}
}
}
}
if ((nWidth > 0) && (nHeight > 0))
{
if (m_pDIBSection) delete m_pDIBSection;
m_pDIBSection = new CDIBSection;
if (m_pDIBSection)
{
m_pDIBSection->Create(nWidth ,nHeight,
24);
unsigned int dibwidth =
m_pDIBSection->GetTotalWidth();
GByte* dest = (GByte*)
m_pDIBSection->GetBits();
GByte* src = (GByte*) pCombinedBytes;
for (int row = 0; row < nHeight; row++)
{
ptr_dest = dest + row * dibwidth
* nSize;
ptr_src = src + row * nWidth *
nSize;
memcpy(ptr_dest, ptr_src, nWidth
* nSize);
}
}
}
if (pBandBuffer)
{
for (int d = 0; d < nSize; d++)
{
CPLFree(pBandBuffer[d]);
}
delete [] pBandBuffer;
}
if (pCombinedBytes)
CPLFree(pCombinedBytes);
CDIBSection* pDIBSection = m_pDIBSection;
return pDIBSection;
}
catch(...)
{
AfxMessageBox("An unspecified error occurred");
}
}
Maybe you can see something in the code above that will help you resolve
your problem.
Martin
-----Original Message-----
From: Clay, Bruce [mailto:bclay at ball.com]
Sent: Wednesday, May 05, 2004 12:45 PM
To: gdal-dev at remotesensing.org
Subject: [Gdal-dev] Parse DEM file problem
I am working on a basic terrain elevation server. During startup I read
the header of all of the DEM files in the specified directory but do not
read the actual data. The GDALDataset pointer returned by GDALOpen is
saved in an STL vector along with the header information.
I am using the GTOPO30 data set.
When I use the data set pointer later to read file from E20N40.DEM
GDALComputeRasterMinMax returns min = -405 max = 4099 which should be
min = -407 max = 4525.
When I use the data set pointer later to read file from E20N90.DEM I get
min = -32768 and max = 32521 which should be min = -137 max = 5483.
Walking through the data seems to give the same min and max values (at
least for what I have done so far).
I see plenty of -9999 which the documentation says is the "no data"
marker so it doesn't seem like it is a byte swap issue.
The questions are:
Does GDAL maintain any state information that I might be overwriting
by reading other headers before reading the data?
Does anyone have any experience with this data set to know if what to
expect from GDAL when it reads the data?
The code I am using comes right out of the test programs / examples and
is as follows:
pBand = pDataSet->GetRasterBand( 1 );
if ((pBand != NULL) && (pHeaderInfo != NULL))
{
pBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
if (mVerbose == TRUE)
{
sprintf(msgBuf, "Block=%dx%d Type=%s,
ColorInterp=%s\n",
nBlockXSize, nBlockYSize,
GDALGetDataTypeName(pBand->GetRasterDataType()),
GDALGetColorInterpretationName(
pBand->GetColorInterpretation()) );
}
adfMinMax[0] = pBand->GetMinimum( &bGotMin );
adfMinMax[1] = pBand->GetMaximum( &bGotMax );
if( ! (bGotMin && bGotMax) )
{
GDALComputeRasterMinMax((GDALRasterBandH)pBand,
TRUE, adfMinMax);
}
Where pDataSet is the stored pointer returned by GDALOpen.
As always any thoughts or suggestions would be greatly appreciated.
Bruce Clay
bclay at ball.com
_______________________________________________
Gdal-dev mailing list
Gdal-dev at remotesensing.org
http://remotesensing.org/mailman/listinfo/gdal-dev
More information about the Gdal-dev
mailing list