From fa7af0b95804fac670a88f602f9cbec922b4d054 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 30 Sep 2016 11:12:06 +0000 Subject: [PATCH] HDF5: correct number of GCPs to avoid dummy trailing (0,0)->(0,0,0) and remove +180 offset applied to GCP longitude. Add instead a heuristics to determine if the product is crossing the antimeridian, and a HDF5_SHIFT_GCPX_BY_180 config option to be able to override the heuristics (fixes #6666) git-svn-id: https://svn.osgeo.org/gdal/branches/2.1@35557 f0d54148-0727-0410-94bb-9a71ac55c965 --- gdal/frmts/hdf5/hdf5imagedataset.cpp | 33 +++++++++++++++++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/gdal/frmts/hdf5/hdf5imagedataset.cpp b/gdal/frmts/hdf5/hdf5imagedataset.cpp index 120a04a..e6211c1 100644 --- a/gdal/frmts/hdf5/hdf5imagedataset.cpp +++ b/gdal/frmts/hdf5/hdf5imagedataset.cpp @@ -825,7 +825,7 @@ CPLErr HDF5ImageDataset::CreateProjections() /* -------------------------------------------------------------------- */ /* Fill the GCPs list. */ /* -------------------------------------------------------------------- */ - nGCPCount = nRasterYSize/nDeltaLat * nRasterXSize/nDeltaLon; + nGCPCount = (nRasterYSize/nDeltaLat) * (nRasterXSize/nDeltaLon); pasGCPList = static_cast( CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) ) ); @@ -835,12 +835,41 @@ CPLErr HDF5ImageDataset::CreateProjections() const int nYLimit = (static_cast(nRasterYSize)/nDeltaLat) * nDeltaLat; const int nXLimit = (static_cast(nRasterXSize)/nDeltaLon) * nDeltaLon; + + // The original code in https://trac.osgeo.org/gdal/changeset/8066 + // always add +180 to the longitudes, but without justification + // I suspect this might be due to handling products crossing the + // antimeridian. Trying to do it just when needed through a heuristics. + bool bHasLonNearMinus180 = false; + bool bHasLonNearPlus180 = false; + bool bHasLonNearZero = false; + for( int j = 0; j < nYLimit; j+=nDeltaLat ) + { + for( int i = 0; i < nXLimit; i+=nDeltaLon ) + { + const int iGCP = j * nRasterXSize + i; + if( Longitude[iGCP] > 170 && Longitude[iGCP] <= 180 ) + bHasLonNearPlus180 = true; + if( Longitude[iGCP] < -170 && Longitude[iGCP] >= -180 ) + bHasLonNearMinus180 = true; + if( fabs(Longitude[iGCP]) < 90 ) + bHasLonNearZero = true; + } + } + const char* pszShiftGCP = + CPLGetConfigOption("HDF5_SHIFT_GCPX_BY_180", NULL); + const bool bAdd180 = (bHasLonNearPlus180 && bHasLonNearMinus180 && + !bHasLonNearZero && pszShiftGCP == NULL) || + (pszShiftGCP != NULL && CPLTestBool(pszShiftGCP)); + for( int j = 0; j < nYLimit; j+=nDeltaLat ) { for( int i = 0; i < nXLimit; i+=nDeltaLon ) { const int iGCP = j * nRasterXSize + i; - pasGCPList[k].dfGCPX = static_cast(Longitude[iGCP])+180.0; + pasGCPList[k].dfGCPX = static_cast(Longitude[iGCP]); + if( bAdd180 ) + pasGCPList[k].dfGCPX += 180.0; pasGCPList[k].dfGCPY = static_cast(Latitude[iGCP]); pasGCPList[k].dfGCPPixel = i + 0.5;