[Gdal-dev] IO Problems on the Edge of Rasters

[Johnson, Donald W MVK Contractor]

 IM working on program which uses very large raster which can not be fully
loaded into memory. Instead chunks of the raster a loaded from and stored to
disk as nessessary. The entire process is hidden by an abstraction class so
that the program logic is identical to working with a large 2D array. This
appeared to work correctly but is failling on the far right hand side of
rasters. For example if the raster is 160 rows and 212 columns and is loaded
in 100x100 blocks, does not correctly right the last 12 columns to disk,
however it returns no error.

Here is a simple test case.

/////////////////////////////////////////////////////////////////////////////
///////////////////

#include "VirtualMatrixGDAL.h"
#include <fstream>
#include <string>

using namespace std;

#include "boost/test/minimal.hpp"

template <class T> void fill_matrix(T& m)
{
	for( unsigned long i = 0; i < m.rows(); ++i )
	{
		long tmp = i * m.cols();
		
		for( unsigned long j = 0; j < m.cols(); ++j )
		{
			m[i][j] = tmp + j;
		}
	}	
}

template <class Stream, class Matrix> void print_matrix(Stream& s,Matrix& m)
{
	for( unsigned long i = 0; i < m.rows(); ++i )
	{	
		for( unsigned long j = 0; j < m.cols(); ++j )
		{
			s << m[i][j] << '\t';

			if ( m[i][j] < 1000 )
			{
				s << '\t';
			}
		}

		s << '\n';
	}	
}

int test_main(int argc, char **argv)
{
	std::ofstream output("results.txt");
	
	virtual_matrix_gdal<long> test_matrix;

	test_matrix.init("test.tif",GDT_Int32,16,24,10,10);

	test_matrix.set_max_memory(5,"MB");

	fill_matrix(test_matrix);
	print_matrix(output,test_matrix);

	test_matrix.store_all_and_release_mem();
	output << '\n';

	print_matrix(output,test_matrix);

	return 0;
}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////

Calls to [][] become calls to data(unsigned long row, unsigned column);

/////////////////////////////////////////////////////////////////////////////
/////////////////////////

template < class T > T& virtual_matrix<T>::data(unsigned long r,unsigned long
c)
{
	unsigned long cr;
	unsigned long cc;
	
	// check for valid coordinates
	if ( !coordinates_valid(r,c) )
	{
		throw coordinate_range_error(r,c,num_rows,num_cols);
	}

	// get the coordinates of the correct cell
	cr = get_cell_row(r);
	cc = get_cell_col(c);
	
	virtual_matrix_cell<T>* target_cell = &cell(cr,cc);

	// load the cell if not active
	if ( target_cell->state() != cell_active )
	{
		load_cell(cr,cc);
	}



	target_cell->dirty(true);

	// get the internal row and column for the cell
	r -= (cell_rows * cr);
	c -= (cell_cols * cc);

	// remove a refrence if over max count
	if ( recent_reference.full() )
	{
		recent_reference.front()->remove_ref();
		recent_reference.pop_front();
	}

	//update the refernce counts
	
	target_cell->add_ref();
	//reference_count[cp] += 1;
	recent_reference.push_back(target_cell);

	return (*target_cell)(r,c); 
}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////

The store_all_and_release_mem simply calls the base class store_all (
store_all was overridden to keep allocated memory when storing )

/////////////////////////////////////////////////////////////////////////////
/////////////////////////

template < class T > void virtual_matrix_gdal<T>::store_all_and_release_mem()
{ super_type::store_all(); }

template < class T > void virtual_matrix<T>::store_all()
{
	for( unsigned long i = 0; i < cell.rows(); ++i )
	{
		for( unsigned long j = 0; j < cell.columns(); ++j )
		{
			if ( cell(i,j).state() == cell_active )
			{
				store_cell(i,j);
			}
			else if ( cell(i,j).state() == cell_null )
			{
				load_cell(i,j);
				store_cell(i,j);
			}
		}
	}
}

