<br><br><div class="gmail_quote">---------- Forwarded message ----------<br>From: <b class="gmail_sendername">Adam Stylinski</b> <span dir="ltr"><<a href="mailto:stylinae@mail.uc.edu">stylinae@mail.uc.edu</a>></span><br>
Date: Mon, Jan 9, 2012 at 9:26 AM<br>Subject: Re: [Liblas-devel] Re: GetProj4() string and NAD83<br>To: Howard Butler <<a href="mailto:hobu.inc@gmail.com">hobu.inc@gmail.com</a>><br><br><br><div><div><div></div><div class="h5">
<div class="gmail_quote">On Wed, Jan 4, 2012 at 2:18 PM, Howard Butler <span dir="ltr"><<a href="mailto:hobu.inc@gmail.com" target="_blank">hobu.inc@gmail.com</a>></span> wrote:</div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Adam,<br>
<br>
Can you give us some code that demonstrates what you're doing?<br>
<br>
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.<br>
<br>
Howard<br>
<div><div></div><div><br>
On Dec 29, 2011, at 12:07 PM, Adam Stylinski wrote:<br>
<br>
> On Thu, Dec 29, 2011 at 12:42 PM, Adam Stylinski <<a href="mailto:stylinae@mail.uc.edu" target="_blank">stylinae@mail.uc.edu</a>> wrote:<br>
> Hello,<br>
><br>
> 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?<br>
><br>
> 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?<br>
><br>
> I should probably mention that I'm failing on this one (though any NAD83 coordinate seems to fail):<br>
> <a href="http://liblas.org/samples/IowaDNR-CloudPeakSoft-1.0-UTM15N.las" target="_blank">http://liblas.org/samples/IowaDNR-CloudPeakSoft-1.0-UTM15N.las</a><br>
> The GetProj4() returns this:<br>
> 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<br>
> And GetWKT():<br>
> GetWKT() = COMPD_CS["unknown",<br>
> PROJCS["NAD83 / UTM zone 15N",<br>
> GEOGCS["NAD83",<br>
> DATUM["North_American_Datum_1983",<br>
> SPHEROID["GRS 1980",6378137,298.2572221010002,<br>
> AUTHORITY["EPSG","7019"]],<br>
> AUTHORITY["EPSG","6269"]],<br>
> PRIMEM["Greenwich",0],<br>
> UNIT["degree",0.0174532925199433],<br>
> AUTHORITY["EPSG","4269"]],<br>
> PROJECTION["Transverse_Mercator"],<br>
> PARAMETER["latitude_of_origin",0],<br>
> PARAMETER["central_meridian",0],<br>
> PARAMETER["scale_factor",1],<br>
> PARAMETER["false_easting",0],<br>
> PARAMETER["false_northing",0],<br>
> UNIT["metre",1,<br>
> AUTHORITY["EPSG","9001"]],<br>
> AUTHORITY["EPSG","26915"]],<br>
> VERT_CS["NAVD88 height",<br>
> VERT_DATUM["North American Vertical Datum 1988",2005,<br>
> AUTHORITY["EPSG","5103"],<br>
> EXTENSION["PROJ4_GRIDS","g2003conus.gtx,g2003alaska.gtx,g2003h01.gtx,g2003p01.gtx"]],<br>
> UNIT["metre",1,<br>
> AUTHORITY["EPSG","9001"]],<br>
> AXIS["Up",UP],<br>
> AUTHORITY["EPSG","5703"]]]<br>
</div></div>> _______________________________________________<br>
> Liblas-devel mailing list<br>
> <a href="mailto:Liblas-devel@lists.osgeo.org" target="_blank">Liblas-devel@lists.osgeo.org</a><br>
> <a href="http://lists.osgeo.org/mailman/listinfo/liblas-devel" target="_blank">http://lists.osgeo.org/mailman/listinfo/liblas-devel</a><br>
<br>
<br>
</blockquote></div></div></div>A pretty simple test case is something like this:<div><br></div><div><div>liblas::ReaderFactory f;</div><div>std::string proj4str;</div><div>liblas::SpatialReference srs;</div><div><span style="white-space:pre-wrap"> </span></div>
<div>liblas::Reader r = f.CreateWithStream(ifs);</div><div>liblas::Header const& header = r.GetHeader();</div><div>numPoints = header.GetPointRecordsCount();</div><div>srs = header.GetSRS();</div><div>proj4str = srs.GetProj4();</div>
<br><div class="gmail_quote">projPJ projection = pj_init_plus(proj4str.c_str());</div><div class="gmail_quote">projPJ latlong = pj_latlong_from_proj(projection);</div><div class="gmail_quote"><br></div><div class="gmail_quote">
double destx = (X from ladar file);</div><div class="gmail_quote">double desty = (Y from ladar file);</div><div class="gmail_quote">double destz = (Z from ladar file);</div><div class="gmail_quote">pj_transform(projection,latlong,1,1,&destx,&desty,&destz);</div>
<div class="gmail_quote">destx *= RAD_TO_DEG;</div><div class="gmail_quote">desty *= RAD_TO_DEG;</div><div class="gmail_quote"><br></div><div class="gmail_quote">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?). </div>
</div></div><div class="gmail_quote"><br></div><div class="gmail_quote">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. </div>
</div><br>