[Liblas-commits] hg: 2 new changesets
liblas-commits at liblas.org
liblas-commits at liblas.org
Tue Jan 12 00:22:46 EST 2010
changeset 098d89056bcc in /Volumes/Data/www/liblas.org/hg
details: http://hg.liblas.orghg?cmd=changeset;node=098d89056bcc
summary: added vertical coordinate system support
changeset 6123def8e21c in /Volumes/Data/www/liblas.org/hg
details: http://hg.liblas.orghg?cmd=changeset;node=6123def8e21c
summary: added python bindings for vertical cs stuff, and very modest tests
diffstat:
apps/las2las.c | 64 ++++++++-
apps/lascommon.c | 2 +-
include/liblas/capi/liblas.h | 6 +-
include/liblas/lasspatialreference.hpp | 31 ++++-
python/liblas/core.py | 9 +
python/liblas/srs.py | 11 +-
python/tests/SRS-GDAL.txt | 12 +-
src/gt_wkt_srs.cpp | 234 +++++++++++++++++++++++++++++++-
src/las_c_api.cpp | 32 ++++-
src/lasspatialreference.cpp | 64 ++++++++-
10 files changed, 445 insertions(+), 20 deletions(-)
diffs (truncated from 711 to 300 lines):
diff -r d4b4cdaefa95 -r 6123def8e21c apps/las2las.c
--- a/apps/las2las.c Mon Jan 11 13:46:53 2010 -0500
+++ b/apps/las2las.c Tue Jan 12 00:17:16 2010 -0500
@@ -140,6 +140,8 @@
LASPointH surviving_point_max = NULL;
double surviving_gps_time_min;
double surviving_gps_time_max;
+ int verticalCSType = -1, verticalDatum = -1, verticalUnits = 9001;
+ const char *verticalCitation = "";
int clipped = 0;
int eliminated_return = 0;
@@ -336,6 +338,44 @@
do_reprojection = TRUE;
}
}
+ else if ( strcmp(argv[i],"--a_srs") == 0 ||
+ strcmp(argv[i],"-a_srs") == 0
+ )
+ {
+ ++i;
+ if (LAS_IsGDALEnabled()) {
+ out_srs = LASSRS_Create();
+ ret = LASSRS_SetFromUserInput(out_srs, argv[i]);
+ if (ret) {
+ LASError_Print("Unable to import output SRS");
+ exit(1);
+ }
+ }
+ }
+ else if ( strcmp(argv[i],"--a_vertcs") == 0 ||
+ strcmp(argv[i],"-a_vertcs") == 0
+ )
+ {
+ ++i;
+ verticalCSType = atoi(argv[i]);
+ ++i;
+ if( i < argc && argv[i][0] != '-' )
+ {
+ verticalCitation = argv[i];
+ ++i;
+
+ if( i < argc && argv[i][0] != '-' )
+ {
+ verticalDatum = atoi(argv[i]);
+ ++i;
+ if( i < argc && argv[i][0] != '-' )
+ {
+ verticalUnits = atoi(argv[i]);
+ ++i;
+ }
+ }
+ }
+ }
else if ( strcmp(argv[i],"--scale") == 0 ||
strcmp(argv[i],"-scale") == 0
)
@@ -811,6 +851,22 @@
LASHeader_SetDataOffset(surviving_header, LASHeader_GetDataOffset(surviving_header)+abs(header_pad));
}
+ /* Do we have vertical cs info to set? */
+ if( verticalCSType > 0 )
+ {
+ if( out_srs == NULL )
+ out_srs = LASHeader_GetSRS(surviving_header);
+
+ if( out_srs == NULL )
+ out_srs = LASSRS_Create();
+
+ LASSRS_SetVerticalCS( out_srs,
+ verticalCSType,
+ verticalCitation,
+ verticalDatum,
+ verticalUnits );
+ }
+
if (do_reprojection) {
if (verbose) {
proj4_text = LASSRS_GetProj4(out_srs);
@@ -825,8 +881,14 @@
}
LASHeader_SetSRS(surviving_header, out_srs);
-
}
+
+ /* Are we just assigning an override SRS? (-a_srs) */
+ else if( out_srs != NULL )
+ {
+ LASHeader_SetSRS(surviving_header, out_srs);
+ }
+
if (verbose) {
fprintf(stderr,
"second pass reading %d and writing %d points ...\n",
diff -r d4b4cdaefa95 -r 6123def8e21c apps/lascommon.c
--- a/apps/lascommon.c Mon Jan 11 13:46:53 2010 -0500
+++ b/apps/lascommon.c Tue Jan 12 00:17:16 2010 -0500
@@ -435,7 +435,7 @@
pSRS = LASHeader_GetSRS(header);
pszProj4 = LASSRS_GetProj4(pSRS);
- pszWKT = LASSRS_GetWKT(pSRS);
+ pszWKT = LASSRS_GetWKT_CompoundOK(pSRS);
pGTIF = LASSRS_GetGTIF(pSRS);
nVLR = LASHeader_GetRecordsCount(header);
diff -r d4b4cdaefa95 -r 6123def8e21c include/liblas/capi/liblas.h
--- a/include/liblas/capi/liblas.h Mon Jan 11 13:46:53 2010 -0500
+++ b/include/liblas/capi/liblas.h Tue Jan 12 00:17:16 2010 -0500
@@ -1117,11 +1117,15 @@
LAS_DLL const GTIF* LASSRS_GetGTIF(LASSRSH hSRS);
-LAS_DLL char* LASSRS_GetWKT(LASSRSH hSRS);
+LAS_DLL char* LASSRS_GetWKT(LASSRSH hSRS );
+LAS_DLL char* LASSRS_GetWKT_CompoundOK( LASSRSH hSRS );
LAS_DLL LASError LASSRS_SetWKT(LASSRSH hSRS, const char* value);
LAS_DLL LASError LASSRS_SetFromUserInput(LASSRSH hSRS, const char* value);
LAS_DLL char* LASSRS_GetProj4(LASSRSH hSRS);
LAS_DLL LASError LASSRS_SetProj4(LASSRSH hSRS, const char* value);
+LAS_DLL LASError LASSRS_SetVerticalCS(LASSRSH hSRS, int verticalCSType,
+ const char *citation, int verticalDatum,
+ int verticalUnits );
LAS_DLL LASSRSH LASHeader_GetSRS(const LASHeaderH hHeader);
LAS_DLL LASError LASHeader_SetSRS(LASHeaderH hHeader, const LASSRSH hSRS);
LAS_DLL void LASSRS_Destroy(LASSRSH hSRS);
diff -r d4b4cdaefa95 -r 6123def8e21c include/liblas/lasspatialreference.hpp
--- a/include/liblas/lasspatialreference.hpp Mon Jan 11 13:46:53 2010 -0500
+++ b/include/liblas/lasspatialreference.hpp Tue Jan 12 00:17:16 2010 -0500
@@ -90,6 +90,11 @@
class LASSpatialReference
{
public:
+ enum WKTModeFlag
+ {
+ eHorizontalOnly = 1,
+ eCompoundOK = 2
+ };
/// Default constructor.
LASSpatialReference();
@@ -116,13 +121,35 @@
/// Returns the OGC WKT describing Spatial Reference System.
/// If GDAL is linked, it uses GDAL's operations and methods to determine
/// the WKT. If GDAL is not linked, no WKT is returned.
- std::string GetWKT() const;
-
+ /// \param mode_flag May be eHorizontalOnly indicating the WKT will not
+ /// include vertical coordinate system info (the default), or
+ /// eCompoundOK indicating the the returned WKT may be a compound
+ /// coordinate system if there is vertical coordinate system info
+ /// available.
+ std::string GetWKT( WKTModeFlag mode_flag = eHorizontalOnly ) const;
+
/// Sets the SRS using GDAL's OGC WKT. If GDAL is not linked, this
/// operation has no effect.
/// \param v - a string containing the WKT string.
void SetWKT(std::string const& v);
+ /// Sets the vertical coordinate system using geotiff key values.
+ /// This operation should normally be done after setting the horizontal
+ /// portion of the coordinate system with something like SetWKT(),
+ /// SetProj4(), SetGTIF() or SetFromUserInput()
+ /// \param verticalCSType - An EPSG vertical coordinate system code,
+ /// normally in the range 5600 to 5799, or -1 if one is not available.
+ /// \param citation - a textual description of the vertical coordinate
+ /// system or an empty string if nothing is available.
+ /// \param verticalDatum - the EPSG vertical datum code, often in the
+ /// range 5100 to 5299 - implied by verticalCSType if that is provided, or
+ /// -1 if no value is available.
+ /// \param verticalUnits - the EPSG vertical units code, often 9001 for Metre.
+ void SetVerticalCS(int verticalCSType,
+ std::string const& citation = "",
+ int verticalDatum = -1,
+ int verticalUnits = 9001 );
+
/// Sets the SRS using GDAL's SetFromUserInput function. If GDAL is not linked, this
/// operation has no effect.
/// \param v - a string containing the definition (filename, proj4, wkt, etc).
diff -r d4b4cdaefa95 -r 6123def8e21c python/liblas/core.py
--- a/python/liblas/core.py Mon Jan 11 13:46:53 2010 -0500
+++ b/python/liblas/core.py Tue Jan 12 00:17:16 2010 -0500
@@ -628,6 +628,15 @@
las.LASSRS_SetProj4.argtypes = [ctypes.c_void_p, ctypes.c_char_p]
las.LASSRS_SetProj4.errcheck = check_return
+las.LASSRS_SetVerticalCS.restype = ctypes.c_int
+las.LASSRS_SetVerticalCS.argtypes = [ctypes.c_void_p, ctypes.c_int,
+ ctypes.c_char_p, ctypes.c_int,
+ ctypes.c_int]
+las.LASSRS_SetVerticalCS.errcheck = check_return
+
+las.LASSRS_GetWKT_CompoundOK.argtypes = [ctypes.c_void_p]
+las.LASSRS_GetWKT_CompoundOK.errcheck = check_value_free
+las.LASSRS_GetWKT_CompoundOK.restype = ctypes.c_char_p
las.LASSRS_GetWKT.argtypes = [ctypes.c_void_p]
las.LASSRS_GetWKT.errcheck = check_value_free
las.LASSRS_GetWKT.restype = ctypes.c_char_p
diff -r d4b4cdaefa95 -r 6123def8e21c python/liblas/srs.py
--- a/python/liblas/srs.py Mon Jan 11 13:46:53 2010 -0500
+++ b/python/liblas/srs.py Tue Jan 12 00:17:16 2010 -0500
@@ -55,9 +55,18 @@
if self.owned:
if self.handle and core:
core.las.LASSRS_Destroy(self.handle)
-
+
+ def set_verticalcs( self, verticalCSType, citation = '',
+ verticalDatum = -1,
+ verticalUnits = 9001 ):
+ return core.las.LASSRS_SetVerticalCS( self.handle, verticalCSType,
+ citation, verticalDatum,
+ verticalUnits )
+
def get_wkt(self):
return core.las.LASSRS_GetWKT(self.handle)
+ def get_wkt_compoundok(self):
+ return core.las.LASSRS_GetWKT_CompoundOK(self.handle)
def set_wkt(self, value):
return core.las.LASSRS_SetWKT(self.handle, value)
wkt = property(get_wkt, set_wkt)
diff -r d4b4cdaefa95 -r 6123def8e21c python/tests/SRS-GDAL.txt
--- a/python/tests/SRS-GDAL.txt Mon Jan 11 13:46:53 2010 -0500
+++ b/python/tests/SRS-GDAL.txt Tue Jan 12 00:17:16 2010 -0500
@@ -115,4 +115,14 @@
>>> os.remove('junk_srs_project.las')
-
\ No newline at end of file
+ >>> f = file.File('../test/data/srs_vertcs.las',mode='r')
+ >>> s = f.header.srs
+ >>> s.get_wkt_compoundok() == """COMPD_CS["unknown",PROJCS["WGS 84 / UTM zone 17N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32617"]],VERT_CS["North American Vertical Datum of 1988",VERT_DATUM["North American Vertical Datum 1988",2005,AUTHORITY["EPSG","5103"]],UNIT["metre",1.0,AUTHORITY["EPSG","9001"]],AXIS["Up",UP],AUTHORITY["EPSG","5703"]]]"""
+ True
+
+ >>> s2 = srs.SRS()
+ >>> s2.wkt = """PROJCS["WGS 84 / UTM zone 17N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32617"]]"""
+ >>> s2.set_verticalcs( 5703, 'abc', 5103, 9001 )
+ True
+ >>> s2.get_wkt_compoundok()
+ 'COMPD_CS["unknown",PROJCS["WGS 84 / UTM zone 17N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32617"]],VERT_CS["abc",VERT_DATUM["North American Vertical Datum 1988",2005,AUTHORITY["EPSG","5103"]],UNIT["metre",1.0,AUTHORITY["EPSG","9001"]],AXIS["Up",UP],AUTHORITY["EPSG","5703"]]]'
diff -r d4b4cdaefa95 -r 6123def8e21c src/gt_wkt_srs.cpp
--- a/src/gt_wkt_srs.cpp Mon Jan 11 13:46:53 2010 -0500
+++ b/src/gt_wkt_srs.cpp Tue Jan 12 00:17:16 2010 -0500
@@ -1,5 +1,5 @@
/******************************************************************************
- * $Id: gt_wkt_srs.cpp 16723 2009-04-06 16:15:17Z hobu $
+ * $Id: gt_wkt_srs.cpp 18490 2010-01-09 05:44:49Z warmerdam $
*
* Project: GeoTIFF Driver
* Purpose: Implements translation between GeoTIFF normalized projection
@@ -42,7 +42,6 @@
#include "xtiffio.h"
#include "cpl_multiproc.h"
-
CPL_C_START
void GTiffOneTimeInit();
int CPL_DLL VSIFCloseL( FILE * );
@@ -122,6 +121,8 @@
pszDatum = CPLStrdup(*ppszDatum);
GTIFFreeMemory( *ppszDatum );
*ppszDatum = pszDatum;
+ if (pszDatum[0] == '\0')
+ return;
/* -------------------------------------------------------------------- */
/* Translate non-alphanumeric values to underscores. */
@@ -177,11 +178,11 @@
/************************************************************************/
/* For example:
- GTCitationGeoKey (Ascii,215): "IMAGINE GeoTIFF Support\nCopyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 16723 $ $Date: 2009-04-06 11:15:17 -0500 (Mon, 06 Apr 2009) $\nProjection Name = UTM\nUnits = meters\nGeoTIFF Units = meters"
+ GTCitationGeoKey (Ascii,215): "IMAGINE GeoTIFF Support\nCopyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 18490 $ $Date: 2010-01-09 00:44:49 -0500 (Sat, 09 Jan 2010) $\nProjection Name = UTM\nUnits = meters\nGeoTIFF Units = meters"
- GeogCitationGeoKey (Ascii,267): "IMAGINE GeoTIFF Support\nCopyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 16723 $ $Date: 2009-04-06 11:15:17 -0500 (Mon, 06 Apr 2009) $\nUnable to match Ellipsoid (Datum) to a GeographicTypeGeoKey value\nEllipsoid = Clarke 1866\nDatum = NAD27 (CONUS)"
+ GeogCitationGeoKey (Ascii,267): "IMAGINE GeoTIFF Support\nCopyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 18490 $ $Date: 2010-01-09 00:44:49 -0500 (Sat, 09 Jan 2010) $\nUnable to match Ellipsoid (Datum) to a GeographicTypeGeoKey value\nEllipsoid = Clarke 1866\nDatum = NAD27 (CONUS)"
- PCSCitationGeoKey (Ascii,214): "IMAGINE GeoTIFF Support\nCopyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 16723 $ $Date: 2009-04-06 11:15:17 -0500 (Mon, 06 Apr 2009) $\nUTM Zone 10N\nEllipsoid = Clarke 1866\nDatum = NAD27 (CONUS)"
+ PCSCitationGeoKey (Ascii,214): "IMAGINE GeoTIFF Support\nCopyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 18490 $ $Date: 2010-01-09 00:44:49 -0500 (Sat, 09 Jan 2010) $\nUTM Zone 10N\nEllipsoid = Clarke 1866\nDatum = NAD27 (CONUS)"
*/
@@ -229,6 +230,15 @@
OGRSpatialReference oSRS;
/* -------------------------------------------------------------------- */
+/* Handle non-standard coordinate systems where GTModelTypeGeoKey */
+/* is not defined, but ProjectedCSTypeGeoKey is defined (ticket #3019) */
+/* -------------------------------------------------------------------- */
+ if( psDefn->Model == KvUserDefined && psDefn->PCS != KvUserDefined)
+ {
+ psDefn->Model = ModelTypeProjected;
+ }
+
+/* -------------------------------------------------------------------- */
/* Handle non-standard coordinate systems as LOCAL_CS. */
/* -------------------------------------------------------------------- */
More information about the Liblas-commits
mailing list