Skip to content

Commit 51a07a6

Browse files
author
rblazek
committedMar 10, 2011
alignment better
git-svn-id: http://svn.osgeo.org/qgis/trunk@15422 c8812cc2-4d05-0410-92ff-de0c093fc19c
·
ltr-3_40Before-merge-new_vector_api
1 parent 8990dc5 commit 51a07a6

File tree

1 file changed

+71
-11
lines changed

1 file changed

+71
-11
lines changed
 

‎src/providers/gdal/qgsgdalprovider.cpp

Lines changed: 71 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -532,6 +532,11 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
532532
QgsDebugMsg( "thePixelHeight = " + QString::number( thePixelHeight ) );
533533
QgsDebugMsg( "theExtent: " + theExtent.toString() );
534534

535+
for ( int i = 0 ; i < 6; i++ )
536+
{
537+
QgsDebugMsg( QString( "transform : %1" ).arg( mGeoTransform[i] ) );
538+
}
539+
535540
// TODO: fill block with no data values
536541

537542
QgsRectangle myRasterExtent = theExtent.intersect( &mExtent );
@@ -589,46 +594,57 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
589594
double center, centerRaster;
590595
int rasterTop, rasterBottom, rasterLeft, rasterRight;
591596

597+
// We have to extend destination grid for GDALRasterIO to keep resolution:
598+
double topSpace, bottomSpace, leftSpace, rightSpace;
599+
592600
// top
593601
center = myRasterExtent.yMaximum() - yRes/2;
594602
// center in raster space
595603
// GDAL states that mGeoTransform[3] is top, but probably it can be also bottom
596-
// if mGeoTransform[5] is negative
597-
if ( mGeoTransform[5] > 0 )
604+
// if mGeoTransform[5] is negative ??? No, mGeoTransform[5] is negative normaly
605+
// - vice versa?
606+
//if ( mGeoTransform[5] > 0 )
607+
if ( mGeoTransform[5] < 0 )
598608
{
599-
centerRaster = ( mGeoTransform[3] - center ) / mGeoTransform[5];
609+
centerRaster = -1. * ( mGeoTransform[3] - center ) / mGeoTransform[5];
600610
}
601611
else
602612
{
603613
centerRaster = ( center - mGeoTransform[3] ) / mGeoTransform[5];
604614
}
605615
rasterTop = static_cast<int> ( floor ( centerRaster ) );
616+
topSpace = (mGeoTransform[3] + rasterTop*mGeoTransform[5])- myRasterExtent.yMaximum();
606617

607618
// bottom
608619
center = myRasterExtent.yMinimum() + yRes/2;
609-
if ( mGeoTransform[5] > 0 )
620+
//if ( mGeoTransform[5] > 0 )
621+
if ( mGeoTransform[5] < 0 )
610622
{
611-
centerRaster = ( mGeoTransform[3] - center ) / mGeoTransform[5];
623+
centerRaster = -1. * ( mGeoTransform[3] - center ) / mGeoTransform[5];
612624
}
613625
else
614626
{
615627
centerRaster = ( center - mGeoTransform[3] ) / mGeoTransform[5];
616628
}
617629
rasterBottom = static_cast<int> ( floor ( centerRaster ) );
630+
bottomSpace = myRasterExtent.yMinimum() - ( mGeoTransform[3] + (rasterBottom+1)*mGeoTransform[5] );
618631

619632
// left
620633
center = myRasterExtent.xMinimum() + xRes/2;
621634
centerRaster = ( center - mGeoTransform[0] ) / mGeoTransform[1];
622635
rasterLeft = static_cast<int> ( floor ( centerRaster ) );
636+
leftSpace = myRasterExtent.xMinimum() - (mGeoTransform[0] + rasterLeft * mGeoTransform[1] );
623637

624638
// right
625639
center = myRasterExtent.xMaximum() - xRes/2;
626640
centerRaster = ( center - mGeoTransform[0] ) / mGeoTransform[1];
627641
rasterRight = static_cast<int> ( floor ( centerRaster ) );
642+
rightSpace = (mGeoTransform[0] + (rasterRight+1)*mGeoTransform[1]) - myRasterExtent.xMaximum();
628643

629644
QgsDebugMsg( QString("rasterTop = %1 rasterBottom = %2 rasterLeft = %3 rasterRight = %4").arg(rasterTop).arg(rasterBottom).arg(rasterLeft).arg(rasterRight) );
630645

