Yes, there is. its called "Convert format". Go to processing toolbox and run it from there by right click then run as batch processing.
UPDATE :
put as follow in creation option column:
-dsco AltitudeMode=absolute -dsco NameField=Name -dsco DescriptionField=Description
For example, in my case will be
-dsco AltitudeMode=clampToGround -dsco NameField=NAME_1 -dsco DescriptionField=Description
TEST RESULT
Be Careful, when insert the creation option attribute (No space at the end allowed, seems it affecting the result):
Wrong → 
Right → 
This extra description value can be modified as you need like AltitudeMode can be clampToGround or relativeToGround like the one that you see within "Save as" method.
UPDATE 2 - Alternative
-dsco altitudeModeEnum=clampToGround -dsco NameField=Name -dsco DescriptionField=description
Alternative option has been tested.
Source : GDAL.org and KML Guide for Google Developer
You can interpolate the two axes independently like:
Code
def find_xy(p1, p2, z):
x1, y1, z1 = p1
x2, y2, z2 = p2
if z2 < z1:
return find_xy(p2, p1, z)
x = np.interp(z, (z1, z2), (x1, x2))
y = np.interp(z, (z1, z2), (y1, y2))
return x, y
Test Code:
import numpy as np
# arbitrary z
z1, x1, y1 = 0.822, 174269.8332, 608470.8717
z2, x2, y2 = 0.692, 174269.5191, 608480.8668
print(find_xy((x1, y1, z1), (x2, y2, z2), .7))
print(find_xy((x1, y1, z1), (x2, y2, z2), .8))
# z = 0
z1, x1, y1 = 0.022, 174267.9486, 608530.8421
z2, x2, y2 = -0.262, 174267.6345, 608540.8372
print(find_xy((x1, y1, z1), (x2, y2, z2), 0))
Results:
(174269.53842923077, 608480.2517169231)
(174269.78004461539, 608472.5631784615)
(174267.92426830987, 608531.6163683098)
Best Answer
Your data seems to be equally spaced, so I suggest to load it as a XYZ Raster file. You have to obey the expected sorting order: http://www.gdal.org/frmt_xyz.html
If qgis fails on the file size, try to do the conversion with pure GDAL. It gives better error messages, and saves the valuable RAM that the QGIS GUI needs.