Skip to content
Permalink
Browse files

ENVI: add read support for reading GCPs (OSGeo#1528), and fix off-by-…

…one offset on line,pixel on reading GCP
  • Loading branch information...
rouault committed May 10, 2019
1 parent 3709876 commit f2f29cd3a7708a4a9553f127b2d19b3cd72f9180
Showing with 106 additions and 3 deletions.
  1. +26 −0 autotest/gdrivers/envi.py
  2. +75 −3 gdal/frmts/raw/envidataset.cpp
  3. +5 −0 gdal/frmts/raw/envidataset.h
@@ -351,6 +351,32 @@ def test_envi_truncated():

assert cs == 2315

###############################################################################
# Test writing & reading GCPs (#1528)


def test_envi_gcp():

filename = '/vsimem/test_envi_gcp.dat'
ds = gdal.GetDriverByName('ENVI').Create(filename, 1, 1)
gcp = gdal.GCP()
gcp.GCPPixel = 1
gcp.GCPLine = 2
gcp.GCPX = 3
gcp.GCPY = 4
ds.SetGCPs([gcp], None)
ds = None
gdal.Unlink(filename + ".aux.xml")

ds = gdal.Open(filename)
assert ds.GetGCPCount() == 1
gcps = ds.GetGCPs()
assert len(gcps) == 1
gcp = gcps[0]
ds = None
assert gcp.GCPPixel == 1
assert gcp.GCPLine == 2
assert gcp.GCPX == 3
assert gcp.GCPY == 4

gdal.GetDriverByName('ENVI').Delete(filename)
@@ -305,6 +305,10 @@ ENVIDataset::~ENVIDataset()
CPLError(CE_Failure, CPLE_FileIO, "I/O error");
}
}
if( !m_asGCPs.empty() )
{
GDALDeinitGCPs(static_cast<int>(m_asGCPs.size()), m_asGCPs.data());
}
CPLFree(pszProjection);
CPLFree(pszHDRFilename);
}
@@ -1106,7 +1110,7 @@ bool ENVIDataset::WritePseudoGcpInfo()
// Write out gcps into the envi header
// returns 0 if the gcps are not present.

const int iNum = GetGCPCount();
const int iNum = std::min(GetGCPCount(), 4);
if (iNum == 0)
return false;

@@ -1120,9 +1124,11 @@ bool ENVIDataset::WritePseudoGcpInfo()
bool bRet = VSIFPrintfL(fp, "geo points = {\n") >= 0;
for( int iR = 0; iR < iNum; iR++ )
{
// Add 1 to pixel and line for ENVI convention
bRet &= VSIFPrintfL(
fp, " %#0.4f, %#0.4f, %#0.8f, %#0.8f",
pGcpStructs[iR].dfGCPPixel, pGcpStructs[iR].dfGCPLine,
1 + pGcpStructs[iR].dfGCPPixel,
1 + pGcpStructs[iR].dfGCPLine,
pGcpStructs[iR].dfGCPY, pGcpStructs[iR].dfGCPX) >= 0;
if( iR < iNum - 1 )
bRet &= VSIFPrintfL(fp, ",\n") >= 0;
@@ -1756,6 +1762,65 @@ void ENVIDataset::ProcessRPCinfo( const char *pszRPCinfo,
CSLDestroy(papszFields);
}

/************************************************************************/
/* GetGCPCount() */
/************************************************************************/

int ENVIDataset::GetGCPCount()
{
int nGCPCount = RawDataset::GetGCPCount();
if( nGCPCount )
return nGCPCount;
return static_cast<int>(m_asGCPs.size());
}

/************************************************************************/
/* GetGCPs() */
/************************************************************************/

const GDAL_GCP *ENVIDataset::GetGCPs()
{
int nGCPCount = RawDataset::GetGCPCount();
if( nGCPCount )
return RawDataset::GetGCPs();
if( !m_asGCPs.empty() )
return m_asGCPs.data();
return nullptr;
}

/************************************************************************/
/* ProcessGeoPoints() */
/* */
/* Extract GCPs */
/************************************************************************/

void ENVIDataset::ProcessGeoPoints( const char *pszGeoPoints )
{
char **papszFields = SplitList(pszGeoPoints);
const int nCount = CSLCount(papszFields);

if( (nCount % 4) != 0 )
{
CSLDestroy(papszFields);
return;
}
m_asGCPs.resize(nCount / 4);
if( !m_asGCPs.empty() )
{
GDALInitGCPs(static_cast<int>(m_asGCPs.size()), m_asGCPs.data());
}
for( int i = 0; i < static_cast<int>(m_asGCPs.size()); i++ )
{
// Substract 1 to pixel and line for ENVI convention
m_asGCPs[i].dfGCPPixel = CPLAtof( papszFields[i * 4 + 0] ) - 1;
m_asGCPs[i].dfGCPLine = CPLAtof( papszFields[i * 4 + 1] ) - 1;
m_asGCPs[i].dfGCPY = CPLAtof( papszFields[i * 4 + 2] );
m_asGCPs[i].dfGCPX = CPLAtof( papszFields[i * 4 + 3] );
m_asGCPs[i].dfGCPZ = 0;
}
CSLDestroy(papszFields);
}

void ENVIDataset::ProcessStatsFile()
{
osStaFilename = CPLResetExtension(pszHDRFilename, "sta");
@@ -2505,14 +2570,21 @@ GDALDataset *ENVIDataset::Open( GDALOpenInfo *poOpenInfo, bool bFileSizeCheck )
pszMapInfo));
}

// Look for RPC mapinfo.
// Look for RPC.
const char* pszRPCInfo = poDS->m_aosHeader["rpc_info"];
if( !poDS->bFoundMapinfo && pszRPCInfo != nullptr )
{
poDS->ProcessRPCinfo(pszRPCInfo,
poDS->nRasterXSize, poDS->nRasterYSize);
}

// Look for geo_points / GCP
const char* pszGeoPoints = poDS->m_aosHeader["geo_points"];
if( !poDS->bFoundMapinfo && pszGeoPoints != nullptr )
{
poDS->ProcessGeoPoints(pszGeoPoints);
}

// Initialize any PAM information.
poDS->SetDescription(poOpenInfo->pszFilename);
poDS->TryLoadXML();
@@ -80,9 +80,12 @@ class ENVIDataset final: public RawDataset

CPLString osStaFilename{};

std::vector<GDAL_GCP> m_asGCPs{};

bool ReadHeader( VSILFILE * );
bool ProcessMapinfo( const char * );
void ProcessRPCinfo( const char *, int, int);
void ProcessGeoPoints( const char* );
void ProcessStatsFile();
static int byteSwapInt(int);
static float byteSwapFloat(float);
@@ -130,6 +133,8 @@ class ENVIDataset final: public RawDataset
const char *pszDomain = "" ) override;
CPLErr SetGCPs( int nGCPCount, const GDAL_GCP *pasGCPList,
const OGRSpatialReference* poSRS ) override;
int GetGCPCount() override;
const GDAL_GCP *GetGCPs() override;

static GDALDataset *Open( GDALOpenInfo * );
static GDALDataset *Open( GDALOpenInfo *, bool bFileSizeCheck );

0 comments on commit f2f29cd

Please sign in to comment.
You can’t perform that action at this time.