Commit c874430e authored by Julien Osman's avatar Julien Osman
Browse files

BUG: Don't rely on GDAL to read TerrasarX GCP

parent 684814b1
Pipeline #7550 failed with stages
in 99 minutes and 54 seconds
SensorID SAR
Mission TSX-1
ProductType SSC____SM_S
Polarization HH
Mode SM
OrbitDirection DESCENDING
TileHintX 15328
TileHintY 1
DataType 11
OrbitNumber 4091
CalScale 1.04345e-05
CalFactor 1.04345e-05
PRF 3832.52
RSF 1.09886e+08
......@@ -17,10 +19,9 @@ RangeTimeFirstPixel 0.00422717
RangeTimeLastPixel 0.00436665
AcquisitionStartTime 2008-03-10T13:32:20.083162Z
AcquisitionStopTime 2008-03-10T13:32:28.617747Z
SAR <SARParam>
SARCalib <SARCalib>
GCP <GCPParam>
Polarization HH
SAR <SARParam>
{"Projection": "GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]",
[{"GCP_Id": "1", "GCP_Info": "", "GCP_Row": "1", "GCP_Col": "1", "GCP_X": "-111.665", "GCP_Y": "36.2693", "GCP_Z": "0", },
{"GCP_Id": "2", "GCP_Info": "", "GCP_Row": "1", "GCP_Col": "232", "GCP_X": "-111.671", "GCP_Y": "36.2701", "GCP_Z": "0", },
......
This diff is collapsed.
......@@ -28,7 +28,7 @@
#include "otbImageKeywordlist.h"
#include "otbXMLMetadataSupplier.h"
#include "otbSARMetadata.h"
#include "otbSpatialReference.h"
#include <boost/filesystem.hpp>
//TODO C++ 17 : use std::optional instead
......@@ -953,6 +953,54 @@ void TerraSarImageMetadataInterface::PrintSelf(std::ostream& os, itk::Indent ind
Superclass::PrintSelf(os, indent);
}
void ReadGeorefGCP(const XMLMetadataSupplier & geoXmlMS, ImageMetadata& imd)
{
Projection::GCPParam gcp;
std::stringstream ss;
// Get the ellipsoid and semi-major, semi-minor axes
if(geoXmlMS.HasValue("geoReference.referenceFrames.sphere.ellipsoidID"))
{
auto ellipsoidID = geoXmlMS.GetAs<std::string>("geoReference.referenceFrames.sphere.ellipsoidID");
auto minor_axis = geoXmlMS.GetAs<double>(0, "geoReference.referenceFrames.sphere.semiMinorAxis");
auto major_axis = geoXmlMS.GetAs<double>(0, "geoReference.referenceFrames.sphere.semiMajorAxis");
if(ellipsoidID.empty() || minor_axis == 0 || major_axis == 0)
{
otbGenericExceptionMacro(MissingMetadataException, << "Cannot read GCP's spatial reference");
}
else if(ellipsoidID == "WGS84")
{
gcp.GCPProjection = SpatialReference::FromWGS84().ToWkt();
}
else
{
gcp.GCPProjection = SpatialReference::FromGeogCS("", "", ellipsoidID, major_axis,
major_axis/(major_axis - minor_axis)).ToWkt();
}
}
auto GCPCount = geoXmlMS.GetAs<unsigned int>(0, "geoReference.geolocationGrid.numberOfGridPoints.total");
// Count the gcps if the given count value is invalid
if(GCPCount == 0)
GCPCount = geoXmlMS.GetNumberOf("geoReference.geolocationGrid.gridPoint");
// Put some reasonable limits of the number of gcps
if(GCPCount > 5000)
GCPCount = 5000;
for(unsigned int i = 1 ; i <= GCPCount ; ++i)
{
ss.str("");
ss << "geoReference.geolocationGrid.gridPoint_" << i << ".";
gcp.GCPs.emplace_back(std::to_string(i), // id
"", // info
geoXmlMS.GetAs<double>(ss.str() + "col"), // col
geoXmlMS.GetAs<double>(ss.str() + "row"), // row
geoXmlMS.GetAs<double>(ss.str() + "lon"), // px
geoXmlMS.GetAs<double>(ss.str() + "lat"), // py
0); // pz ("height" in the xml file, but GDAL doesn't read it, so we do the same)
}
imd.Add(MDGeom::GCP, gcp);
}
void ReadSARSensorModel(const XMLMetadataSupplier & xmlMS,
const XMLMetadataSupplier & geoXmlMS,
......@@ -1123,20 +1171,15 @@ void TerraSarImageMetadataInterface::ParseGdal(ImageMetadata &imd)
LoadRadiometricCalibrationData(sarCalib, MainXMLFileMetadataSupplier, imd, polarization);
imd.Add(MDGeom::SARCalib, sarCalib);
assert(m_MetadataSupplierInterface->GetNbBands() == imd.Bands.size());
SARParam sarParam;
// Open the georef file containing GCPs
XMLMetadataSupplier GCPXMLFileMS(MainDirectory + "/ANNOTATION/GEOREF.xml");
if(imd.Has(MDGeom::GCP))
ReadSARSensorModel(MainXMLFileMetadataSupplier, GCPXMLFileMS, imd.GetGCPParam(), polarization, sarParam);
else
{
Projection::GCPParam gcp;
ReadSARSensorModel(MainXMLFileMetadataSupplier, GCPXMLFileMS, gcp, polarization, sarParam);
}
// Fetch the GCP
if(!imd.Has(MDGeom::GCP))
ReadGeorefGCP(GCPXMLFileMS, imd);
SARParam sarParam;
ReadSARSensorModel(MainXMLFileMetadataSupplier, GCPXMLFileMS, imd.GetGCPParam(), polarization, sarParam);
imd.Add(MDGeom::SAR, sarParam);
}
......@@ -1155,6 +1198,7 @@ void TerraSarImageMetadataInterface::ParseGeom(ImageMetadata & imd)
Fetch(MDNum::PRF, imd, "sensor_params.prf");
Fetch(MDNum::CalFactor, imd, "calibration.calibrationConstant.calFactor");
Fetch(MDNum::CalScale, imd, "calibration.calibrationConstant.calFactor");
Fetch(MDStr::Polarization, imd, "acquisitionInfo.polarisationList[0]");
// Main XML file
std::string MainFilePath = ExtractXMLFiles::GetResourceFile(itksys::SystemTools::GetFilenamePath(m_MetadataSupplierInterface->GetResourceFile("")),
......@@ -1166,13 +1210,18 @@ void TerraSarImageMetadataInterface::ParseGeom(ImageMetadata & imd)
imd.Add(MDStr::OrbitDirection, MainXMLFileMS.GetAs<std::string>("level1Product.productInfo.missionInfo.orbitDirection"));
imd.Add(MDNum::OrbitNumber, MainXMLFileMS.GetAs<double>("level1Product.productInfo.missionInfo.absOrbit"));
imd.Add(MDNum::RSF, MainXMLFileMS.GetAs<double>("level1Product.productSpecific.complexImageInfo.commonRSF"));
}
assert(m_MetadataSupplierInterface->GetNbBands() == imd.Bands.size());
// Open the georef file containing GCPs
XMLMetadataSupplier GCPXMLFileMS(itksys::SystemTools::GetParentDirectory(MainFilePath) + "/ANNOTATION/GEOREF.xml");
SARParam sarParam;
Fetch(MDStr::Polarization, imd, "acquisitionInfo.polarisationList[0]");
imd.Add(MDGeom::SAR, sarParam);
// Fetch the GCP
if(!imd.Has(MDGeom::GCP))
ReadGeorefGCP(GCPXMLFileMS, imd);
SARParam sarParam;
ReadSARSensorModel(MainXMLFileMS, GCPXMLFileMS, imd.GetGCPParam(), imd[MDStr::Polarization], sarParam);
imd.Add(MDGeom::SAR, sarParam);
}
SARCalib sarCalib;
LoadRadiometricCalibrationData(sarCalib, *m_MetadataSupplierInterface, imd);
......
......@@ -453,7 +453,7 @@ set(ikonos_filename LARGEINPUT{IKONOS/BLOSSEVILLE/po_2619900_blu_0000000.tif})
set(worldview2_filename LARGEINPUT{WORLDVIEW2/ROME/WV-2_standard_8band_bundle_16bit/052298844010_01_P001_MUL/09DEC10103019-M2AS-052298844010_01_P001.TIF})
set(quickbird_filename LARGEINPUT{QUICKBIRD/CEVENNES/06FEB12104912-P1BS-005533998070_01_P001.TIF})
set(Radarsat2_filename LARGEINPUT{RADARSAT2/ALTONA/Fine_Quad-Pol_Dataset/PK6621_DK406_FQ9_20080405_124900_HH_VV_HV_VH_SLC_Altona})
set(TerraSar_filename LARGEINPUT{TERRASARX/2008-03-10_GrandCanyon_SSC/TSX1_SAR__SSC______SM_S_SRA_20080310T133220_20080310T133228})
set(TerraSar_filename LARGEINPUT{TERRASARX/2008-03-10_GrandCanyon_SSC/TSX1_SAR__SSC______SM_S_SRA_20080310T133220_20080310T133228/IMAGEDATA/IMAGE_HH_SRA_strip_011.cos})
set(CosmoSkyMed_filename LARGEINPUT{COSMOSKYMED/Toulouse_spotlight/CSKS3_SCS_B_S2_08_HH_RD_SF_20110418180325_20110418180332.h5})
set(Sentinel1_filename LARGEINPUT{SENTINEL1/S1A_S6_SLC__1SSV_20150619T195043})
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment