Skip to content

Commit

Permalink
GTiff: read GCPs in ESRI <GeodataXform> .aux.xml
Browse files Browse the repository at this point in the history 
git-svn-id: https://svn.osgeo.org/gdal/trunk@37188 f0d54148-0727-0410-94bb-9a71ac55c965
  • Loading branch information
@rouault
rouault committed Jan 19, 2017
1 parent 564092a commit 7025e35
Show file tree
Hide file tree
Showing 5 changed files with 137 additions and 20 deletions.
Binary file added autotest/gcore/data/arcgis93_geodataxform_gcp.tif
View file
Binary file not shown.
1 change: 1 addition & 0 deletions autotest/gcore/data/arcgis93_geodataxform_gcp.tif.aux.xml
View file
@@ -0,0 +1 @@
<?xml version="1.0" encoding="utf-8" ?><GeodataXform xsi:type='typens:PolynomialXform' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance' xmlns:xs='http://www.w3.org/2001/XMLSchema' xmlns:typens='http://www.esri.com/schemas/ArcGIS/9.3'><PolynomialOrder>1</PolynomialOrder><SpatialReference xsi:type='typens:ProjectedCoordinateSystem'><WKT>PROJCS[&quot;NAD_1927_UTM_Zone_12N&quot;,GEOGCS[&quot;GCS_North_American_1927&quot;,DATUM[&quot;D_North_American_1927&quot;,SPHEROID[&quot;Clarke_1866&quot;,6378206.4,294.9786982]],PRIMEM[&quot;Greenwich&quot;,0.0],UNIT[&quot;Degree&quot;,0.0174532925199433]],PROJECTION[&quot;Transverse_Mercator&quot;],PARAMETER[&quot;False_Easting&quot;,500000.0],PARAMETER[&quot;False_Northing&quot;,0.0],PARAMETER[&quot;Central_Meridian&quot;,-111.0],PARAMETER[&quot;Scale_Factor&quot;,0.9996],PARAMETER[&quot;Latitude_Of_Origin&quot;,0.0],UNIT[&quot;Meter&quot;,1.0],AUTHORITY[&quot;EPSG&quot;,26712]]</WKT><HighPrecision>true</HighPrecision><WKID>26712</WKID></SpatialReference><SourceGCPs xsi:type='typens:ArrayOfDouble'><Double>1.88333333333333</Double><Double>3.8233333333333301</Double><Double>10.516666666666699</Double><Double>3.8966666666666701</Double><Double>19.129999999999999</Double><Double>3.96</Double><Double>27.776666666666699</Double><Double>4.0300000000000002</Double><Double>1.8133333333333299</Double><Double>15.543333333333299</Double><Double>10.436666666666699</Double><Double>15.58</Double><Double>19.043333333333301</Double><Double>15.626666666666701</Double><Double>27.6666666666667</Double><Double>15.686666666666699</Double><Double>1.79666666666667</Double><Double>27.27</Double><Double>10.373333333333299</Double><Double>27.283333333333299</Double><Double>18.969999999999999</Double><Double>27.313333333333301</Double><Double>27.579999999999998</Double><Double>27.3533333333333</Double><Double>1.8200000000000001</Double><Double>38.990000000000002</Double><Double>10.383333333333301</Double><Double>38.9866666666667</Double><Double>18.940000000000001</Double><Double>38.993333333333297</Double><Double>27.523333333333301</Double><Double>39.026666666666699</Double></SourceGCPs><TargetGCPs xsi:type='typens:ArrayOfDouble'><Double>500000</Double><Double>4705078.79016612</Double><Double>513694.73894985497</Double><Double>4705092.24672045</Double><Double>527389.48823190096</Double><Double>4705132.6166430097</Double><Double>541084.25817557296</Double><Double>4705199.9007125003</Double><Double>500000</Double><Double>4723585.2528382801</Double><Double>513658.31203038001</Double><Double>4723598.71614639</Double><Double>527316.63369301497</Double><Double>4723639.1063283999</Double><Double>540974.97461739404</Double><Double>4723706.4241574099</Double><Double>500000</Double><Double>4742092.2619709503</Double><Double>513621.76888738602</Double><Double>4742105.7315773396</Double><Double>527243.54671029805</Double><Double>4742146.1406523297</Double><Double>540865.34240148403</Double><Double>4742213.4899633499</Double><Double>500000</Double><Double>4760599.8178412998</Double><Double>513585.10982275498</Double><Double>4760613.2932902602</Double><Double>527170.22788755305</Double><Double>4760653.7198910601</Double><Double>540755.362433651</Double><Double>4760721.0984054497</Double></TargetGCPs></GeodataXform>
23 changes: 23 additions & 0 deletions autotest/gcore/tiff_read.py
View file
Expand Up @@ -2826,6 +2826,28 @@ def tiff_read_unit_from_srs():