631-
// Allocate temporary block
646+
QgsDebugMsg( QString("topSpace = %1 bottomSpace = %2 leftSpace = %3 rightSpace = %4").arg(topSpace).arg(bottomSpace).arg(leftSpace).arg(rightSpace) );
647+
632648
int width = right - left + 1;
633649
int height = bottom - top + 1;
634650

@@ -637,6 +653,29 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
637653

638654
QgsDebugMsg( QString("width = %1 height = %2 rasterWidth = %3 rasterHeight = %4").arg(width).arg(height).arg(rasterWidth).arg(rasterHeight) );
639655

656+
657+
double rasterXRes = extent().width() / xSize();
658+
double rasterYRes = extent().height() / ySize();
659+
660+
// TODO: what is better floor/ceil, can be negative?
661+
// should be similar
662+
//double xAdd = rasterWidth*rasterXRes - width*xRes;
663+
double xAdd = leftSpace + rightSpace;
664+
int xAddPixels = static_cast<int> ( round( xAdd / xRes ) );
665+
int leftAddPixels = static_cast<int> ( round( leftSpace / xRes ) );
666+
667+
//double leftAdd = rasterWidth*rasterXRes - width*xRes;
668+
double yAdd = topSpace + bottomSpace;
669+
int yAddPixels = static_cast<int> ( round( yAdd / yRes ) );
670+
int topAddPixels = static_cast<int> ( round( topSpace / yRes ) );
671+
672+
QgsDebugMsg( QString("xAddPixels = %1 yAddPixels = %2 leftAddPixels = %3 topAddPixels = %4").arg(xAddPixels).arg(yAddPixels).arg(leftAddPixels).arg(topAddPixels) );
673+
674+
int totalWidth = width + xAddPixels;
675+
int totalHeight = height + yAddPixels;
676+
677+
QgsDebugMsg( QString("totalWidth = %1 totalHeight = %2").arg(totalWidth).arg(totalHeight) );
678+
640679
int size = dataTypeSize(theBandNo) / 8;
641680

642681
// fill with null values
@@ -650,7 +689,11 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
650689
}
651690

652691
GDALRasterBandH gdalBand = GDALGetRasterBand( mGdalDataset, theBandNo );
692+
GDALDataType type = ( GDALDataType )mGdalDataType[theBandNo-1];
693+
CPLErrorReset();
653694

695+
// This can be probably used if xAddPixels and yAddPixels are 0 to avoid memcpy
696+
#if 0
654697
// Calc beginnig of data if raster does not start at top
655698
block = (char *) theBlock;
656699
if ( top != 0 )
@@ -663,23 +706,40 @@ void QgsGdalProvider::readBlock( int theBandNo, QgsRectangle const & theExtent,
663706
if ( left != 0 ) {
664707
block += size * left;
665708
}
666-
667-
GDALDataType type = ( GDALDataType )mGdalDataType[theBandNo-1];
668-
669-
CPLErrorReset();
670709
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
671710
rasterLeft, rasterTop, rasterWidth, rasterHeight,
672711
(void *)block,
673712
width, height, type,
674713
0, nLineSpace );
714+
#endif
715+
716+
// Allocate temporary block
717+
void *tmpBlock = malloc( size * totalWidth * totalHeight );
718+
719+
CPLErrorReset();
720+
CPLErr err = GDALRasterIO( gdalBand, GF_Read,
721+
rasterLeft, rasterTop, rasterWidth, rasterHeight,
722+
(void *)tmpBlock,
723+
totalWidth, totalHeight, type,
724+
0, 0 );
675725

676726
if ( err != CPLE_None )
677727
{
678728
QgsLogger::warning( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
679729
QgsDebugMsg ( "RasterIO error: " + QString::fromUtf8( CPLGetLastErrorMsg() ) );
730+
free ( tmpBlock );
731+
return;
732+
}
733+
734+
for ( int i = 0; i < height; i++ ) {
735+
int r = i + topAddPixels;
736+
char *src = (char *)tmpBlock + size*r*totalWidth + size*leftAddPixels;
737+
char *dst = (char *)theBlock + size*(top+i)*thePixelWidth + size*(left);
738+
memcpy ( dst, src, size*width );
680739
}
681-
682740

741+
free ( tmpBlock );
742+
return;
683743
}
684744

685745
// this is old version which was using GDALWarpOperation, unfortunately

0 commit comments

Comments
 (0)
Please sign in to comment.