max-zilla
commented
6 years ago Current progres, ran into same error @ZongyangLi observed with LAS overflows trying to convert to lat/long. So instead I'm converting to UTM with an offset. This generates the TIF but it's coarse (1m resolution) and not in the correct spot yet. Changing -dx and -dy to 0.1 (or smaller) makes the pixels smaller but does not increase sampling.
My test data is 2017-08-10__00-01-19-561. metadata.json is in raw_data/scanner3DTop, plyfiles are in Level_1/scanner3DTop, there is also raster using old code in Level_1/laser3d_heightmap to compare to.
def generate_las_from_ply(inp, out, pco):
"""
:param inp: list of input PLY files or single file path
:param out: output LAS file
:param pco: point cloud origin from metadata e.g. md['sensor_variable_metadata']['point_cloud_origin_m']['east']
"""
if not isinstance(inp, list):
inp = [inp]
# Create concatenated list of vertices to generate one merged LAS file
first = True
for plyf in inp:
plydata = PlyData.read(plyf)
merged_x = plydata['vertex']['x'] if first else numpy.concatenate([merged_x, plydata['vertex']['x']])
merged_y = plydata['vertex']['y'] if first else numpy.concatenate([merged_y, plydata['vertex']['y']])
merged_z = plydata['vertex']['z'] if first else numpy.concatenate([merged_z, plydata['vertex']['z']])
first = False
# PCO: {'y': 48.458000000000006, 'x': 127.55299999999998, 'z': 2.8050000000000006}
# merged x y z
# [-550.38256836 -549.65881348 -540.06091309 ..., 569.70611572 570.31878662 570.88928223]
# [-4115.94775391 -4115.02734375 -4128.1796875 ..., -25900.09570312 -25900.234375 -25900.1796875 ]
# [-55.74633789 -54.41064453 -73.07043457 ..., -8.83422852 -9.02001953 -8.92834473]
# Convert scanner coords to UTM
utm_x, utm_y = scanalyzer_to_mac(
(merged_x * 0.01) + pco['x'],
(merged_y * 0.01) + pco['y']
)
utm_z = (merged_z * 0.01) + pco['z']
# utm x y z
# [ 409006.03125 409006. 409006.15625 ..., 409223.65625 409223.65625 409223.65625]
# [ 3660097. 3660097. 3660097.25 ..., 3660106.5 3660106.5 3660106.5 ]
# [ 2.24753666 2.26089358 2.07429576 ..., 2.71665788 2.71479988 2.7157166 ]
# Create header and populate with scale and UTM-12 georeference offset
w = laspy.base.Writer(out, 'w', laspy.header.Header())
# https://gis.stackexchange.com/questions/185131/converting-between-coordinate-systems-for-las-file-using-laspy
xmin = numpy.floor(numpy.min(utm_x))
ymin = numpy.floor(numpy.min(utm_y))
zmin = numpy.floor(numpy.min(utm_z))
xmax = numpy.ceil(numpy.max(utm_x))
ymax = numpy.ceil(numpy.max(utm_y))
zmax = numpy.ceil(numpy.max(utm_z))
xrg = xmax - xmin
yrg = ymax - ymin
zrg = zmax - zmin
w.set_header_property("x_offset", xmin)
w.set_header_property("y_offset", ymin)
w.set_header_property("z_offset", zmin)
safety_factor = 1
maxval = 2e9
w.header.scale = [safety_factor*xrg/maxval, safety_factor*yrg/maxval, safety_factor*zrg/maxval]
print(utm_x)
print(utm_y)
print(utm_z)
print(w.header.offset)
print(w.header.scale)
w.set_x(utm_x, True)
w.set_y(utm_y, True)
w.set_z(utm_z, True)
print(w.get_points())
w.header.update_min_max(True)
w.close()
return
def generate_tif_from_las(inp, out, bounds, mode='max'):
"""
Create a raster (e.g. Digital Surface Map) from LAS pointcloud.
:param inp: input LAS file
:param out: output GeoTIFF file
:param bounds: (lat min, lat max, long min, long max) (lowerY, upperY, rightX, leftX)
:param mode: max | min | mean | median | sum (http://www.nongnu.org/pktools/html/md_pklas2img.html)
generally:
max = highest pixels in cell, usually canopy - equates to a DSM (Digital Surface Map)
min = lowest pixel in a cell, usually soil - equates to DTM (Digital Terrain Map)
"""
print('pklas2img -i '+inp+' -o '+out+' -comp '+mode+' -n z -a_srs epsg:32612')
subprocess.call(['pklas2img -i '+inp+' -o '+out+' -comp '+mode+' -n z -a_srs epsg:32612'], shell=True)
# Pretend we parsed this from the metadata
with open(md) as jmd:
md = clean_metadata(json.load(jmd), 'scanner3DTop')
pcoe = md['sensor_variable_metadata']['point_cloud_origin_m']['east']
pcow = md['sensor_variable_metadata']['point_cloud_origin_m']['west']
bbxe = md['spatial_metadata']['east']['bounding_box']
bbxw = md['spatial_metadata']['west']['bounding_box']
bounds_e = geojson_to_tuples(bbxe)
bounds_w = geojson_to_tuples(bbxw)
bounds_merge = (
bounds_e[0] if bounds_e[0] < bounds_w[0] else bounds_w[0],
bounds_e[1] if bounds_e[1] > bounds_w[1] else bounds_w[1],
bounds_e[2] if bounds_e[2] < bounds_w[2] else bounds_w[2],
bounds_e[3] if bounds_e[3] > bounds_w[3] else bounds_w[3],
)
if CONVERT:
print("east single")
generate_las_from_ply(inpe, eastout, pcoe)
generate_tif_from_las(eastout, eastout.replace(".las", ".tif"), bounds_e)
print("west single")
generate_las_from_ply(inpw, westout, pcow)
generate_tif_from_las(westout, westout.replace(".las", ".tif"), bounds_w)
print("merge!")
generate_las_from_ply([inpe, inpw], out, pcoe)
generate_tif_from_las(out, out.replace(".las", ".tif"), bounds_merge)
...until I divided the PCO Y offset by 2. I think the offset referred to the other end of the swath when we actually wanted the center? Either way I applied that change in desperation and now things line up perfectly. Here's the two 'west' scans:
And the two 'east' scans:
(we could almost get away with just using every other scan because there's a lot of overlap)
...and if you set it to, e.g. 0.1m resolution using dx/dy args for pktools you get this:
smaller pixels but not MORE pixels. Bizarre.
dlebauer
https://github.com/terraref/laser3d/pull/2
In the new implementation using plyfile + laspy + pktools, the LAS file headers include offset and scale that are supposed to georeference the LAS. However, pklas2img creates a geotiff that is not correctly georeferenced.
Currently in that PR, my code uses same transformation we used for Solmaz's BMP files to rotate and georeference the TIF file to a new TIF file. But if we can properly fix LAS header and/or pklas2img conversion call, we should be able to directly generate a georeferenced GeoTIFF from LAS file and not have to worry about other georeferencing.
Ideas:
If we can solve this issue for this sprint, we should be pretty much ready to deploy updated extractors for ply2las and heightmap.
Will post a bit more in the morning with sample files and images.