return 'success'

###############################################################################
# Test reading ArcGIS 9.3 .aux.xml

def tiff_read_arcgis93_geodataxform_gcp():

ds = gdal.Open('data/arcgis93_geodataxform_gcp.tif')
if ds.GetGCPProjection().find('26712') < 0:
gdaltest.post_reason('fail')
return 'fail'
if ds.GetGCPCount() != 16:
gdaltest.post_reason('fail')
return 'fail'
gcp = ds.GetGCPs()[0]
if abs(gcp.GCPPixel - 565) > 1e-5 or \
abs(gcp.GCPLine - 11041) > 1e-5 or \
abs(gcp.GCPX - 500000) > 1e-5 or \
abs(gcp.GCPY - 4705078.79016612) > 1e-5 or \
abs(gcp.GCPZ - 0) > 1e-5:
gdaltest.post_reason('fail')
return 'fail'
return 'success'

###############################################################################

for item in init_list:
Expand Down Expand Up @@ -2918,6 +2940,7 @@ def tiff_read_unit_from_srs():
gdaltest_list.append( (tiff_read_ycbcr_lzw) )

gdaltest_list.append( (tiff_read_unit_from_srs) )
gdaltest_list.append( (tiff_read_arcgis93_geodataxform_gcp) )

# gdaltest_list = [ tiff_read_ycbcr_lzw ]

Expand Down
74 changes: 74 additions & 0 deletions gdal/frmts/gtiff/geotiff.cpp
View file
Expand Up @@ -11877,6 +11877,80 @@ void GTiffDataset::ApplyPamInfo()
bLookedForProjection = true;
}

if( m_nPAMGeorefSrcIndex >= 0 && nGCPCount == 0 )
{
CPLXMLNode *psValueAsXML = NULL;
CPLXMLNode *psGeodataXform = NULL;
char** papszXML = oMDMD.GetMetadata( "xml:ESRI" );
if (CSLCount(papszXML) == 1)
{
psValueAsXML = CPLParseXMLString( papszXML[0] );
if( psValueAsXML )
psGeodataXform = CPLGetXMLNode(psValueAsXML, "=GeodataXform");
}

const char* pszTIFFTagResUnit = GetMetadataItem("TIFFTAG_RESOLUTIONUNIT");
const char* pszTIFFTagXRes = GetMetadataItem("TIFFTAG_XRESOLUTION");
const char* pszTIFFTagYRes = GetMetadataItem("TIFFTAG_YRESOLUTION");
if (psGeodataXform && pszTIFFTagResUnit &&pszTIFFTagXRes &&
pszTIFFTagYRes && atoi(pszTIFFTagResUnit) == 2 )
{
CPLXMLNode* psSourceGCPs = CPLGetXMLNode(psGeodataXform,
"SourceGCPs");
CPLXMLNode* psTargetGCPs = CPLGetXMLNode(psGeodataXform,
"TargetGCPs");
if( psSourceGCPs && psTargetGCPs )
{
std::vector<double> adfSourceGCPs, adfTargetGCPs;
for( CPLXMLNode* psIter = psSourceGCPs->psChild;
psIter != NULL;
psIter = psIter->psNext )
{
if( psIter->eType == CXT_Element &&
EQUAL(psIter->pszValue, "Double") )
{
adfSourceGCPs.push_back(
CPLAtof( CPLGetXMLValue(psIter, NULL, "") ) );
}
}
for( CPLXMLNode* psIter = psTargetGCPs->psChild;
psIter != NULL;
psIter = psIter->psNext )
{
if( psIter->eType == CXT_Element &&
EQUAL(psIter->pszValue, "Double") )
{
adfTargetGCPs.push_back(
CPLAtof( CPLGetXMLValue(psIter, NULL, "") ) );
}
}
if( adfSourceGCPs.size() == adfTargetGCPs.size() &&
(adfSourceGCPs.size() % 2) == 0 )
{
nGCPCount = static_cast<int>(
adfSourceGCPs.size() / 2);
pasGCPList = static_cast<GDAL_GCP *>(
CPLCalloc(sizeof(GDAL_GCP), nGCPCount) );
for( int i = 0; i < nGCPCount; ++i )
{
pasGCPList[i].pszId = CPLStrdup("");
pasGCPList[i].pszInfo = CPLStrdup("");
// The origin used is the bottom left corner,
// and raw values are in inches !
pasGCPList[i].dfGCPPixel = adfSourceGCPs[2*i] *
CPLAtof(pszTIFFTagXRes);
pasGCPList[i].dfGCPLine = nRasterYSize -
adfSourceGCPs[2*i+1] * CPLAtof(pszTIFFTagYRes);
pasGCPList[i].dfGCPX = adfTargetGCPs[2*i];
pasGCPList[i].dfGCPY = adfTargetGCPs[2*i+1];
}
}
}
}
if( psValueAsXML )
CPLDestroyXMLNode(psValueAsXML);
}