template < class T > void virtual_matrix<T>::store_cell(unsigned long r,
unsigned long c)
{
	if ( r < num_cell_rows && c < num_cell_cols )
	{
		cell(r,c).state( cell_disk);
		cell(r,c).dirty( false);

		write_cell_to_file(r,c);

		cell(r,c).resize(0,0);

		current_memory_allocation -= cell_memory_requirements(r,c);
	}
}

template < class T > void virtual_matrix<T>::load_cell(unsigned long r,
unsigned long c)
{
	virtual_matrix_cell<T>& target_cell = cell(r,c);
	
	if ( memory_avalable(r,c) )
	{		
		if ( target_cell.state() == cell_null )
		{
			target_cell.resize(cell_rows,cell_cols);
			std::fill_n( target_cell.data(), target_cell.size(),
default_val );
		}
		else if ( cell(r,c).state() == cell_disk )
		{
			target_cell.resize(cell_rows,cell_cols);
			load_cell_from_file(r,c);
		}

		current_memory_allocation += cell_memory_requirements(r,c);
	}
	else
	{
		
		// find the cell with the least recent accesses
		find_least_used_cell();

		virtual_matrix_cell<T>& swap_cell = cell(t_row,t_col);

		// write the choosen cell back to disk
		if ( swap_cell.dirty() )
		{
			write_cell_to_file(t_row,t_col);
		}

		// swap the target cell with the cell to be loaded
		target_cell.swap( swap_cell );

		// mark the swaped file as being on disk
		swap_cell.state(cell_disk);

		// load the new cell
		if ( target_cell.state() == cell_null )
		{
			std::fill_n( target_cell.data(), target_cell.size(),
default_val );
		}
		else if ( target_cell.state() == cell_disk )
		{
			load_cell_from_file(r,c);	
		}
	}

	target_cell.state( cell_active );
	target_cell.dirty(false);
}

/////////////////////////////////////////////////////////////////////////////
///////////////////////

The only function call that uses gdal in this are the calls to
write_cell_to_file and load_cell_from file.
Here is ther implementation and where I assume the error is. 

/////////////////////////////////////////////////////////////////////////////
///////////////////////

template <class T> void virtual_matrix_gdal<T>::load_cell_from_file(unsigned
long r, unsigned long c)
{
	GDALRasterBand* pBand = storageBand;
	CPLErr err = CE_None; 

	int r_offset = r * cell_rows;
	int c_offset = c * cell_cols;

	int x_size = pBand->GetXSize();
	int y_size = pBand->GetYSize();
	int read_rows = ( r_offset + cell_rows <= y_size ) ? cell_rows :
y_size - r_offset;
	int read_cols = ( c_offset + cell_cols <= x_size ) ? cell_cols :
x_size - c_offset;

	T* buffer =  cell(r,c).data();
	err =
pBand->RasterIO(GF_Read,c_offset,r_offset,read_cols,read_rows,buffer,read_col
s,read_rows,data_type,0,0);
}

template <class T> void virtual_matrix_gdal<T>::write_cell_to_file(unsigned
long r, unsigned long c)
{
	GDALRasterBand* pBand = storageBand;
	CPLErr err = CE_None;

	unsigned long r_offset = r * cell_rows;
	unsigned long c_offset = c * cell_cols;

	int x_size = pBand->GetXSize();
	int y_size = pBand->GetYSize();
	int write_rows = ( r_offset + cell_rows <= y_size ) ? cell_rows :
y_size - r_offset;
	int write_cols = ( c_offset + cell_cols <= x_size ) ? cell_cols :
x_size - c_offset;

	T* buffer =  cell(r,c).data();
	err =
pBand->RasterIO(GF_Write,c_offset,r_offset,write_cols,write_rows,buffer,write
_cols,write_rows,data_type,0,0);
}

/////////////////////////////////////////////////////////////////////////////
///////////////////////

The problem is I have step though these function when writing or reading the
cells that get corrupted and 
RasterIO() does't give any errors.

Does anyone see what im doing wrong. If someone wants the complete test
project ill gladly send the files, I didn't want to include the entire
implementation of virtual_matrix because it is a large file. This is on
Windows XP using Visual Studio 2003. The version of GDAL is 1.2.3

Thanks
 Donald J