-
Christof Kaufmann authoredChristof Kaufmann authored
geoinfo.cpp 33.01 KiB
#include <cstring>
#include <cmath>
#include <cassert>
#include "filesystem.h"
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <gdal_utils.h>
#include <cpl_string.h>
#include "geoinfo.h"
#include "exceptions.h"
namespace {
std::vector<imagefusion::Coordinate> convertProj(imagefusion::GeoInfo const& giSrc, imagefusion::GeoInfo const& giDst, std::vector<imagefusion::Coordinate> const& coords, bool srcImg, bool dstImg) {
// collect coordinates
int nSamplePoints = coords.size();
std::vector<double> x;
std::vector<double> y;
std::vector<double> z(nSamplePoints, 0.0);
std::vector<int> succ(nSamplePoints);
x.reserve(nSamplePoints);
y.reserve(nSamplePoints);
for (imagefusion::Coordinate const& c : coords) {
x.push_back(c.x);
y.push_back(c.y);
}
// build the transformer object
char* pszSrcWKT = nullptr;
giSrc.geotransSRS.exportToWkt(&pszSrcWKT);
if (std::strlen(pszSrcWKT) <= 0) {
IF_THROW_EXCEPTION(imagefusion::invalid_argument_error(
"Please provide a valid source projection reference system for reprojection."));
}
char* pszDstWKT = nullptr;
giDst.geotransSRS.exportToWkt(&pszDstWKT);
if (std::strlen(pszDstWKT) <= 0) {
if (pszSrcWKT)
CPLFree(pszSrcWKT);
IF_THROW_EXCEPTION(imagefusion::invalid_argument_error(
"Please provide a valid destination projection reference system for reprojection."));
}
double const* padfSrcGeoTransform = srcImg ? giSrc.geotrans.values.data() : nullptr;
double const* padfDstGeoTransform = dstImg ? giDst.geotrans.values.data() : nullptr;
void* pTransformer = GDALCreateGenImgProjTransformer3(pszSrcWKT, padfSrcGeoTransform,
pszDstWKT, padfDstGeoTransform);
// transform in place
constexpr int fromDstToSrc = FALSE;
GDALGenImgProjTransform(pTransformer, fromDstToSrc, nSamplePoints,
x.data(), y.data(), z.data(), succ.data());
CPLFree(pszSrcWKT);
CPLFree(pszDstWKT);
GDALDestroyGenImgProjTransformer(pTransformer);
// convert to coordinates again
std::vector<imagefusion::Coordinate> ret;
ret.reserve(nSamplePoints);
for (int i = 0; i < nSamplePoints; ++i)
ret.emplace_back(x.at(i), y.at(i));
return ret;
}
std::vector<imagefusion::Coordinate> convertLongLat(imagefusion::GeoInfo const& gi, std::vector<imagefusion::Coordinate> const& coords, bool imgCoords /* otherwise proj coords */, bool reverse /* from long lat? */) {
// collect coordinates
int nSamplePoints = coords.size();
std::vector<double> x;
std::vector<double> y;
x.reserve(nSamplePoints);
y.reserve(nSamplePoints);
for (imagefusion::Coordinate c : coords) {
if (imgCoords && !reverse)
c = gi.geotrans.imgToProj(c);
x.push_back(c.x);
y.push_back(c.y);
}
// prepare long/lat target CRS, see https://gdal.org/tutorials/osr_api_tut.html#coordinate-transformation
OGRSpatialReference* poLongLat = gi.geotransSRS.CloneGeogCS();
#if GDAL_VERSION_MAJOR >= 3
poLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
#endif
// build the transformer object
OGRCoordinateTransformation* poTransform;
if (reverse) {
poTransform = OGRCreateCoordinateTransformation(poLongLat,
const_cast<OGRSpatialReference*>(&gi.geotransSRS));
}
else {
poTransform = OGRCreateCoordinateTransformation(const_cast<OGRSpatialReference*>(&gi.geotransSRS),
poLongLat);
}
if (poTransform == nullptr)
IF_THROW_EXCEPTION(imagefusion::invalid_argument_error(
"Cannot make the transformation to or from latitude / longitude. Please provide a valid projection reference system."));
// transform in place
if (!poTransform->Transform(nSamplePoints, x.data(), y.data())) {
OGRCoordinateTransformation::DestroyCT(poTransform);
IF_THROW_EXCEPTION(imagefusion::invalid_argument_error(
"Could not transform to or from latitude / longitude. Don't know why, sorry! Maybe the points are given in the wrong coordinate space and thus are not in the expected range."));
}
OGRCoordinateTransformation::DestroyCT(poTransform);
OGRSpatialReference::DestroySpatialReference(poLongLat);
// convert to coordinates again
std::vector<imagefusion::Coordinate> ret;
ret.reserve(nSamplePoints);
if (imgCoords && reverse)
for (int i = 0; i < nSamplePoints; ++i)
ret.emplace_back(gi.geotrans.projToImg({x.at(i), y.at(i)}));
else
for (int i = 0; i < nSamplePoints; ++i)
ret.emplace_back(x.at(i), y.at(i));
return ret;
}
} /* anonymous namespace */
namespace imagefusion {
void GeoInfo::readFrom(std::string const& filename) {
GDALDataset* img = (GDALDataset*) GDALOpen(filename.c_str(), GA_ReadOnly);
if (!img)
IF_THROW_EXCEPTION(runtime_error("Could not open image '" + filename + "' with GDAL to read its GeoInfo. "
"Either the file does not exist or GDAL could not find an appropriate driver to read the image"))
<< boost::errinfo_file_name(filename);
readFrom(img);
GDALClose(img);
}
void GeoInfo::readFrom(GDALDataset const* img_arg) {
if (!img_arg)
IF_THROW_EXCEPTION(runtime_error("The GDALDataset from which the GeoInfo should be read is nullptr."));
GDALDataset* img = const_cast<GDALDataset*>(img_arg);
// filename
char** fileList = img->GetFileList();
if (CSLCount(fileList) > 0)
filename = CSLGetField(fileList, 0);
if (fileList)
CSLDestroy(fileList);
// determine size, channels, baseType, colorTable
size.width = img->GetRasterXSize();
size.height = img->GetRasterYSize();
channels = img->GetRasterCount();
GDALColorTable* ct = nullptr;
if (channels > 0) {
GDALDataType gdal_type = img->GetRasterBand(1)->GetRasterDataType();
baseType = toBaseType(gdal_type);
for (int i = 2; i <= channels; ++i) {
if (img->GetRasterBand(i)->GetRasterDataType() != gdal_type) {
baseType = Type::invalid;
break;
}
}
ct = img->GetRasterBand(1)->GetColorTable();
if (ct) {
if (channels != 1)
IF_THROW_EXCEPTION(file_format_error("Color indexed images are only supported for single-index images. This image has " + std::to_string(channels) + " indices."))
<< boost::errinfo_file_name(filename);
if (baseType != Type::uint8)
IF_THROW_EXCEPTION(file_format_error("This is a color indexed image, which is not of type uint8 (or GDT_Byte). I did not know this is possible. Please help me to fix this!"))
<< boost::errinfo_file_name(filename);
int ctsize = ct->GetColorEntryCount();
colorTable.clear();
colorTable.reserve(ctsize);
for (int i = 0; i < ctsize; ++i) {
GDALColorEntry const* ce = ct->GetColorEntry(i);
colorTable.push_back(std::array<short,4>{ce->c1, ce->c2, ce->c3, ce->c4});
}
}
}
else
baseType = Type::invalid;
// no data values
clearNodataValues();
bool hasAnyNodataVals = false;
for (int i = 1; i <= channels; ++i) {
int hasNodataVal;
double nodataVal = img->GetRasterBand(i)->GetNoDataValue(&hasNodataVal);
if (hasNodataVal) {
hasAnyNodataVals = true;
setNodataValue(nodataVal, i-1);
}
}
if (hasAnyNodataVals)
nodataValues.resize(channels, std::numeric_limits<double>::quiet_NaN());
// geo transform and projection
CPLErr errorGeotransform = img->GetGeoTransform(geotrans.values.data());
const char* geotransProjStr = img->GetProjectionRef();
if (errorGeotransform == CE_None && std::strlen(geotransProjStr) > 0) {
geotransSRS.SetFromUserInput(geotransProjStr);
// geotransSRS.importFromWkt(const_cast<char**>(&geotransProjStr)); // same but const cast required!?
}
// ground control points and projection
int gcp_count = img->GetGCPCount();
gcps.clear();
if (gcp_count > 0) {
GDAL_GCP const* gdal_gcps = img->GetGCPs();
for (int i = 0; i < gcp_count; ++i) {
gcps.push_back(GCP{gdal_gcps->pszId,
gdal_gcps->pszInfo,
gdal_gcps->dfGCPPixel,
gdal_gcps->dfGCPLine,
gdal_gcps->dfGCPX,
gdal_gcps->dfGCPY,
gdal_gcps->dfGCPZ});
gdal_gcps++;
}
const char* gcpProjStr = img->GetGCPProjection();
if (std::strlen(gcpProjStr) > 0)
gcpSRS.SetFromUserInput(gcpProjStr);
}
// remaining meta data
char** domains = img->GetMetadataDomainList();
if (domains != nullptr) {
for (char** dom = domains; *dom != nullptr; ++dom) {
char** metaitems = img->GetMetadata(*dom);
if (metaitems != nullptr) {
for (char** item = metaitems; *item != nullptr; ++item) {
char* key;
const char* value = CPLParseNameValue(*item, &key);
metadata[std::string(*dom)][std::string(key)] = std::string(value);
VSIFree(key);
}
}
}
CSLDestroy(domains);
}
}
void GeoInfo::addTo(std::string const& filename) const {
if (!fs::exists(filename))
IF_THROW_EXCEPTION(not_found_error("Could not find any file at path " + filename + " to add GeoInfo to it."))
<< boost::errinfo_file_name(filename);
GDALDataset* img = (GDALDataset*) GDALOpen(filename.c_str(), GA_Update);
if (!img)
IF_THROW_EXCEPTION(runtime_error("The corresponding GDAL driver does not support update of metadata for this image type."));
addTo(img);
bool shouldHaveWrittenColorTable = !colorTable.empty() && img->GetRasterCount() == 1 && img->GetRasterBand(1)->GetRasterDataType() == GDT_Byte;
GDALClose(img);
// check for changed color table and warn
if (shouldHaveWrittenColorTable) {
GeoInfo test{filename};
compareColorTables(test, false);
}
}
void GeoInfo::addTo(GDALDataset* img) const {
if (!img)
IF_THROW_EXCEPTION(runtime_error("The GDALDataset to which the GeoInfo should be added is nullptr."));
// nodata values
if (!nodataValues.empty()) {
std::size_t rastercount = img->GetRasterCount();
for (std::size_t channel = 0; channel < rastercount; ++channel) {
std::size_t ndChan = std::min(channel, nodataValues.size() - 1);
double nodataVal = getNodataValue(ndChan);
if (!std::isnan(nodataVal))
img->GetRasterBand(channel+1)->SetNoDataValue(nodataVal);
}
}
// write color table only for uint8x1 images
// GDAL GTiff driver DOES CHANGE THE ALPHA CHANNEL TO 255 FOR ALL INDICES EXCEPT NODATA WHERE IT WILL BE CHANGED TO 0 ==> cannot write an arbitrary color table
if (!colorTable.empty() && img->GetRasterCount() == 1 && img->GetRasterBand(1)->GetRasterDataType() == GDT_Byte) {
GDALColorTable ct;
for (unsigned int i = 0; i < colorTable.size(); ++i) {
std::array<short, 4> const& col = colorTable.at(i);
GDALColorEntry* ce = const_cast<GDALColorEntry*>(reinterpret_cast<GDALColorEntry const*>(&col));
ct.SetColorEntry(i, ce);
}
img->GetRasterBand(1)->SetColorTable(&ct);
}
// geo transform (will only be used if it is not an identity transform)
bool useGeotransform = hasGeotransform();
if (useGeotransform) {
img->SetGeoTransform(const_cast<double*>(geotrans.values.data())); // const cast is fine, SetGeoTransform will not change geotransform
char* projStr;
geotransSRS.exportToWkt(&projStr);
if (std::strlen(projStr) > 0)
img->SetProjection(projStr);
CPLFree(projStr);
}
// ground control points (set only, if not using geotransform, since geotransform is more common)
bool useGCPs = !gcps.empty() && const_cast<OGRSpatialReference&>(gcpSRS).Validate() == OGRERR_NONE;
if (useGCPs && !useGeotransform) {
std::vector<GDAL_GCP> gdal_gcps;
for (GCP const& gcp : gcps) {
gdal_gcps.push_back(GDAL_GCP{const_cast<char*>(gcp.id.c_str()),
const_cast<char*>(gcp.info.c_str()),
gcp.pixel - 1, // workaround for incrementing dfGCPPixel on write
gcp.line - 1, // workaround for incrementing dfGCPLine on write
gcp.x,
gcp.y,
gcp.z});
}
if (!gdal_gcps.empty()) {
char* projStr;
gcpSRS.exportToWkt(&projStr);
img->SetGCPs(gdal_gcps.size(), gdal_gcps.data(), projStr);
CPLFree(projStr);
}
}
// remaining meta data
for (auto const& mp : metadata) {
const char* domain = mp.first.c_str();
for (auto const& item : mp.second) {
const char* key = item.first.c_str();
const char* val = item.second.c_str();
img->SetMetadataItem(key, val, domain);
}
}
}
bool GeoInfo::compareColorTables(GeoInfo const& other, bool quiet) const {
if (colorTable.empty())
return true;
if (!colorTable.empty() && other.colorTable.empty()) {
if (!quiet)
std::cerr << "Warning: The color table has been removed by the driver!" << std::endl;
return false;
}
if (colorTable.size() > other.colorTable.size()) {
if (!quiet)
std::cerr << "Warning: Not all entries of the color tables were written by the driver!" << std::endl;
return false;
}
// check which channels have been changed and if alpha has been changed to 255 except for nodata index
std::array<bool, 4> hasChanged = {false, false, false, false};
bool allNonNodataAlpha255Before = true;
bool allNonNodataAlpha255After = true;
for (unsigned int i = 0; i < colorTable.size(); ++i) {
std::array<short, 4> const& before = colorTable.at(i);
std::array<short, 4> const& after = other.colorTable.at(i);
if (!hasNodataValue() || i != getNodataValue())
if (before.back() != 255)
allNonNodataAlpha255Before = false;
if (!other.hasNodataValue() || i != other.getNodataValue())
if (after.back() != 255)
allNonNodataAlpha255After = false;
for (unsigned int c = 0; c < 4; ++c) {
hasChanged.at(c) = before.at(c) != after.at(c);
// if we don't need a warning we can quit directly
if (quiet && hasChanged.at(c))
return false;
}
}
bool nodataAlpha0Before = true;
bool nodataAlpha0After = true;
if (hasNodataValue() && colorTable.at(getNodataValue()).back() != 0)
nodataAlpha0Before = false;
if (other.hasNodataValue() && other.colorTable.at(other.getNodataValue()).back() != 0)
nodataAlpha0After = false;
// find and print reasons
bool areCompatible = true;
std::string warning = "Warning: The driver changed the color table. ";
for (unsigned int c = 0; c < 4; ++c) {
if (hasChanged.at(c)) {
if (areCompatible)
warning = "Channel(s) ";
warning += std::to_string(c) + ", ";
areCompatible = false;
}
}
if (!areCompatible) {
warning.pop_back(); warning.pop_back();
warning += " has/have changed. ";
}
if (!nodataAlpha0Before && nodataAlpha0After)
warning += "The nodata entry's alpha value has been changed to 0. ";
if (!allNonNodataAlpha255Before && allNonNodataAlpha255After)
warning += "All non-nodata entries' alpha values have been changed to 255.";
std::cerr << warning << std::endl;
return areCompatible;
}
GeoInfo GeoInfo::subdatasetGeoInfo(unsigned int idx) const {
if (!hasMetadataDomain("SUBDATASETS"))
IF_THROW_EXCEPTION(file_format_error("The file does not contain any subdatasets. "
"Hence, cannot collect the GeoInfo of subdataset "
+ std::to_string(idx)));
auto const& metaSDS = getMetadataItems("SUBDATASETS");
unsigned int numSDS = metaSDS.size() / 2; // for each subdataset there is a name and a description item
if (idx >= numSDS)
IF_THROW_EXCEPTION(invalid_argument_error("The file contains " + std::to_string(numSDS) + " subdatasets. "
"You tried to obtain subdataset " + std::to_string(idx) +
", which does not exist."));
std::string const& name = metaSDS.at("SUBDATASET_" + std::to_string(idx + 1) + "_NAME");
GDALDataset* sds = (GDALDataset*) GDALOpen(name.c_str(), GA_ReadOnly);
if (!sds)
IF_THROW_EXCEPTION(runtime_error("GDAL could not open the subdataset with index " + std::to_string(idx) + " and virtual filename " + name));
GeoInfo gi;
gi.readFrom(sds);
return gi;
}
GDALDataset* GeoInfo::openVrtGdalDataset(std::vector<int> const& selChans, InterpMethod interp) const {
// collect subdataset names
auto sds = subdatasets();
if (sds.empty())
IF_THROW_EXCEPTION(image_type_error("The file does not contain any subdatasets."));
char** papszSrcDSNames = nullptr;
int bandCount = selChans.empty() ? sds.size() : selChans.size();
if (selChans.empty())
for (unsigned int c = 0; c < sds.size(); ++c)
papszSrcDSNames = CSLAddString(papszSrcDSNames, sds.at(c).first.c_str());
else {
if (selChans.size() == 1) {
int c = selChans.front();
if (c >= static_cast<int>(sds.size()) || c < 0)
IF_THROW_EXCEPTION(image_type_error("You acquired subdataset " + std::to_string(c)
+ " that does not exist. The image only has "
+ std::to_string(sds.size()) + " subdatasets."));
return (GDALDataset*) GDALOpen(sds.at(c).first.c_str(), GA_ReadOnly);
}
for (int c : selChans) {
if (c >= static_cast<int>(sds.size()) || c < 0)
IF_THROW_EXCEPTION(image_type_error("You acquired subdataset " + std::to_string(c)
+ " that does not exist. The image only has "
+ std::to_string(sds.size()) + " subdatasets."));
papszSrcDSNames = CSLAddString(papszSrcDSNames, sds.at(c).first.c_str());
}
}
// set some default options
char** papszVRTopts = nullptr;
papszVRTopts = CSLAddString(papszVRTopts, "-resolution");
papszVRTopts = CSLAddString(papszVRTopts, "highest");
papszVRTopts = CSLAddString(papszVRTopts, "-separate");
papszVRTopts = CSLAddString(papszVRTopts, "-r");
if (interp == InterpMethod::nearest)
papszVRTopts = CSLAddString(papszVRTopts, "nearest");
else if (interp == InterpMethod::bilinear)
papszVRTopts = CSLAddString(papszVRTopts, "bilinear");
else if (interp == InterpMethod::cubic)
papszVRTopts = CSLAddString(papszVRTopts, "cubic");
else // interp == InterpMethod::cubicspline
papszVRTopts = CSLAddString(papszVRTopts, "cubicspline");
GDALBuildVRTOptions* vrtOpts = GDALBuildVRTOptionsNew(papszVRTopts, nullptr);
CSLDestroy(papszVRTopts);
int bUsageError = FALSE;
GDALDataset* poImg = (GDALDataset*) GDALBuildVRT(/* pszDest */ nullptr, bandCount, /* pahSrcDS */ nullptr, papszSrcDSNames, vrtOpts, &bUsageError);
CSLDestroy(papszSrcDSNames);
GDALBuildVRTOptionsFree(vrtOpts);
if (bUsageError) {
GDALClose(poImg);
IF_THROW_EXCEPTION(runtime_error("Some error occurred while combinding the subdatasets. Possible reasons "
"include they have different data types, projection coordinate systems, etc. "
"Maybe there is another error or warning message from GDAL, which known more."));
}
return poImg;
}
GeoInfo::GeoInfo(std::string const& filename, std::vector<int> const& selChans, Rectangle crop, bool flipH, bool flipV, bool recurseSubdatasets)
: GeoInfo()
{
readFrom(filename);
bool doRecurseSubdatasets = hasSubdatasets() && (recurseSubdatasets || !selChans.empty());
if (doRecurseSubdatasets) {
channels = subdatasetsCount();
}
if (!selChans.empty()) {
for (int const& c : selChans)
if (c >= static_cast<int>(channels) || c < 0)
IF_THROW_EXCEPTION(image_type_error("You acquired a channel (" + std::to_string(c)
+ ") that does not exist. The image only has "
+ std::to_string(channels) + " channels."))
<< boost::errinfo_file_name(filename);
channels = selChans.size();
}
if (doRecurseSubdatasets) {
GDALDataset* poVrtImg = openVrtGdalDataset(selChans);
*this = GeoInfo{};
readFrom(poVrtImg);
GDALClose(poVrtImg);
}
crop.x = std::min(std::max(crop.x, 0), width());
crop.y = std::min(std::max(crop.y, 0), height());
if (crop.width == 0)
crop.width = width() - crop.x;
if (crop.height == 0)
crop.height = height() - crop.y;
size.width = crop.width;
size.height = crop.height;
if (hasGeotransform()) {
geotrans.translateImage(crop.x, crop.y);
geotrans.flipImage(flipH, flipV, size);
}
}
CoordRectangle intersectRect(GeoInfo const& refGI, CoordRectangle const& refRect, GeoInfo const& otherGI, CoordRectangle const& otherRect, unsigned int numPoints) {
if (refRect.width <= 0 || refRect.height <= 0 || otherRect.width <= 0 || otherRect.height <= 0 )
return {};
if (numPoints <= 0)
return {};
// collect source boundaries
std::vector<Coordinate> boundariesRef = detail::makeRectBoundaryCoords(refRect, numPoints);
// transform to other projection coordinate space and restrict by other rectangle
std::vector<Coordinate> boundariesOther = refGI.projToProj(boundariesRef, otherGI);
Coordinate otherBR = otherRect.br();
for (Coordinate& c : boundariesOther) {
// note: this corresponds to a projection onto the otherRect boundary and can give points that are out of refRect!
c.x = std::min(std::max(c.x, otherRect.x), otherBR.x);
c.y = std::min(std::max(c.y, otherRect.y), otherBR.y);
}
// transform back
boundariesRef = otherGI.projToProj(boundariesOther, refGI);
// intersect (enclose all points by a rectangle)
CoordRectangle projectedLimitation = detail::getRectFromBoundaryCoords(boundariesRef);
// the new bounds can be larger than refRect, so take refRect as limit
return projectedLimitation & refRect;
}
void GeoInfo::intersect(GeoInfo const& other, unsigned int numCoords, bool shrink) {
// find intersection
CoordRectangle inter = intersectRect(*this, projRect(), other, other.projRect(), numCoords);
if (inter.area() <= 0) {
size = Size{};
return;
}
setExtents(inter, shrink);
}
void GeoInfo::setExtents(CoordRectangle const& ex, bool shrink) {
assert(geotrans.XToY == 0 && geotrans.YToX == 0 && "Only simple transformations supported.");
// get top left and bottom right corners
Coordinate projCornerTL = ex.tl();
Coordinate projCornerBR = ex.br();
if (geotrans.XToX < 0)
std::swap(projCornerTL.x, projCornerBR.x);
if (geotrans.YToY < 0)
std::swap(projCornerTL.y, projCornerBR.y);
Coordinate imgCornerTL = geotrans.projToImg(projCornerTL);
Coordinate imgCornerBR = geotrans.projToImg(projCornerBR);
// find new size
constexpr double abstol = 1e-11;
double newWidth = imgCornerBR.x - imgCornerTL.x;
double newHeight = imgCornerBR.y - imgCornerTL.y;
if (shrink) {
newWidth = std::floor(newWidth + abstol);
newHeight = std::floor(newHeight + abstol);
}
else {
newWidth = std::ceil(newWidth - abstol);
newHeight = std::ceil(newHeight - abstol);
}
// if one corner did not change, try to preserve it, if all changed take top left
bool alignLeftX = false; // otherwise to right x
if (std::abs(imgCornerTL.x) < abstol || std::abs(imgCornerBR.x - width()) >= abstol)
alignLeftX = true; // default, when both changed
else
imgCornerTL.x = imgCornerBR.x - newWidth;
bool alignTopY = false; // otherwise to bottom y
if (std::abs(imgCornerTL.y) < abstol || std::abs(imgCornerBR.y - height()) >= abstol)
alignTopY = true; // default, when both changed
else
imgCornerTL.y = imgCornerBR.y - newHeight;
if (!alignLeftX || !alignTopY)
projCornerTL = geotrans.imgToProj(imgCornerTL);
// save extents
size.width = newWidth;
size.height = newHeight;
geotrans.offsetX = projCornerTL.x;
geotrans.offsetY = projCornerTL.y;
}
GeoTransform& GeoTransform::clear() {
offsetX = 0;
offsetY = 0;
XToX = 1;
YToX = 0;
XToY = 0;
YToY = 1;
return *this;
}
Coordinate GeoInfo::imgToProj(Coordinate const& c_i, GeoInfo const& to) const {
return imgToProj(std::vector<Coordinate>{c_i}, to).front();
}
Coordinate GeoInfo::projToImg(Coordinate const& c_p, GeoInfo const& to) const {
return projToImg(std::vector<Coordinate>{c_p}, to).front();
}
Coordinate GeoInfo::imgToImg(Coordinate const& c_i, GeoInfo const& to) const {
return imgToImg(std::vector<Coordinate>{c_i}, to).front();
}
Coordinate GeoInfo::projToProj(Coordinate const& c_p, GeoInfo const& to) const {
return projToProj(std::vector<Coordinate>{c_p}, to).front();
}
std::vector<Coordinate> GeoInfo::imgToProj(std::vector<Coordinate> const& c_i, GeoInfo const& to) const {
return convertProj(*this, to, c_i, true, false);
}
std::vector<Coordinate> GeoInfo::projToImg(std::vector<Coordinate> const& c_p, GeoInfo const& to) const {
return convertProj(*this, to, c_p, false, true);
}
std::vector<Coordinate> GeoInfo::imgToImg(std::vector<Coordinate> const& c_i, GeoInfo const& to) const {
return convertProj(*this, to, c_i, true, true);
}
std::vector<Coordinate> GeoInfo::projToProj(std::vector<Coordinate> const& c_p, GeoInfo const& to) const {
return convertProj(*this, to, c_p, false, false);
}
Coordinate GeoInfo::imgToLongLat(Coordinate const& c_i) const {
return imgToLongLat(std::vector<Coordinate>{c_i}).front();
}
Coordinate GeoInfo::projToLongLat(Coordinate const& c_p) const {
return projToLongLat(std::vector<Coordinate>{c_p}).front();
}
Coordinate GeoInfo::longLatToProj(Coordinate const& c_l) const {
return longLatToProj(std::vector<Coordinate>{c_l}).front();
}
Coordinate GeoInfo::longLatToImg(Coordinate const& c_l) const {
return longLatToImg(std::vector<Coordinate>{c_l}).front();
}
std::vector<Coordinate> GeoInfo::imgToLongLat(std::vector<Coordinate> const& c_i) const {
return convertLongLat(*this, c_i, true, false);
}
std::vector<Coordinate> GeoInfo::projToLongLat(std::vector<Coordinate> const& c_p) const {
return convertLongLat(*this, c_p, false, false);
}
std::vector<Coordinate> GeoInfo::longLatToProj(std::vector<Coordinate> const& c_l) const {
return convertLongLat(*this, c_l, false, true);
}
std::vector<Coordinate> GeoInfo::longLatToImg(std::vector<Coordinate> const& c_l) const {
return convertLongLat(*this, c_l, true, true);
}
Coordinate GeoTransform::imgToProj(Coordinate const& c_i) const {
return Coordinate{
offsetX + XToX * (c_i.x) + YToX * (c_i.y),
offsetY + XToY * (c_i.x) + YToY * (c_i.y)
};
}
Coordinate GeoTransform::projToImg(Coordinate const& c_p) const {
double const x = c_p.x - offsetX;
double const y = c_p.y - offsetY;
double const det = XToX * YToY - XToY * YToX;
return Coordinate{
( YToY * x - YToX * y) / det,
(-XToY * x + XToX * y) / det
};
}
CoordRectangle GeoTransform::imgToProj(CoordRectangle const& r_i) const {
Coordinate projTL = imgToProj(r_i.tl());
Coordinate projBR = imgToProj(r_i.br());
return CoordRectangle{std::min(projTL.x, projBR.x), std::min(projTL.y, projBR.y),
std::abs(projTL.x - projBR.x), std::abs(projTL.y - projBR.y)};
}
CoordRectangle GeoTransform::projToImg(CoordRectangle const& r_p) const {
Coordinate imgTL = projToImg(r_p.tl());
Coordinate imgBR = projToImg(r_p.br());
return CoordRectangle{std::min(imgTL.x, imgBR.x), std::min(imgTL.y, imgBR.y),
std::abs(imgTL.x - imgBR.x), std::abs(imgTL.y - imgBR.y)};
}
GeoTransform& GeoTransform::scaleProjection(double xscale, double yscale) {
offsetX *= xscale;
offsetY *= yscale;
XToX *= xscale;
YToX *= xscale;
XToY *= yscale;
YToY *= yscale;
return *this;
}
GeoTransform& GeoTransform::scaleImage(double xscale, double yscale) {
XToX *= xscale;
YToX *= yscale;
XToY *= xscale;
YToY *= yscale;
return *this;
}
GeoTransform& GeoTransform::rotateProjection(double angle) {
double cos = std::cos(angle * M_PI / 180);
double sin = std::sin(angle * M_PI / 180);
auto rotate = [cos,sin] (double xin, double yin) {
return std::array<double,2>{
cos * xin - sin * yin, // xout
sin * xin + cos * yin // yout
};
};
constexpr unsigned int X = 0;
constexpr unsigned int Y = 1;
std::array<double,2> rotated;
// rotate origin
rotated = rotate(offsetX, offsetY);
offsetX = rotated[X];
offsetY = rotated[Y];
// rotate first column of A
rotated = rotate(XToX, XToY);
XToX = rotated[X];
XToY = rotated[Y];
// rotate second column of A
rotated = rotate(YToX, YToY);
YToX = rotated[X];
YToY = rotated[Y];
return *this;
}
GeoTransform& GeoTransform::rotateImage(double angle) {
double cos = std::cos(angle * M_PI / 180);
double sin = std::sin(angle * M_PI / 180);
// get first column of A
double temp_XToX = XToX * cos + YToX * sin;
double temp_XToY = XToY * cos + YToY * sin;
// get second column of A
double temp_YToX = -XToX * sin + YToX * cos;
double temp_YToY = -XToY * sin + YToY * cos;
XToX = temp_XToX;
XToY = temp_XToY;
YToX = temp_YToX;
YToY = temp_YToY;
return *this;
}
GeoTransform& GeoTransform::flipImage(bool flipH, bool flipV, Size s) {
translateImage(flipH ? s.width : 0,
flipV ? s.height : 0);
scaleImage(flipH ? -1 : 1,
flipV ? -1 : 1);
return *this;
}
GeoTransform& GeoTransform::shearXProjection(double factor) {
offsetX += offsetY * factor;
XToX += XToY * factor;
YToX += YToY * factor;
return *this;
}
GeoTransform& GeoTransform::shearXImage(double factor) {
YToX += XToX * factor;
YToY += XToY * factor;
return *this;
}
GeoTransform& GeoTransform::shearYProjection(double factor) {
offsetY += offsetX * factor;
XToY += XToX * factor;
YToY += YToX * factor;
return *this;
}
GeoTransform& GeoTransform::shearYImage(double factor) {
XToX += YToX * factor;
XToY += YToY * factor;
return *this;
}
GeoTransform& GeoTransform::translateProjection(double xoffset, double yoffset) {
offsetX += xoffset;
offsetY += yoffset;
return *this;
}
GeoTransform& GeoTransform::translateImage(double xoffset, double yoffset) {
offsetX += XToX * xoffset + YToX * yoffset;
offsetY += XToY * xoffset + YToY * yoffset;
return *this;
}
void GeoTransform::set(double topLeft_x, double topLeft_y, double x_to_x, double y_to_x, double x_to_y, double y_to_y) {
offsetX = topLeft_x;
offsetY = topLeft_y;
XToX = x_to_x;
YToY = y_to_y;
YToX = y_to_x;
XToY = x_to_y;
}
bool GeoTransform::isIdentity() const {
return offsetX == 0
&& offsetY == 0
&& XToX == 1
&& YToX == 0
&& XToY == 0
&& YToY == 1;
}
bool GeoInfo::hasGeotransform() const {
return !geotrans.isIdentity() && const_cast<OGRSpatialReference&>(geotransSRS).Validate() == OGRERR_NONE;
}
bool GeoInfo::hasGCPs() const {
return gcps.size() >= 3 && const_cast<OGRSpatialReference&>(gcpSRS).Validate() == OGRERR_NONE;
}
std::map<std::string, std::string> const& GeoInfo::getMetadataItems(std::string const& domain) const {
auto dom_it = metadata.find(domain);
if (dom_it != metadata.end())
return dom_it->second;
IF_THROW_EXCEPTION(not_found_error("Could not find the requested metadata domain " + domain + "!"));
}
double GeoInfo::getNodataValue(unsigned int channel) const {
return nodataValues.at(std::min(static_cast<std::size_t>(channel), nodataValues.size() - 1));
}
} /* namespace imagefusion */