/* -------------------------------------------------------------------- */
/* Copy any PAM metadata into our GeoTIFF context, and with */
/* the PAM info overriding the GeoTIFF context. */
Expand Down
59 changes: 39 additions & 20 deletions gdal/gcore/gdalpamdataset.cpp
View file
Expand Up @@ -414,37 +414,56 @@ CPLErr GDALPamDataset::XMLInit( CPLXMLNode *psTree, const char *pszUnused )
oMDMD.XMLInit( psTree, TRUE );

/* -------------------------------------------------------------------- */
/* Try loading ESRI xml encoded projection */
/* Try loading ESRI xml encoded GeodataXform. */
/* -------------------------------------------------------------------- */
if (psPam->pszProjection == NULL)
{
char** papszXML = oMDMD.GetMetadata( "xml:ESRI" );
if (CSLCount(papszXML) == 1)
// ArcGIS 9.3: GeodataXform as a root element
CPLXMLNode* psGeodataXform = CPLGetXMLNode(psTree, "=GeodataXform");
CPLXMLNode *psValueAsXML = NULL;
if( psGeodataXform != NULL )
{
CPLXMLNode *psValueAsXML = CPLParseXMLString( papszXML[0] );
if (psValueAsXML)
char* apszMD[2];
apszMD[0] = CPLSerializeXMLTree(psGeodataXform);
apszMD[1] = NULL;
oMDMD.SetMetadata( apszMD, "xml:ESRI");
CPLFree(apszMD[0]);
}
else
{
// ArcGIS 10: GeodataXform as content of xml:ESRI metadata domain.
char** papszXML = oMDMD.GetMetadata( "xml:ESRI" );
if (CSLCount(papszXML) == 1)
{
psValueAsXML = CPLParseXMLString( papszXML[0] );
if( psValueAsXML )
psGeodataXform = CPLGetXMLNode(psValueAsXML, "=GeodataXform");
}
}

if (psGeodataXform)
{
const char* pszESRI_WKT = CPLGetXMLValue(psGeodataXform,
"SpatialReference.WKT", NULL);
if (pszESRI_WKT)
{
const char* pszESRI_WKT = CPLGetXMLValue(psValueAsXML,
"=GeodataXform.SpatialReference.WKT", NULL);
if (pszESRI_WKT)
OGRSpatialReference* poSRS = new OGRSpatialReference(NULL);
char* pszTmp = (char*)pszESRI_WKT;
if (poSRS->importFromWkt(&pszTmp) == OGRERR_NONE &&
poSRS->morphFromESRI() == OGRERR_NONE)
{
OGRSpatialReference* poSRS = new OGRSpatialReference(NULL);
char* pszTmp = (char*)pszESRI_WKT;
if (poSRS->importFromWkt(&pszTmp) == OGRERR_NONE &&
poSRS->morphFromESRI() == OGRERR_NONE)
char* pszWKT = NULL;
if (poSRS->exportToWkt(&pszWKT) == OGRERR_NONE)
{
char* pszWKT = NULL;
if (poSRS->exportToWkt(&pszWKT) == OGRERR_NONE)
{
psPam->pszProjection = CPLStrdup(pszWKT);
}
CPLFree(pszWKT);
psPam->pszProjection = CPLStrdup(pszWKT);
}
delete poSRS;
CPLFree(pszWKT);
}
CPLDestroyXMLNode(psValueAsXML);
delete poSRS;
}
}
if( psValueAsXML )
CPLDestroyXMLNode(psValueAsXML);
}

/* -------------------------------------------------------------------- */
Expand Down

0 comments on commit 7025e35

Please sign in to comment.