As in my previous post, you can see how a generic approach to taking irregularly spaced text format data may be applied to a wide range of contexts. In this post I show a particular example of taking some freely available LIDAR result data and processing it with the same approach.
Get the Data
To follow along, you can grab the same data I did from OpenTopography.org. You can find your way through their search tools for an area of interest.
My resulting data is in CSV/ASCII format, with no gridding applied. Grab the dataset and metadata from here.
Unzipping the file produces a 138MB text file called points.txt. This is already in a format that OGR can read, but to help it along just rename the file to points.csv and we don’t have to mess around changing the delimiter (as in previous post). Note the first line contains column headings too, so this is a perfect example to be using:
x,y,z,classification,gpstime,scan_angle,intensity,number_of_returns...
330000.05,7001998.62,1339.69,2,0.000000,0,171,0,0,3,0,0,0,0,0,0
330000.13,7001999.41,1340.28,1,0.000000,0,115,0,0,3,0,0,0,0,0,0
330000.32,7001999.60,1340.12,2,0.000000,0,131,0,0,3,0,0,0,0,0,0
330000.22,7001998.78,1339.83,2,0.000000,0,184,0,0,3,0,0,0,0,0,0
Create VRT file
In case it wasn’t clear previously, we create a VRT file so that gdal_grid (or any OGR application) knows how to find the X/Y coordinates in the CSV file. So we just create and tweak like in the previous example but include the field names that are already in points.txt, etc. Here is our new points.vrt file:
<OGRVRTDataSource>
<OGRVRTLayer name="points">
<SrcDataSource>points.csv</SrcDataSource>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>EPSG:32607</LayerSRS>
<GeometryField encoding="PointFromColumns" x="x" y="y" z="z"/>
</OGRVRTLayer>
</OGRVRTDataSource>
Test results by running against ogrinfo:
ogrinfo points.vrt -al -so
INFO: Open of `points.vrt'
using driver `VRT' successful.
Layer name: points
Geometry: Point
Feature Count: 2267774
Extent: (330000.000000, 7001270.110000) - (330589.740000, 7001999.990000)
Layer SRS WKT:
PROJCS["WGS 84 / UTM zone 7N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-141],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
AUTHORITY["EPSG","32607"],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]
x: String (0.0)
y: String (0.0)
z: String (0.0)
classification: String (0.0)
gpstime: String (0.0)
scan_angle: String (0.0)
intensity: String (0.0)
number_of_returns: String (0.0)
return_number: String (0.0)
point_source_ID: String (0.0)
edge_of_flight_line_flag: String (0.0)
direction_of_scan_flag: String (0.0)
user_data: String (0.0)
red_channel: String (0.0)
green_channel: String (0.0)
blue_channel: String (0.0)
Grid the Data
Here’s the part that takes a while. With > 2 million rows it takes a bit of time to process, depending on what algorithm settings you use. To just use the default (inverse distance to a power) run it without an algorithm option:
gdal_grid -zfield z -l points points.vrt points.tif
Grid data type is "Float64"
Grid size = (256 256).
Corner coordinates = (329998.848164 7002001.415547)-(330590.891836 7001268.684453).
Grid cell size = (2.303672 2.851094).
Source point count = 2267774.
Algorithm name: "invdist".
Options are "power=2.000000:smoothing=0.000000:radius1=0.000000:
radius2=0.000000:angle=0.000000:max_points=0:min_points=0:nodata=0.000000"
I also got good results from the using the nearest neighbour options as shown in the gdal_grid docs. Note it creates a 256×256 raster by default. In the second example it shows how to set the output size (-outsize 512 512). The second example took about an hour to process on my one triple core AMD 64 desktop machine with 8GB RAM:
gdal_grid -zfield z -outsize 512 512 \
-a nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0 \
-l points points.vrt points_nearest_lrg.tif
Grid data type is "Float64"
Grid size = (512 512).
Corner coordinates = (329999.424082 7002000.702773)-(330590.315918 7001269.397227).
Grid cell size = (1.151836 1.425547).
Source point count = 2267774.
Algorithm name: "nearest".
Options are "radius1=0.000000:radius2=0.000000:angle=0.000000:nodata=0.000000"
