Fwd: [Liblas-devel] Re: GetProj4() string and NAD83

[Adam Stylinski]

---------- Forwarded message ----------
From: Adam Stylinski <stylinae at mail.uc.edu>
Date: Mon, Jan 9, 2012 at 9:26 AM
Subject: Re: [Liblas-devel] Re: GetProj4() string and NAD83
To: Howard Butler <hobu.inc at gmail.com>


On Wed, Jan 4, 2012 at 2:18 PM, Howard Butler <hobu.inc at gmail.com> wrote:

> Adam,
>
> Can you give us some code that demonstrates what you're doing?
>
> You can use ReprojectionTransform in liblas/transform.hpp to handle
> reprojecting the data via GDAL. It also has support for rescaling the data
> if you give it a new liblas::Header with appropriate scale/offset in the
> case where you are going from something like NAD83 with 0.01 scale to
> something like WGS84 where you might want 0.0000001 or something.
>
> Howard
>
> On Dec 29, 2011, at 12:07 PM, Adam Stylinski wrote:
>
> > On Thu, Dec 29, 2011 at 12:42 PM, Adam Stylinski <stylinae at mail.uc.edu>
> wrote:
> > Hello,
> >
> > I'm having issues with LADAR files which do not use the WGS84 coordinate
> system.  I am translating coordinates into longitude and latitude and
> testing with google earth.  Whenever a LAS file has NAD83 coordinate it
> seems to have a tendency to be completely wrong (most of which end up being
> somewhere in the ocean).  Is there something I'm doing that doesn't work?
> >
> > I use the GetProj4() to get the proj4 string.  I then do a
> pj_init_plus(proj4str) to get the represented coordinate system, and then
> do a pj_latlong_from_proj() on that coordinate system to get the lat/long
> coordinate system.  Then to translate to lat/long I do a pj_transform with
> the src being the original coordinate system, the destination being the
> lat/long coordinate system, the next two parameters being 1,1, and the x
> and y being from a point pulled from GetX(), GetY().  Am I doing this wrong?
> >
> > I should probably mention that I'm failing on this one (though any NAD83
> coordinate seems to fail):
> > http://liblas.org/samples/IowaDNR-CloudPeakSoft-1.0-UTM15N.las
> > The GetProj4() returns this:
> > GetProj4() = +proj=tmerc +lat_0=0 +lon_0=0 +k=1 +x_0=0 +y_0=0
> +datum=NAD83 +units=m
> +geoidgrids=g2003conus.gtx,g2003alaska.gtx,g2003h01.gtx,g2003p01.gtx
> +no_defs
> > And GetWKT():
> > GetWKT() = COMPD_CS["unknown",
> >     PROJCS["NAD83 / UTM zone 15N",
> >         GEOGCS["NAD83",
> >             DATUM["North_American_Datum_1983",
> >                 SPHEROID["GRS 1980",6378137,298.2572221010002,
> >                     AUTHORITY["EPSG","7019"]],
> >                 AUTHORITY["EPSG","6269"]],
> >             PRIMEM["Greenwich",0],
> >             UNIT["degree",0.0174532925199433],
> >             AUTHORITY["EPSG","4269"]],
> >         PROJECTION["Transverse_Mercator"],
> >         PARAMETER["latitude_of_origin",0],
> >         PARAMETER["central_meridian",0],
> >         PARAMETER["scale_factor",1],
> >         PARAMETER["false_easting",0],
> >         PARAMETER["false_northing",0],
> >         UNIT["metre",1,
> >             AUTHORITY["EPSG","9001"]],
> >         AUTHORITY["EPSG","26915"]],
> >     VERT_CS["NAVD88 height",
> >         VERT_DATUM["North American Vertical Datum 1988",2005,
> >             AUTHORITY["EPSG","5103"],
> >
> EXTENSION["PROJ4_GRIDS","g2003conus.gtx,g2003alaska.gtx,g2003h01.gtx,g2003p01.gtx"]],
> >         UNIT["metre",1,
> >             AUTHORITY["EPSG","9001"]],
> >         AXIS["Up",UP],
> >         AUTHORITY["EPSG","5703"]]]
> > _______________________________________________
> > Liblas-devel mailing list
> > Liblas-devel at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/liblas-devel
>
>
> A pretty simple test case is something like this:

liblas::ReaderFactory f;
std::string proj4str;
liblas::SpatialReference srs;
 liblas::Reader r = f.CreateWithStream(ifs);
liblas::Header const& header = r.GetHeader();
numPoints = header.GetPointRecordsCount();
srs = header.GetSRS();
proj4str = srs.GetProj4();

projPJ projection = pj_init_plus(proj4str.c_str());
projPJ latlong = pj_latlong_from_proj(projection);

double destx = (X from ladar file);
double desty =  (Y from ladar file);
double destz = (Z from ladar file);
pj_transform(projection,latlong,1,1,&destx,&desty,&destz);
destx *= RAD_TO_DEG;
desty *= RAD_TO_DEG;

Then I go through some correctly behaving functions to convert decimal
degrees to DDMMSS.  pj_transform's return value I'm checking and it's zero.
 Some of the NAD83 ladar files are giving me what looks like HUGE_VAL
which, according to the proj API documentation, means the points are
failing to transform (yet somehow still pj_transform returns zero?).

I am not a geologist so I'm sorry for my limited experience and knowledge
when it comes to world coordinate systems and these APIs.  I'm probably
doing something completely wrong, though obviously for only use cases that
aren't WGS84.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/liblas-devel/attachments/20120110/f4a2d15d/attachment.html