1
/******************************************************************************
2
* $Id: lasspatialreference.cpp 1104 2009-03-17 04:15:19Z hobu $
4
* Project: libLAS - http://liblas.org - A BSD library for LAS format data.
5
* Purpose: LAS Spatial Reference class
6
* Author: Howard Butler, hobu.inc@gmail.com
7
* Frank Warmerdam, warmerdam@pobox.com
9
******************************************************************************
10
* Copyright (c) 2009, Howard Butler
12
* All rights reserved.
14
* Redistribution and use in source and binary forms, with or without
15
* modification, are permitted provided that the following
18
* * Redistributions of source code must retain the above copyright
19
* notice, this list of conditions and the following disclaimer.
20
* * Redistributions in binary form must reproduce the above copyright
21
* notice, this list of conditions and the following disclaimer in
22
* the documentation and/or other materials provided
23
* with the distribution.
24
* * Neither the name of the Martin Isenburg or Iowa Department
25
* of Natural Resources nor the names of its contributors may be
26
* used to endorse or promote products derived from this software
27
* without specific prior written permission.
29
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
30
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
31
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
32
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
33
* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
34
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
35
* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
36
* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
37
* AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
38
* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
39
* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
41
****************************************************************************/
47
// Supress inclusion of cpl_serv.h per #194, perhaps remove one day
48
// when libgeotiff 1.4.0+ is widely used
49
#define CPL_SERV_H_INCLUDED
51
#include <ogr_srs_api.h>
53
#include <geo_normalize.h>
54
#include <geovalues.h>
55
#include <ogr_spatialref.h>
58
#include <cpl_multiproc.h>
62
#ifdef HAVE_LIBGEOTIFF
64
#include <geo_simpletags.h>
65
#include <geo_normalize.h>
66
#include <geo_simpletags.h>
67
#include <geovalues.h>
68
#endif // HAVE_LIBGEOTIFF
70
#include <liblas/spatialreference.hpp>
71
#include <liblas/detail/private_utility.hpp>
72
#include <liblas/external/property_tree/xml_parser.hpp>
74
#include <boost/concept_check.hpp>
75
#include <boost/cstdint.hpp>
82
# include "cpl_conv.h"
85
using namespace boost;
89
SpatialReference::SpatialReference()
97
SpatialReference::SpatialReference(std::vector<VariableRecord> const& vlrs)
105
SpatialReference::SpatialReference(SpatialReference const& other)
110
SetVLRs(other.GetVLRs());
114
SpatialReference& SpatialReference::operator=(SpatialReference const& rhs)
118
SetVLRs(rhs.GetVLRs());
125
bool SpatialReference::operator==(const SpatialReference& input) const
129
OGRSpatialReferenceH current = OSRNewSpatialReference(GetWKT(eCompoundOK, false).c_str());
130
OGRSpatialReferenceH other = OSRNewSpatialReference(input.GetWKT(eCompoundOK, false).c_str());
132
int output = OSRIsSame(current, other);
134
OSRDestroySpatialReference( current );
135
OSRDestroySpatialReference( other );
137
return bool(output == 1);
140
boost::ignore_unused_variable_warning(input);
141
throw std::runtime_error ("SpatialReference equality testing not available without GDAL+libgeotiff support");
146
SpatialReference::~SpatialReference()
148
#ifdef HAVE_LIBGEOTIFF
162
/// Keep a copy of the VLRs that are related to GeoTIFF SRS information.
163
void SpatialReference::SetVLRs(std::vector<VariableRecord> const& vlrs)
166
std::string const uid("LASF_Projection");
168
// Wipe out any existing VLRs that might exist on the SpatialReference
171
// We only copy VLR records from the list which are related to GeoTIFF keys.
172
// They must have an id of "LASF_Projection" and a record id that's related.
173
std::vector<VariableRecord>::const_iterator it;
174
for (it = vlrs.begin(); it != vlrs.end(); ++it)
176
VariableRecord const& vlr = *it;
179
m_vlrs.push_back(vlr);
184
void SpatialReference::AddVLR(VariableRecord const& vlr)
188
m_vlrs.push_back(vlr);
192
bool SpatialReference::IsGeoVLR(VariableRecord const& vlr) const
194
std::string const las_projid("LASF_Projection");
195
std::string const liblas_id("liblas");
197
// GTIFF_GEOKEYDIRECTORY == 34735
198
if (las_projid == vlr.GetUserId(true).c_str() && 34735 == vlr.GetRecordId())
203
// GTIFF_DOUBLEPARAMS == 34736
204
if (las_projid == vlr.GetUserId(true).c_str() && 34736 == vlr.GetRecordId())
209
// GTIFF_ASCIIPARAMS == 34737
210
if (las_projid == vlr.GetUserId(true).c_str() && 34737 == vlr.GetRecordId())
216
if (liblas_id == vlr.GetUserId(true).c_str() && 2112 == vlr.GetRecordId())
224
std::vector<VariableRecord> SpatialReference::GetVLRs() const
229
void SpatialReference::ClearVLRs( GeoVLRType eType )
232
std::vector<VariableRecord>::iterator it;
233
std::string const liblas_id("liblas");
235
for (it = m_vlrs.begin(); it != m_vlrs.end(); )
237
VariableRecord const& vlr = *it;
240
// for now we can assume all m_vlrs records are LASF_Projection.
241
if( eType == eOGRWKT &&
242
2112 == vlr.GetRecordId() &&
243
liblas_id == vlr.GetUserId(true).c_str() )
246
else if( eType == eGeoTIFF
247
&& (34735 == vlr.GetRecordId()
248
|| 34736 == vlr.GetRecordId()
249
|| 34737 == vlr.GetRecordId()) )
253
it = m_vlrs.erase( it );
258
if( eType == eOGRWKT )
260
else if( eType == eGeoTIFF )
262
#ifdef HAVE_LIBGEOTIFF
277
void SpatialReference::ResetVLRs()
282
#ifdef HAVE_LIBGEOTIFF
299
throw std::invalid_argument("m_tiff was null, cannot reset VLRs without m_tiff");
302
throw std::invalid_argument("m_gtiff was null, cannot reset VLRs without m_gtiff");
304
//GTIFF_GEOKEYDIRECTORY == 34735
305
ret = ST_GetKey(m_tiff, 34735, &kcount, &ktype, (void**)&kdata);
308
VariableRecord record;
310
record.SetRecordId(34735);
311
record.SetUserId("LASF_Projection");
312
record.SetDescription("GeoTIFF GeoKeyDirectoryTag");
313
std::vector<uint8_t> data;
315
// Shorts are 2 bytes in length
316
uint16_t length = 2 * static_cast<uint16_t>(kcount);
317
record.SetRecordLength(length);
319
// Copy the data into the data vector
320
for (i = 0; i < kcount; i++)
324
uint8_t* v = reinterpret_cast<uint8_t*>(&kvalue);
326
data.push_back(v[0]);
327
data.push_back(v[1]);
330
record.SetData(data);
331
m_vlrs.push_back(record);
334
// GTIFF_DOUBLEPARAMS == 34736
335
ret = ST_GetKey(m_tiff, 34736, &dcount, &dtype, (void**)&ddata);
338
VariableRecord record;
340
record.SetRecordId(34736);
341
record.SetUserId("LASF_Projection");
342
record.SetDescription("GeoTIFF GeoDoubleParamsTag");
344
std::vector<uint8_t> data;
346
// Doubles are 8 bytes in length
347
uint16_t length = 8 * static_cast<uint16_t>(dcount);
348
record.SetRecordLength(length);
350
// Copy the data into the data vector
351
for (i=0; i<dcount;i++)
355
uint8_t* v = reinterpret_cast<uint8_t*>(&dvalue);
357
data.push_back(v[0]);
358
data.push_back(v[1]);
359
data.push_back(v[2]);
360
data.push_back(v[3]);
361
data.push_back(v[4]);
362
data.push_back(v[5]);
363
data.push_back(v[6]);
364
data.push_back(v[7]);
366
record.SetData(data);
367
m_vlrs.push_back(record);
370
// GTIFF_ASCIIPARAMS == 34737
371
ret = ST_GetKey(m_tiff, 34737, &acount, &atype, (void**)&adata);
374
VariableRecord record;
376
record.SetRecordId(34737);
377
record.SetUserId("LASF_Projection");
378
record.SetDescription("GeoTIFF GeoAsciiParamsTag");
379
std::vector<uint8_t> data;
381
// whoa. If the returned count was 0, it is because there
382
// was a bug in libgeotiff that existed before r1531 where it
383
// didn't calculate the string length for an ASCII geotiff tag.
384
// We need to throw an exception in this case because we're
385
// screwed, and if we don't we'll end up writing bad GeoTIFF keys.
388
throw std::runtime_error("GeoTIFF ASCII key with no returned size. "
389
"Upgrade your libgeotiff to a version greater "
390
"than r1531 (libgeotiff 1.2.6)");
393
uint16_t length = static_cast<uint16_t>(acount);
394
record.SetRecordLength(length);
396
// Copy the data into the data vector
397
for (i=0; i<acount;i++)
400
uint8_t* v = reinterpret_cast<uint8_t*>(&avalue);
401
data.push_back(v[0]);
403
record.SetData(data);
405
if (data.size() > (std::numeric_limits<boost::uint16_t>::max()))
407
std::ostringstream oss;
408
std::vector<uint8_t>::size_type overrun = data.size() - static_cast<std::vector<uint8_t>::size_type>(std::numeric_limits<boost::uint16_t>::max());
409
oss << "The size of the GeoTIFF GeoAsciiParamsTag, " << data.size() << ", is " << overrun
410
<< " bytes too large to fit inside the maximum size of a VLR which is "
411
<< (std::numeric_limits<boost::uint16_t>::max()) << " bytes.";
412
throw std::runtime_error(oss.str());
416
m_vlrs.push_back(record);
418
#endif // ndef HAVE_LIBGEOTIFF
422
m_wkt = GetWKT( eCompoundOK );
424
// Add a WKT VLR if we have a WKT definition.
427
VariableRecord wkt_record;
428
std::vector<uint8_t> data;
429
const uint8_t* wkt_bytes = reinterpret_cast<const uint8_t*>(m_wkt.c_str());
431
wkt_record.SetRecordId( 2112 );
432
wkt_record.SetUserId("liblas");
433
wkt_record.SetDescription( "OGR variant of OpenGIS WKT SRS" );
435
// Would you be surprised if this remarshalling of bytes
436
// was annoying to me? FrankW
437
while( *wkt_bytes != 0 )
438
data.push_back( *(wkt_bytes++) );
440
data.push_back( '\0' );
442
if (data.size() > std::numeric_limits<boost::uint16_t>::max())
444
std::ostringstream oss;
445
std::vector<uint8_t>::size_type overrun = data.size() - static_cast<std::vector<uint8_t>::size_type>(std::numeric_limits<boost::uint16_t>::max());
446
oss << "The size of the wkt, " << data.size() << ", is " << overrun
447
<< " bytes too large to fit inside the maximum size of a VLR which is "
448
<< std::numeric_limits<boost::uint16_t>::max() << " bytes.";
449
throw std::runtime_error(oss.str());
454
wkt_record.SetRecordLength( static_cast<boost::uint16_t>(data.size()) );
455
wkt_record.SetData(data);
457
// not to speak of this additional copy!
458
m_vlrs.push_back( wkt_record );
462
void SpatialReference::SetGTIF(GTIF* pgtiff, ST_TIFF* ptiff)
464
m_gtiff = (GTIF*)pgtiff;
465
m_tiff = (ST_TIFF*)ptiff;
470
const GTIF* SpatialReference::GetGTIF()
472
#ifndef HAVE_LIBGEOTIFF
476
// If we already have m_gtiff and m_tiff, that is because we have
477
// already called GetGTIF once before. VLRs ultimately drive how the
478
// SpatialReference is defined, not the GeoTIFF keys.
491
m_tiff = ST_Create();
492
std::string const uid("LASF_Projection");
494
// Nothing is going to happen here if we don't have any VLRs describing
495
// SRS information on the SpatialReference.
496
for (uint16_t i = 0; i < m_vlrs.size(); ++i)
498
VariableRecord record = m_vlrs[i];
499
std::vector<uint8_t> data = record.GetData();
500
if (uid == record.GetUserId(true).c_str() && 34735 == record.GetRecordId())
502
int count = data.size()/sizeof(int16_t);
503
short *data_s = (short *) &(data[0]);
505
// discard invalid "zero" geotags some software emits.
507
&& data_s[count-1] == 0
508
&& data_s[count-2] == 0
509
&& data_s[count-3] == 0
510
&& data_s[count-4] == 0 )
516
ST_SetKey(m_tiff, record.GetRecordId(), count, STT_SHORT, data_s);
519
if (uid == record.GetUserId(true).c_str() && 34736 == record.GetRecordId())
521
int count = data.size() / sizeof(double);
522
ST_SetKey(m_tiff, record.GetRecordId(), count, STT_DOUBLE, &(data[0]));
525
if (uid == record.GetUserId(true).c_str() && 34737 == record.GetRecordId())
527
int count = data.size()/sizeof(uint8_t);
528
ST_SetKey(m_tiff, record.GetRecordId(), count, STT_ASCII, &(data[0]));
532
m_gtiff = GTIFNewSimpleTags(m_tiff);
534
throw std::runtime_error("The geotiff keys could not be read from VLR records");
540
std::string SpatialReference::GetWKT( WKTModeFlag mode_flag) const
542
return GetWKT(mode_flag, false);
545
/// Fetch the SRS as WKT
546
std::string SpatialReference::GetWKT(WKTModeFlag mode_flag , bool pretty) const
549
boost::ignore_unused_variable_warning(mode_flag);
550
boost::ignore_unused_variable_warning(pretty);
552
// we don't have a way of making this pretty, or of stripping the compound wrapper.
556
// If we already have Well Known Text then try return it, possibly
557
// after some preprocessing.
560
std::string result_wkt = m_wkt;
562
if( (mode_flag == eHorizontalOnly
563
&& strstr(result_wkt.c_str(),"COMPD_CS") != NULL)
566
OGRSpatialReference* poSRS = (OGRSpatialReference*) OSRNewSpatialReference(result_wkt.c_str());
569
if( mode_flag == eHorizontalOnly )
570
poSRS->StripVertical();
573
poSRS->exportToPrettyWkt(&pszWKT, FALSE );
575
poSRS->exportToWkt( &pszWKT );
577
OSRDestroySpatialReference( poSRS );
586
// Otherwise build WKT from GeoTIFF VLRs.
592
return std::string();
595
if (GTIFGetDefn(m_gtiff, &sGTIFDefn))
597
pszWKT = GTIFGetOGISDefn( m_gtiff, &sGTIFDefn );
600
OGRSpatialReference* poSRS = (OGRSpatialReference*) OSRNewSpatialReference(NULL);
601
char *pszOrigWKT = pszWKT;
602
poSRS->importFromWkt( &pszOrigWKT );
606
poSRS->exportToPrettyWkt(&pszWKT, false);
607
OSRDestroySpatialReference( poSRS );
611
// save this for future calls, etc.
612
// m_wkt = std::string( pszWKT );
614
// Older versions of GDAL lack StripVertical(), but should never
615
// actually return COMPD_CS anyways.
616
#if (GDAL_VERSION_NUM >= 1700) && (GDAL_RELEASE_DATE >= 20100110)
618
&& mode_flag == eHorizontalOnly
619
&& strstr(pszWKT,"COMPD_CS") != NULL )
621
OGRSpatialReference* poSRS = (OGRSpatialReference*) OSRNewSpatialReference(NULL);
622
char *pszOrigWKT = pszWKT;
623
poSRS->importFromWkt( &pszOrigWKT );
628
poSRS->StripVertical();
630
poSRS->exportToPrettyWkt(&pszWKT, false);
632
poSRS->exportToWkt( &pszWKT );
634
OSRDestroySpatialReference( poSRS );
637
boost::ignore_unused_variable_warning(mode_flag);
642
std::string tmp(pszWKT);
647
return std::string();
651
void SpatialReference::SetFromUserInput(std::string const& v)
656
const char* input = v.c_str();
658
// OGRSpatialReference* poSRS = (OGRSpatialReference*) OSRNewSpatialReference(NULL);
659
OGRSpatialReference srs(NULL);
660
if (OGRERR_NONE != srs.SetFromUserInput(const_cast<char *> (input)))
662
throw std::invalid_argument("could not import coordinate system into OSRSpatialReference SetFromUserInput");
665
srs.exportToWkt(&poWKT);
667
std::string tmp(poWKT);
672
boost::ignore_unused_variable_warning(v);
673
throw std::runtime_error("GDAL is not available, SpatialReference could not be set from WKT");
677
void SpatialReference::SetWKT(std::string const& v)
688
ret = GTIFSetFromOGISDefn( m_gtiff, v.c_str() );
691
throw std::invalid_argument("could not set m_gtiff from WKT");
694
ret = GTIFWriteKeys(m_gtiff);
697
throw std::runtime_error("The geotiff keys could not be written");
700
boost::ignore_unused_variable_warning(v);
706
void SpatialReference::SetVerticalCS(boost::int32_t verticalCSType,
707
std::string const& citation,
708
boost::int32_t verticalDatum,
709
boost::int32_t verticalUnits)
716
#ifdef HAVE_LIBGEOTIFF
717
if( verticalCSType != KvUserDefined && verticalCSType > 0 )
718
GTIFKeySet( m_gtiff, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
722
GTIFKeySet( m_gtiff, VerticalCitationGeoKey, TYPE_ASCII, 0,
725
if( verticalDatum > 0 && verticalDatum != KvUserDefined )
726
GTIFKeySet( m_gtiff, VerticalDatumGeoKey, TYPE_SHORT, 1,
729
if( verticalUnits > 0 && verticalUnits != KvUserDefined )
730
GTIFKeySet( m_gtiff, VerticalUnitsGeoKey, TYPE_SHORT, 1,
733
int ret = GTIFWriteKeys(m_gtiff);
736
throw std::runtime_error("The geotiff keys could not be written");
739
// Clear WKT so it gets regenerated
740
m_wkt = std::string("");
744
boost::ignore_unused_variable_warning(citation);
745
boost::ignore_unused_variable_warning(verticalUnits);
746
boost::ignore_unused_variable_warning(verticalDatum);
747
boost::ignore_unused_variable_warning(verticalCSType);
748
#endif /* def HAVE_LIBGEOTIFF */
751
std::string SpatialReference::GetProj4() const
755
std::string wkt = GetWKT(eCompoundOK);
756
const char* poWKT = wkt.c_str();
758
OGRSpatialReference srs(NULL);
759
if (OGRERR_NONE != srs.importFromWkt(const_cast<char **> (&poWKT)))
761
return std::string();
765
srs.exportToProj4(&proj4);
766
std::string tmp(proj4);
772
// if we have libgeotiff but not GDAL, we'll use the
773
// simple method in libgeotiff
774
#if defined(HAVE_LIBGEOTIFF) && !defined(HAVE_GDAL)
778
if (m_gtiff && GTIFGetDefn(m_gtiff, &defn))
780
char* proj4def = GTIFGetProj4Defn(&defn);
781
std::string tmp(proj4def);
782
GTIFFreeMemory( proj4def );
788
// By default or if we have neither GDAL nor proj.4, we can't do squat
789
return std::string();
793
void SpatialReference::SetProj4(std::string const& v)
803
const char* poProj4 = v.c_str();
805
OGRSpatialReference srs(NULL);
806
if (OGRERR_NONE != srs.importFromProj4(const_cast<char *>(poProj4)))
808
throw std::invalid_argument("could not import proj4 into OSRSpatialReference SetProj4");
811
srs.exportToWkt(&poWKT);
813
std::string tmp(poWKT);
817
ret = GTIFSetFromOGISDefn( m_gtiff, tmp.c_str() );
820
throw std::invalid_argument("could not set m_gtiff from Proj4");
823
ret = GTIFWriteKeys(m_gtiff);
826
throw std::runtime_error("The geotiff keys could not be written");
831
if (m_gtiff && GTIFGetDefn(m_gtiff, &defn))
833
char* proj4def = GTIFGetProj4Defn(&defn);
834
std::string tmp(proj4def);
835
GTIFFreeMemory( proj4def );
838
boost::ignore_unused_variable_warning(v);
841
// if we have libgeotiff but not GDAL, we'll use the
842
// simple method in libgeotiff
843
#if defined(HAVE_LIBGEOTIFF) && !defined(HAVE_GDAL)
846
ret = GTIFSetFromProj4( m_gtiff, v.c_str());
849
throw std::invalid_argument("PROJ.4 string is invalid or unsupported");
852
ret = GTIFWriteKeys(m_gtiff);
855
throw std::runtime_error("The geotiff keys could not be written");
862
liblas::property_tree::ptree SpatialReference::GetPTree( ) const
864
using liblas::property_tree::ptree;
867
#if defined(HAVE_GDAL)
868
srs.put("proj4", GetProj4());
869
srs.put("prettywkt", GetWKT(liblas::SpatialReference::eHorizontalOnly, true));
870
srs.put("wkt", GetWKT(liblas::SpatialReference::eHorizontalOnly, false));
871
srs.put("compoundwkt", GetWKT(eCompoundOK, false));
872
srs.put("prettycompoundwkt", GetWKT(eCompoundOK, true));
873
srs.put("gtiff", GetGTIFFText());
876
#if defined(HAVE_LIBGEOTIFF) && !defined(HAVE_GDAL)
878
std::string message("Reference defined, but GDAL is not available for WKT support");
879
srs.put("proj4", GetProj4());
880
srs.put("prettywkt", message);
881
srs.put("wkt", message);
882
srs.put("compoundwkt", message);
883
srs.put("prettycompoundwkt", message);
884
srs.put("gtiff", GetGTIFFText());
887
#if !defined(HAVE_LIBGEOTIFF) && !defined(HAVE_GDAL)
890
if (m_vlrs.size() > 0 && m_wkt.size() == 0)
892
message = "Reference defined with VLR keys, but GeoTIFF and GDAL support are not available to produce definition";
893
} else if (m_wkt.size() > 0)
895
message = "Reference defined with WKT, but GeoTIFF and GDAL support are not available to produce definition";
901
srs.put("proj4", message);
902
srs.put("prettywkt", message);
903
srs.put("wkt", message);
904
srs.put("compoundwkt", message);
905
srs.put("prettycompoundwkt", message);
906
srs.put("gtiff", message);
913
std::string SpatialReference::GetGTIFFText() const
915
#ifndef HAVE_LIBGEOTIFF
916
return std::string("");
919
if( m_gtiff == NULL )
920
return std::string("");
922
detail::geotiff_dir_printer geotiff_printer;
923
GTIFPrint(m_gtiff, detail::libLASGeoTIFFPrint, &geotiff_printer);
924
return geotiff_printer.output();
929
} // namespace liblas
931
std::ostream& operator<<(std::ostream& ostr, const liblas::SpatialReference& srs)
935
liblas::property_tree::ptree tree;
936
std::string name ("spatialreference");
937
tree.put_child(name, srs.GetPTree());
938
liblas::property_tree::write_xml(ostr, tree);
942
boost::ignore_unused_variable_warning(ostr);
943
boost::ignore_unused_variable_warning(srs);
944
throw std::runtime_error("SpatialReference io operator<< is not available without GDAL+libgeotiff support");
added
removed
src/spatialreference.cpp
