1 /*
2 * Copyright 2012-2020 Kai Pastor
3 *
4 * This file is part of OpenOrienteering.
5 *
6 * OpenOrienteering is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * OpenOrienteering is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with OpenOrienteering. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20
21 #include "georeferencing.h"
22
23 #include <algorithm>
24 #include <cmath> // IWYU pragma: keep
25 #include <iterator>
26 #include <utility>
27
28 #include <QtGlobal>
29 #include <QtMath>
30 #include <QByteArray>
31 #include <QDebug>
32 #include <QDir> // IWYU pragma: keep
33 #include <QFileInfo>
34 #include <QLatin1String>
35 #include <QLocale>
36 #include <QPoint>
37 #include <QSignalBlocker>
38 #include <QStandardPaths> // IWYU pragma: keep
39 #include <QStringRef>
40 #include <QXmlStreamReader>
41 #include <QXmlStreamWriter>
42
43 #ifdef ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
44 # include <proj_api.h>
45 #else
46 # include <proj.h>
47 #endif
48
49 #include "core/crs_template.h"
50 #include "fileformats/file_format.h"
51 #include "fileformats/xml_file_format.h"
52 #include "util/xml_stream_util.h"
53
54
55 // ### A namespace which collects various string constants of type QLatin1String. ###
56
57 namespace literal
58 {
59 static const QLatin1String georeferencing("georeferencing");
60
61 static const QLatin1String scale("scale");
62 static const QLatin1String grid_scale_factor{"grid_scale_factor"};
63 static const QLatin1String auxiliary_scale_factor{"auxiliary_scale_factor"};
64 static const QLatin1String declination("declination");
65 static const QLatin1String grivation("grivation");
66
67 static const QLatin1String ref_point("ref_point");
68 static const QLatin1String ref_point_deg("ref_point_deg");
69 static const QLatin1String projected_crs("projected_crs");
70 static const QLatin1String geographic_crs("geographic_crs");
71 static const QLatin1String spec("spec");
72 static const QLatin1String parameter("parameter");
73
74 static const QLatin1String x("x");
75 static const QLatin1String y("y");
76
77 static const QLatin1String id("id");
78 static const QLatin1String language("language");
79
80 static const QLatin1String lat("lat");
81 static const QLatin1String lon("lon");
82
83 static const QLatin1String proj_4("PROJ.4");
84 static const QLatin1String geographic_coordinates("Geographic coordinates");
85 }
86
87
88
89 namespace OpenOrienteering {
90
91 namespace
92 {
93 #ifdef Q_OS_ANDROID
94 /**
95 * Provides a cache directory for RROJ resource files.
96 */
projCacheDirectory()97 const QDir& projCacheDirectory()
98 {
99 static auto const dir = []() -> QDir {
100 auto cache = QDir(QStandardPaths::writableLocation(QStandardPaths::CacheLocation));
101 auto proj = QStringLiteral("proj");
102 if (!cache.exists(proj))
103 cache.mkpath(proj);
104 if (!cache.exists(proj))
105 qDebug("Could not create a cache directory for PROJ resource files");
106 cache.cd(proj);
107 return cache;
108 }();
109 return dir;
110 }
111
112 /**
113 * Ensures PROJ resource file availability in the cache directory.
114 *
115 * This function implements the interface required by
116 * proj_context_set_file_finder().
117 *
118 * This functions checks if the requested file name is available in the
119 * cache directory. If not, it tries to copy the file from the proj
120 * subfolder of the assets folder to this temporary directory.
121 *
122 * It always returns nullptr, letting PROJ find the file in the cache
123 * directory via the configured search path.
124 */
125 extern "C"
projFileHelperAndroid(PJ_CONTEXT *,const char * name,void *)126 const char* projFileHelperAndroid(PJ_CONTEXT* /*ctx*/, const char* name, void* /*user_data*/)
127 {
128 auto const& cache = projCacheDirectory();
129 if (cache.exists())
130 {
131 auto const name_string = QString::fromUtf8(name);
132 auto const path = cache.filePath(name_string);
133 if (!QFile::exists(path)
134 && !QFile::copy(QLatin1String("assets:/proj/") + name_string, path))
135 {
136 qDebug("Could not provide projection data file '%s'", name);
137 }
138 }
139 return nullptr;
140 }
141 #endif
142
143
144 /** Helper for PROJ initialization.
145 *
146 * To be used as static object in the right place.
147 */
148 class ProjSetup
149 {
150 public:
ProjSetup()151 ProjSetup()
152 {
153 #ifdef ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
154 auto proj_data = QFileInfo(QLatin1String("data:/proj"));
155 if (proj_data.exists())
156 {
157 static auto const location = proj_data.absoluteFilePath().toLocal8Bit();
158 static auto* data = location.constData();
159 pj_set_searchpath(1, &data);
160 }
161
162 // no Android file finder support for deprecated PROJ API
163 #else
164 proj_context_use_proj4_init_rules(PJ_DEFAULT_CTX, 1);
165
166 #if defined(Q_OS_ANDROID)
167 // Register file finder function needed by Proj.4
168 proj_context_set_file_finder(nullptr, &projFileHelperAndroid, nullptr);
169 auto proj_data = QFileInfo(projCacheDirectory().path());
170 #else
171 auto proj_data = QFileInfo(QLatin1String("data:/proj"));
172 #endif // defined(Q_OS_ANDROID)
173 if (proj_data.exists())
174 {
175 static auto const location = proj_data.absoluteFilePath().toLocal8Bit();
176 static auto* const data = location.constData();
177 proj_context_set_search_paths(nullptr, 1, &data);
178 }
179 #endif // ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
180 }
181 };
182
183
184 /**
185 * List of substitutions for specifications which are known to be broken in PROJ.
186 */
187 const char* spec_substitutions[][2] = {
188 // #542, S-JTSK (Greenwich) / Krovak East North
189 { "+init=epsg:5514", "+proj=krovak"
190 " +lat_0=49.5 +lon_0=24.83333333333333"
191 " +alpha=30.28813972222222 +k=0.9999"
192 " +x_0=0 +y_0=0"
193 " +ellps=bessel"
194 " +pm=greenwich"
195 " +towgs84=542.5,89.2,456.9,5.517,2.275,5.516,6.96"
196 " +units=m"
197 " +no_defs"
198 },
199 };
200
201
202 } // namespace
203
204
205
206 //### ProjTransform ###
207
208 #ifdef ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
209
ProjTransform(ProjTransformData * pj)210 ProjTransform::ProjTransform(ProjTransformData* pj) noexcept
211 : pj{pj}
212 {}
213
ProjTransform(ProjTransform && other)214 ProjTransform::ProjTransform(ProjTransform&& other) noexcept
215 {
216 operator=(std::move(other));
217 }
218
ProjTransform(const QString & crs_spec)219 ProjTransform::ProjTransform(const QString& crs_spec)
220 {
221 auto spec_latin1 = crs_spec.toLatin1();
222 if (!spec_latin1.contains("+no_defs"))
223 spec_latin1.append(" +no_defs");
224
225 *pj_get_errno_ref() = 0;
226 pj = pj_init_plus(spec_latin1);
227 }
228
~ProjTransform()229 ProjTransform::~ProjTransform()
230 {
231 if (pj)
232 pj_free(pj);
233 }
234
operator =(ProjTransform && other)235 ProjTransform& ProjTransform::operator=(ProjTransform&& other) noexcept
236 {
237 std::swap(pj, other.pj);
238 return *this;
239 }
240
isValid() const241 bool ProjTransform::isValid() const noexcept
242 {
243 return pj != nullptr;
244 }
245
isGeographic() const246 bool ProjTransform::isGeographic() const
247 {
248 return isValid() && pj_is_latlong(pj);
249 }
250
forward(const LatLon & lat_lon,bool * ok) const251 QPointF ProjTransform::forward(const LatLon& lat_lon, bool* ok) const
252 {
253 static auto const geographic_crs = ProjTransform(Georeferencing::geographic_crs_spec);
254
255 double easting = qDegreesToRadians(lat_lon.longitude()), northing = qDegreesToRadians(lat_lon.latitude());
256 if (geographic_crs.isValid())
257 {
258 auto ret = pj_transform(geographic_crs.pj, pj, 1, 1, &easting, &northing, nullptr);
259 if (ok)
260 *ok = (ret == 0);
261 }
262 else if (ok)
263 {
264 *ok = false;
265 }
266 return {easting, northing};
267 }
268
inverse(const QPointF & projected_coords,bool * ok) const269 LatLon ProjTransform::inverse(const QPointF& projected_coords, bool* ok) const
270 {
271 static auto const geographic_crs = ProjTransform(Georeferencing::geographic_crs_spec);
272
273 double easting = projected_coords.x(), northing = projected_coords.y();
274 if (geographic_crs.isValid())
275 {
276 auto ret = pj_transform(pj, geographic_crs.pj, 1, 1, &easting, &northing, nullptr);
277 if (ok)
278 *ok = (ret == 0);
279 }
280 else if (ok)
281 {
282 *ok = false;
283 }
284 return LatLon::fromRadiant(northing, easting);
285 }
286
errorText() const287 QString ProjTransform::errorText() const
288 {
289 auto err_no = *pj_get_errno_ref();
290 return (err_no == 0) ? QString() : QString::fromLatin1(pj_strerrno(err_no));
291 }
292
293 #else
294
ProjTransform(ProjTransformData * pj)295 ProjTransform::ProjTransform(ProjTransformData* pj) noexcept
296 : pj{pj}
297 {}
298
ProjTransform(ProjTransform && other)299 ProjTransform::ProjTransform(ProjTransform&& other) noexcept
300 {
301 operator=(std::move(other));
302 }
303
ProjTransform(const QString & crs_spec)304 ProjTransform::ProjTransform(const QString& crs_spec)
305 {
306 auto spec_latin1 = crs_spec.toLatin1();
307 #ifdef PROJ_ISSUE_1573
308 // Cf. https://github.com/OSGeo/PROJ/pull/1573
309 spec_latin1.replace("+datum=potsdam", "+ellps=bessel +nadgrids=@BETA2007.gsb");
310 #endif
311 pj = proj_create_crs_to_crs(PJ_DEFAULT_CTX, Georeferencing::geographic_crs_spec.toLatin1(), spec_latin1, nullptr);
312 if (pj)
313 operator=({proj_normalize_for_visualization(PJ_DEFAULT_CTX, pj)});
314 }
315
~ProjTransform()316 ProjTransform::~ProjTransform()
317 {
318 if (pj)
319 proj_destroy(pj);
320 }
321
operator =(ProjTransform && other)322 ProjTransform& ProjTransform::operator=(ProjTransform&& other) noexcept
323 {
324 std::swap(pj, other.pj);
325 return *this;
326 }
327
isValid() const328 bool ProjTransform::isValid() const noexcept
329 {
330 return pj != nullptr;
331 }
332
isGeographic() const333 bool ProjTransform::isGeographic() const
334 {
335 /// \todo Evaluate proj_get_type() instead
336 return isValid() && proj_angular_output(pj, PJ_FWD);
337 }
338
forward(const LatLon & lat_lon,bool * ok) const339 QPointF ProjTransform::forward(const LatLon& lat_lon, bool* ok) const
340 {
341 proj_errno_reset(pj);
342 auto pj_coord = proj_trans(pj, PJ_FWD, proj_coord(lat_lon.longitude(), lat_lon.latitude(), 0, HUGE_VAL));
343 if (ok)
344 *ok = proj_errno(pj) == 0;
345 return {pj_coord.xy.x, pj_coord.xy.y};
346 }
347
inverse(const QPointF & projected_coords,bool * ok) const348 LatLon ProjTransform::inverse(const QPointF& projected_coords, bool* ok) const
349 {
350 proj_errno_reset(pj);
351 auto pj_coord = proj_trans(pj, PJ_INV, proj_coord(projected_coords.x(), projected_coords.y(), 0, HUGE_VAL));
352 if (ok)
353 *ok = proj_errno(pj) == 0;
354 return {pj_coord.lp.phi, pj_coord.lp.lam};
355 }
356
errorText() const357 QString ProjTransform::errorText() const
358 {
359 auto err_no = proj_errno(pj);
360 return (err_no == 0) ? QString() : QString::fromLatin1(proj_errno_string(err_no));
361 }
362
363 #endif
364
365
366
367 //### Georeferencing ###
368
369 const QString Georeferencing::geographic_crs_spec(QString::fromLatin1("+proj=latlong +datum=WGS84"));
370
Georeferencing()371 Georeferencing::Georeferencing()
372 : state(Local),
373 scale_denominator{1000},
374 combined_scale_factor{1.0},
375 auxiliary_scale_factor{1.0},
376 grid_scale_factor{1.0},
377 declination(0.0),
378 grivation(0.0),
379 grivation_error(0.0),
380 convergence(0.0),
381 map_ref_point(0, 0),
382 projected_ref_point(0, 0)
383 {
384 static ProjSetup run_once;
385
386 updateTransformation();
387
388 projected_crs_id = QString::fromLatin1("Local");
389 }
390
Georeferencing(const Georeferencing & other)391 Georeferencing::Georeferencing(const Georeferencing& other)
392 : QObject(),
393 state(other.state),
394 scale_denominator(other.scale_denominator),
395 combined_scale_factor{other.combined_scale_factor},
396 auxiliary_scale_factor{other.auxiliary_scale_factor},
397 grid_scale_factor(other.grid_scale_factor),
398 declination(other.declination),
399 grivation(other.grivation),
400 grivation_error(other.grivation_error),
401 convergence(other.convergence),
402 map_ref_point(other.map_ref_point),
403 projected_ref_point(other.projected_ref_point),
404 projected_crs_id(other.projected_crs_id),
405 projected_crs_spec(other.projected_crs_spec),
406 projected_crs_parameters(other.projected_crs_parameters),
407 proj_transform(projected_crs_spec),
408 geographic_ref_point(other.geographic_ref_point)
409 {
410 updateTransformation();
411 }
412
413 Georeferencing::~Georeferencing() = default;
414
operator =(const Georeferencing & other)415 Georeferencing& Georeferencing::operator=(const Georeferencing& other)
416 {
417 if (&other == this)
418 return *this;
419
420 state = other.state;
421 scale_denominator = other.scale_denominator;
422 combined_scale_factor = other.combined_scale_factor;
423 auxiliary_scale_factor = other.auxiliary_scale_factor;
424 grid_scale_factor = other.grid_scale_factor;
425 declination = other.declination;
426 grivation = other.grivation;
427 grivation_error = other.grivation_error;
428 convergence = other.convergence;
429 map_ref_point = other.map_ref_point;
430 projected_ref_point = other.projected_ref_point;
431 from_projected = other.from_projected;
432 to_projected = other.to_projected;
433 projected_ref_point = other.projected_ref_point;
434 projected_crs_id = other.projected_crs_id;
435 projected_crs_spec = other.projected_crs_spec;
436 projected_crs_parameters = other.projected_crs_parameters;
437 proj_transform = ProjTransform(other.projected_crs_spec);
438 geographic_ref_point = other.geographic_ref_point;
439
440 emit stateChanged();
441 emit transformationChanged();
442 emit declinationChanged();
443 emit auxiliaryScaleFactorChanged();
444 emit projectionChanged();
445
446 return *this;
447 }
448
449
450
isGeographic() const451 bool Georeferencing::isGeographic() const
452 {
453 return proj_transform.isGeographic();
454 }
455
456
457
load(QXmlStreamReader & xml,bool load_scale_only)458 void Georeferencing::load(QXmlStreamReader& xml, bool load_scale_only)
459 {
460 Q_ASSERT(xml.name() == literal::georeferencing);
461
462 {
463 // Reset to default values
464 const QSignalBlocker block(this);
465 *this = Georeferencing();
466 }
467
468 XmlElementReader georef_element(xml);
469 scale_denominator = georef_element.attribute<int>(literal::scale);
470 if (scale_denominator <= 0)
471 throw FileFormatException(tr("Map scale specification invalid or missing."));
472
473 if (georef_element.hasAttribute(literal::grid_scale_factor))
474 {
475 combined_scale_factor = roundScaleFactor(georef_element.attribute<double>(literal::grid_scale_factor));
476 if (combined_scale_factor <= 0.0)
477 throw FileFormatException(tr("Invalid grid scale factor: %1").arg(QString::number(combined_scale_factor)));
478 }
479 state = Local;
480 if (load_scale_only)
481 {
482 xml.skipCurrentElement();
483 }
484 else
485 {
486 if (georef_element.hasAttribute(literal::auxiliary_scale_factor))
487 {
488 auxiliary_scale_factor = roundScaleFactor(georef_element.attribute<double>(literal::auxiliary_scale_factor));
489 if (auxiliary_scale_factor <= 0.0)
490 throw FileFormatException(tr("Invalid auxiliary scale factor: %1").arg(QString::number(auxiliary_scale_factor)));
491 }
492 if (georef_element.hasAttribute(literal::declination))
493 declination = roundDeclination(georef_element.attribute<double>(literal::declination));
494 if (georef_element.hasAttribute(literal::grivation))
495 {
496 grivation = roundDeclination(georef_element.attribute<double>(literal::grivation));
497 grivation_error = georef_element.attribute<double>(literal::grivation) - grivation;
498 }
499
500 while (xml.readNextStartElement())
501 {
502 if (xml.name() == literal::ref_point)
503 {
504 XmlElementReader ref_point_element(xml);
505 map_ref_point.setX(ref_point_element.attribute<double>(literal::x));
506 map_ref_point.setY(ref_point_element.attribute<double>(literal::y));
507 }
508 else if (xml.name() == literal::projected_crs)
509 {
510 XmlElementReader crs_element(xml);
511 projected_crs_id = crs_element.attribute<QString>(literal::id);
512 while (xml.readNextStartElement())
513 {
514 XmlElementReader current_element(xml);
515 if (xml.name() == literal::spec)
516 {
517 const QString language = current_element.attribute<QString>(literal::language);
518 if (language != literal::proj_4)
519 throw FileFormatException(tr("Unknown CRS specification language: %1").arg(language));
520 projected_crs_spec = xml.readElementText();
521 }
522 else if (xml.name() == literal::parameter)
523 {
524 projected_crs_parameters.push_back(xml.readElementText());
525 }
526 else if (xml.name() == literal::ref_point)
527 {
528 projected_ref_point.setX(current_element.attribute<double>(literal::x));
529 projected_ref_point.setY(current_element.attribute<double>(literal::y));
530 }
531 else // unknown
532 {
533 ; // nothing
534 }
535 }
536 }
537 else if (xml.name() == literal::geographic_crs)
538 {
539 while (xml.readNextStartElement())
540 {
541 XmlElementReader current_element(xml);
542 if (xml.name() == literal::spec)
543 {
544 const QString language = current_element.attribute<QString>(literal::language);
545 if (language != literal::proj_4)
546 throw FileFormatException(tr("Unknown CRS specification language: %1").arg(language));
547 QString geographic_crs_spec = xml.readElementText();
548 if (Georeferencing::geographic_crs_spec != geographic_crs_spec)
549 throw FileFormatException(tr("Unsupported geographic CRS specification: %1").arg(geographic_crs_spec));
550 }
551 else if (xml.name() == literal::ref_point)
552 {
553 // Legacy, latitude/longitude in radiant
554 double latitude = current_element.attribute<double>(literal::lat);
555 double longitude = current_element.attribute<double>(literal::lon);
556 geographic_ref_point = LatLon::fromRadiant(latitude, longitude);
557 }
558 else if (xml.name() == literal::ref_point_deg)
559 {
560 // Legacy, latitude/longitude in degrees
561 double latitude = current_element.attribute<double>(literal::lat);
562 double longitude = current_element.attribute<double>(literal::lon);
563 geographic_ref_point = LatLon(latitude, longitude);
564 }
565 else // unknown
566 {
567 ; // nothing
568 }
569 }
570 }
571 else // unknown
572 {
573 xml.skipCurrentElement();
574 }
575 }
576 }
577
578 emit stateChanged();
579 updateTransformation();
580 emit declinationChanged();
581 if (!projected_crs_spec.isEmpty())
582 {
583 proj_transform = {projected_crs_spec};
584 if (proj_transform.isValid())
585 {
586 state = Normal;
587 }
588 updateGridCompensation();
589 if (!georef_element.hasAttribute(literal::auxiliary_scale_factor))
590 {
591 initAuxiliaryScaleFactor();
592 }
593 emit auxiliaryScaleFactorChanged();
594 if (proj_transform.isValid())
595 {
596 emit stateChanged();
597 }
598 }
599 emit projectionChanged();
600 }
601
save(QXmlStreamWriter & xml) const602 void Georeferencing::save(QXmlStreamWriter& xml) const
603 {
604 XmlElementWriter georef_element(xml, literal::georeferencing);
605 georef_element.writeAttribute(literal::scale, scale_denominator);
606 if (combined_scale_factor != 1.0 || auxiliary_scale_factor != 1.0)
607 {
608 if (combined_scale_factor != 1.0)
609 georef_element.writeAttribute(literal::grid_scale_factor, combined_scale_factor, scaleFactorPrecision());
610 georef_element.writeAttribute(literal::auxiliary_scale_factor, auxiliary_scale_factor, scaleFactorPrecision());
611 }
612 if (!qIsNull(declination))
613 georef_element.writeAttribute(literal::declination, declination, declinationPrecision());
614 if (!qIsNull(grivation))
615 georef_element.writeAttribute(literal::grivation, grivation, declinationPrecision());
616
617 if (!qIsNull(map_ref_point.lengthSquared()))
618 {
619 XmlElementWriter ref_point_element(xml, literal::ref_point);
620 ref_point_element.writeAttribute(literal::x, map_ref_point.x());
621 ref_point_element.writeAttribute(literal::y, map_ref_point.y());
622 }
623
624 {
625 XmlElementWriter crs_element(xml, literal::projected_crs);
626 if (!projected_crs_id.isEmpty())
627 crs_element.writeAttribute(literal::id, projected_crs_id);
628
629 if (!projected_crs_spec.isEmpty())
630 {
631 XmlElementWriter spec_element(xml, literal::spec);
632 spec_element.writeAttribute(literal::language, literal::proj_4);
633 xml.writeCharacters(projected_crs_spec);
634 }
635
636 if (!projected_crs_parameters.empty())
637 {
638 for (const auto& projected_crs_parameter : projected_crs_parameters)
639 {
640 XmlElementWriter parameter_element(xml, literal::parameter);
641 xml.writeCharacters(projected_crs_parameter);
642 Q_UNUSED(parameter_element); // Suppress compiler warnings
643 }
644 }
645
646 if (!qIsNull(projected_ref_point.manhattanLength()))
647 {
648 XmlElementWriter ref_point_element(xml, literal::ref_point);
649 ref_point_element.writeAttribute(literal::x, projected_ref_point.x(), 6);
650 ref_point_element.writeAttribute(literal::y, projected_ref_point.y(), 6);
651 }
652 }
653
654 if (state == Normal)
655 {
656 XmlElementWriter crs_element(xml, literal::geographic_crs);
657 crs_element.writeAttribute(literal::id, literal::geographic_coordinates);
658 {
659 XmlElementWriter spec_element(xml, literal::spec);
660 spec_element.writeAttribute(literal::language, literal::proj_4);
661 xml.writeCharacters(geographic_crs_spec);
662 }
663 if (XMLFileFormat::active_version < 6)
664 {
665 // Legacy compatibility
666 XmlElementWriter ref_point_element(xml, literal::ref_point);
667 ref_point_element.writeAttribute(literal::lat, degToRad(geographic_ref_point.latitude()), 10);
668 ref_point_element.writeAttribute(literal::lon, degToRad(geographic_ref_point.longitude()), 10);
669 }
670 {
671 XmlElementWriter ref_point_element(xml, literal::ref_point_deg);
672 ref_point_element.writeAttribute(literal::lat, geographic_ref_point.latitude(), 8);
673 ref_point_element.writeAttribute(literal::lon, geographic_ref_point.longitude(), 8);
674 }
675 }
676 }
677
678
setState(Georeferencing::State value)679 void Georeferencing::setState(Georeferencing::State value)
680 {
681 if (state != value)
682 {
683 state = value;
684 updateTransformation();
685
686 if (state != Normal)
687 {
688 setProjectedCRS(QStringLiteral("Local"), {});
689 updateGridCompensation();
690 }
691
692 emit stateChanged();
693 }
694 }
695
setScaleDenominator(int value)696 void Georeferencing::setScaleDenominator(int value)
697 {
698 Q_ASSERT(value > 0);
699 if (scale_denominator != (unsigned int)value)
700 {
701 scale_denominator = value;
702 updateTransformation();
703 }
704 }
705
setCombinedScaleFactor(double value)706 void Georeferencing::setCombinedScaleFactor(double value)
707 {
708 Q_ASSERT(value > 0);
709 double combined_scale_factor = roundScaleFactor(value);
710 double auxiliary_scale_factor = roundScaleFactor(value / grid_scale_factor);
711 setScaleFactors(combined_scale_factor, auxiliary_scale_factor);
712 }
713
setAuxiliaryScaleFactor(double value)714 void Georeferencing::setAuxiliaryScaleFactor(double value)
715 {
716 Q_ASSERT(value > 0);
717 double auxiliary_scale_factor = roundScaleFactor(value);
718 double combined_scale_factor = roundScaleFactor(value * grid_scale_factor);
719 setScaleFactors(combined_scale_factor, auxiliary_scale_factor);
720 }
721
setScaleFactors(double combined_scale_factor,double auxiliary_scale_factor)722 void Georeferencing::setScaleFactors(double combined_scale_factor, double auxiliary_scale_factor)
723 {
724 bool combined_change = combined_scale_factor != this->combined_scale_factor;
725 bool auxiliary_change = auxiliary_scale_factor != this->auxiliary_scale_factor;
726 if (combined_change || auxiliary_change)
727 {
728 this->auxiliary_scale_factor = auxiliary_scale_factor;
729 this->combined_scale_factor = combined_scale_factor;
730 if (combined_change)
731 updateTransformation();
732 if (auxiliary_change)
733 emit auxiliaryScaleFactorChanged();
734 }
735 }
736
setDeclination(double value)737 void Georeferencing::setDeclination(double value)
738 {
739 double declination = roundDeclination(value);
740 double grivation = roundDeclination(value - convergence);
741 setDeclinationAndGrivation(declination, grivation);
742 }
743
setGrivation(double value)744 void Georeferencing::setGrivation(double value)
745 {
746 double grivation = roundDeclination(value);
747 double declination = roundDeclination(value + convergence);
748 setDeclinationAndGrivation(declination, grivation);
749 }
750
setDeclinationAndGrivation(double declination,double grivation)751 void Georeferencing::setDeclinationAndGrivation(double declination, double grivation)
752 {
753 bool declination_change = declination != this->declination;
754 bool grivation_change = grivation != this->grivation;
755 if (declination_change || grivation_change)
756 {
757 this->declination = declination;
758 this->grivation = grivation;
759 this->grivation_error = 0.0;
760 if (grivation_change)
761 updateTransformation();
762
763 if (declination_change)
764 emit declinationChanged();
765 }
766 }
767
setMapRefPoint(const MapCoord & point)768 void Georeferencing::setMapRefPoint(const MapCoord& point)
769 {
770 if (map_ref_point != point)
771 {
772 map_ref_point = point;
773 updateTransformation();
774 }
775 }
776
setProjectedRefPoint(const QPointF & point,bool update_grivation,bool update_scale_factor)777 void Georeferencing::setProjectedRefPoint(const QPointF& point, bool update_grivation, bool update_scale_factor)
778 {
779 if (projected_ref_point != point || state == Normal)
780 {
781 projected_ref_point = point;
782 bool ok = {};
783 LatLon new_geo_ref_point;
784
785 switch (state)
786 {
787 default:
788 qWarning("Undefined georeferencing state");
789 Q_FALLTHROUGH();
790 case Local:
791 break;
792 case Normal:
793 new_geo_ref_point = toGeographicCoords(point, &ok);
794 if (ok && new_geo_ref_point != geographic_ref_point)
795 {
796 geographic_ref_point = new_geo_ref_point;
797 updateGridCompensation();
798 if (update_grivation)
799 updateGrivation();
800 if (update_scale_factor)
801 updateCombinedScaleFactor();
802 emit projectionChanged();
803 }
804 }
805 updateTransformation();
806 }
807 }
808
getProjectedCRSName() const809 QString Georeferencing::getProjectedCRSName() const
810 {
811 QString name = tr("Local");
812 if (auto temp = CRSTemplateRegistry().find(projected_crs_id))
813 name = temp->name();
814
815 return name;
816 }
817
getProjectedCoordinatesName() const818 QString Georeferencing::getProjectedCoordinatesName() const
819 {
820 QString name = tr("Local coordinates");
821 if (auto temp = CRSTemplateRegistry().find(projected_crs_id))
822 name = temp->coordinatesName(projected_crs_parameters);
823
824 return name;
825 }
826
getConvergence() const827 double Georeferencing::getConvergence() const
828 {
829 return convergence;
830 }
831
updateGridCompensation()832 void Georeferencing::updateGridCompensation()
833 {
834 convergence = 0.0;
835 grid_scale_factor = 1.0;
836
837 if (state != Normal || !isValid())
838 return;
839
840 const double delta = 1000.0; // meters
841
842 QString local_crs_spec = QString::fromLatin1("+proj=sterea +lat_0=%1 +lon_0=%2 +ellps=WGS84 +units=m")
843 .arg(geographic_ref_point.latitude(), 0, 'f')
844 .arg(geographic_ref_point.longitude(), 0, 'f');
845 ProjTransform local_proj_transform(local_crs_spec);
846 if (!local_proj_transform.isValid())
847 return;
848
849 // Determine 1 km baselines west-east and south-north on the ellipsoid.
850 bool ok_east_point = false;
851 const QPointF east_projected_coords {delta/2, 0};
852 const LatLon east_point = local_proj_transform.inverse(east_projected_coords, &ok_east_point);
853
854 bool ok_north_point = false;
855 const QPointF north_projected_coords {0, delta/2};
856 const LatLon north_point = local_proj_transform.inverse(north_projected_coords, &ok_north_point);
857
858 bool ok_west_point = false;
859 const QPointF west_projected_coords {-delta/2, 0};
860 const LatLon west_point = local_proj_transform.inverse(west_projected_coords, &ok_west_point);
861
862 bool ok_south_point = false;
863 const QPointF south_projected_coords {0, -delta/2};
864 const LatLon south_point = local_proj_transform.inverse(south_projected_coords, &ok_south_point);
865
866 if (!(ok_east_point && ok_north_point && ok_west_point && ok_south_point))
867 return;
868
869 // Get projected coordinates on same meridian and on same parallel around reference point.
870 bool ok_east = false;
871 QPointF projected_east = toProjectedCoords(east_point, &ok_east);
872 bool ok_north = false;
873 QPointF projected_north = toProjectedCoords(north_point, &ok_north);
874 bool ok_west = false;
875 QPointF projected_west = toProjectedCoords(west_point, &ok_west);
876 bool ok_south = false;
877 QPointF projected_south = toProjectedCoords(south_point, &ok_south);
878 if (!(ok_east && ok_north && ok_west && ok_south))
879 return;
880
881 // Points on the same meridian
882 const double d_northing_dy = (projected_north.y() - projected_south.y()) / delta;
883 const double d_easting_dy = (projected_north.x() - projected_south.x()) / delta;
884 // Points on the same parallel
885 const double d_northing_dx = (projected_east.y() - projected_west.y()) / delta;
886 const double d_easting_dx = (projected_east.x() - projected_west.x()) / delta;
887
888 // A transform with a tiny (or negative) determinant is nonsense for a map, and
889 // would cause blow-ups.
890 const double determinant = d_easting_dx*d_northing_dy - d_northing_dx*d_easting_dy;
891 if (determinant < 0.00000000001)
892 return;
893
894 // This is the angle between true azimuth and grid azimuth.
895 // In case of deformation, the convergence varies with direction and this is an average.
896 convergence = qRadiansToDegrees(atan2(d_northing_dx - d_easting_dy,
897 d_easting_dx + d_northing_dy));
898
899 // This is the scale factor from true distance to grid distance.
900 // In case of deformation, the scale factor varies with direction and this is an average.
901 grid_scale_factor = sqrt(determinant);
902 }
903
setGeographicRefPoint(LatLon lat_lon,bool update_grivation,bool update_scale_factor)904 void Georeferencing::setGeographicRefPoint(LatLon lat_lon, bool update_grivation, bool update_scale_factor)
905 {
906 bool geo_ref_point_changed = geographic_ref_point != lat_lon;
907 if (geo_ref_point_changed || state == Normal)
908 {
909 geographic_ref_point = lat_lon;
910 if (state != Normal)
911 setState(Normal);
912
913 bool ok = {};
914 QPointF new_projected_ref = toProjectedCoords(lat_lon, &ok);
915 if (ok && new_projected_ref != projected_ref_point)
916 {
917 projected_ref_point = new_projected_ref;
918 updateGridCompensation();
919 if (update_grivation)
920 updateGrivation();
921 if (update_scale_factor)
922 updateCombinedScaleFactor();
923 updateTransformation();
924 emit projectionChanged();
925 }
926 else if (geo_ref_point_changed)
927 {
928 emit projectionChanged();
929 }
930 }
931 }
932
updateTransformation()933 void Georeferencing::updateTransformation()
934 {
935 // Use the grivation and combined scale factor.
936 QTransform transform;
937 transform.translate(projected_ref_point.x(), projected_ref_point.y());
938 transform.rotate(-grivation);
939 double scale = combined_scale_factor * scale_denominator / 1000.0; // to meters
940 transform.scale(scale, -scale);
941 transform.translate(-map_ref_point.x(), -map_ref_point.y());
942
943 if (to_projected != transform)
944 {
945 to_projected = transform;
946 from_projected = transform.inverted();
947 emit transformationChanged();
948 }
949 }
950
updateCombinedScaleFactor()951 void Georeferencing::updateCombinedScaleFactor()
952 {
953 setAuxiliaryScaleFactor(auxiliary_scale_factor);
954 }
955
initAuxiliaryScaleFactor()956 void Georeferencing::initAuxiliaryScaleFactor()
957 {
958 setCombinedScaleFactor(combined_scale_factor);
959 }
960
updateGrivation()961 void Georeferencing::updateGrivation()
962 {
963 setDeclination(declination);
964 }
965
initDeclination()966 void Georeferencing::initDeclination()
967 {
968 setGrivation(grivation);
969 }
970
setTransformationDirectly(const QTransform & transform)971 void Georeferencing::setTransformationDirectly(const QTransform& transform)
972 {
973 if (transform != to_projected)
974 {
975 to_projected = transform;
976 from_projected = to_projected.inverted();
977 emit transformationChanged();
978 }
979 }
980
mapToProjected() const981 const QTransform& Georeferencing::mapToProjected() const
982 {
983 return to_projected;
984 }
985
projectedToMap() const986 const QTransform& Georeferencing::projectedToMap() const
987 {
988 return from_projected;
989 }
990
991
992
setProjectedCRS(const QString & id,QString spec,std::vector<QString> params)993 bool Georeferencing::setProjectedCRS(const QString& id, QString spec, std::vector<QString> params)
994 {
995 // Default return value if no change is necessary
996 bool ok = (state == Normal || projected_crs_spec.isEmpty());
997
998 for (const auto& substitution : spec_substitutions)
999 {
1000 if (QLatin1String(substitution[0]) == spec)
1001 {
1002 spec = QString::fromLatin1(substitution[1]);
1003 break;
1004 }
1005 }
1006
1007 // Changes in params shall already be recorded in spec
1008 if (projected_crs_id != id
1009 || projected_crs_spec != spec
1010 || projected_crs_parameters.size() != params.size()
1011 || !std::equal(begin(params), end(params), begin(projected_crs_parameters))
1012 || (!ok && !spec.isEmpty()) )
1013 {
1014 projected_crs_id = id;
1015 projected_crs_spec = spec;
1016 if (projected_crs_spec.isEmpty())
1017 {
1018 projected_crs_parameters.clear();
1019 proj_transform = {};
1020 ok = (state != Normal);
1021 }
1022 else
1023 {
1024 projected_crs_parameters.swap(params); // params was passed by value!
1025 proj_transform = {projected_crs_spec};
1026 ok = proj_transform.isValid();
1027 if (ok && state != Normal)
1028 setState(Normal);
1029 }
1030 if (!ok)
1031 updateGridCompensation();
1032
1033 emit projectionChanged();
1034 }
1035
1036 return ok;
1037 }
1038
toProjectedCoords(const MapCoord & map_coords) const1039 QPointF Georeferencing::toProjectedCoords(const MapCoord& map_coords) const
1040 {
1041 return to_projected.map(QPointF(map_coords));
1042 }
1043
toProjectedCoords(const MapCoordF & map_coords) const1044 QPointF Georeferencing::toProjectedCoords(const MapCoordF& map_coords) const
1045 {
1046 return to_projected.map(map_coords);
1047 }
1048
toMapCoords(const QPointF & projected_coords) const1049 MapCoord Georeferencing::toMapCoords(const QPointF& projected_coords) const
1050 {
1051 return MapCoord(from_projected.map(projected_coords));
1052 }
1053
toMapCoordF(const QPointF & projected_coords) const1054 MapCoordF Georeferencing::toMapCoordF(const QPointF& projected_coords) const
1055 {
1056 return MapCoordF(from_projected.map(projected_coords));
1057 }
1058
toGeographicCoords(const MapCoordF & map_coords,bool * ok) const1059 LatLon Georeferencing::toGeographicCoords(const MapCoordF& map_coords, bool* ok) const
1060 {
1061 return toGeographicCoords(toProjectedCoords(map_coords), ok);
1062 }
1063
toGeographicCoords(const QPointF & projected_coords,bool * ok) const1064 LatLon Georeferencing::toGeographicCoords(const QPointF& projected_coords, bool* ok) const
1065 {
1066 return proj_transform.isValid() ? proj_transform.inverse(projected_coords, ok) : LatLon{};
1067 }
1068
toProjectedCoords(const LatLon & lat_lon,bool * ok) const1069 QPointF Georeferencing::toProjectedCoords(const LatLon& lat_lon, bool* ok) const
1070 {
1071 return proj_transform.isValid() ? proj_transform.forward(lat_lon, ok) : QPointF{};
1072 }
1073
toMapCoords(const LatLon & lat_lon,bool * ok) const1074 MapCoord Georeferencing::toMapCoords(const LatLon& lat_lon, bool* ok) const
1075 {
1076 return toMapCoords(toProjectedCoords(lat_lon, ok));
1077 }
1078
toMapCoordF(const LatLon & lat_lon,bool * ok) const1079 MapCoordF Georeferencing::toMapCoordF(const LatLon& lat_lon, bool* ok) const
1080 {
1081 return toMapCoordF(toProjectedCoords(lat_lon, ok));
1082 }
1083
toMapCoordF(const Georeferencing * other,const MapCoordF & map_coords,bool * ok) const1084 MapCoordF Georeferencing::toMapCoordF(const Georeferencing* other, const MapCoordF& map_coords, bool* ok) const
1085 {
1086 if (!other)
1087 {
1088 if (ok)
1089 *ok = true;
1090 return map_coords;
1091 }
1092
1093 if (isLocal() || other->isLocal())
1094 {
1095 if (ok)
1096 *ok = true;
1097 return toMapCoordF(other->toProjectedCoords(map_coords));
1098 }
1099
1100 bool ok_forward = proj_transform.isValid();
1101 bool ok_inverse = other->proj_transform.isValid();
1102 QPointF projected = other->toProjectedCoords(map_coords);
1103 if (ok_forward && ok_inverse) {
1104 // Use geographic coordinates as intermediate step to enforce
1105 // that coordinates are assumed to have WGS84 datum if datum is specified in only one CRS spec.
1106 /// \todo Use direct pipeline instead of intermediate WGS84
1107 bool ok_inverse, ok_forward;
1108 projected = proj_transform.forward(other->proj_transform.inverse(projected, &ok_inverse), &ok_forward);
1109 }
1110 if (ok)
1111 *ok = ok_inverse && ok_forward;
1112 return toMapCoordF(projected);
1113 }
1114
getErrorText() const1115 QString Georeferencing::getErrorText() const
1116 {
1117 return proj_transform.errorText();
1118 }
1119
radToDeg(double val)1120 double Georeferencing::radToDeg(double val)
1121 {
1122 return qRadiansToDegrees(val);
1123 }
1124
degToRad(double val)1125 double Georeferencing::degToRad(double val)
1126 {
1127 return qDegreesToRadians(val);
1128 }
1129
degToDMS(double val)1130 QString Georeferencing::degToDMS(double val)
1131 {
1132 qint64 tmp = val * 360000;
1133 int csec = tmp % 6000;
1134 tmp = tmp / 6000;
1135 int min = tmp % 60;
1136 int deg = tmp / 60;
1137 QString ret = QString::fromUtf8("%1°%2'%3\"").arg(deg).arg(min).arg(QLocale().toString(csec/100.0,'f',2));
1138 return ret;
1139 }
1140
operator <<(QDebug dbg,const Georeferencing & georef)1141 QDebug operator<<(QDebug dbg, const Georeferencing &georef)
1142 {
1143 dbg.nospace()
1144 << "Georeferencing(1:" << georef.scale_denominator
1145 << " " << georef.combined_scale_factor
1146 << " " << georef.declination
1147 << " " << georef.grivation
1148 << "deg, " << georef.projected_crs_id
1149 << " (" << georef.projected_crs_spec
1150 << ") " << QString::number(georef.projected_ref_point.x(), 'f', 8) << "," << QString::number(georef.projected_ref_point.y(), 'f', 8);
1151 if (georef.isLocal())
1152 dbg.nospace() << ", local)";
1153 else
1154 dbg.nospace() << ", geographic)";
1155
1156 return dbg.space();
1157 }
1158
1159
1160 } // namespace OpenOrienteering
1161