This tutorial is partialy based on Yuriy’s Czoli article ‘Processing LiDAR to extract building heights’.
On Mac OSX
- Install Postgres, PostGIS and OSM-PostGIS tools
brew install postgis brew install protobuf-c brew install osm2pgsql --with-protobuf-c
Prepair the data
Download the OSM Metro Extracts data in
OSM XMLformat and the style file
Translate the LIDar data to
las2SM.py data.las SMdata.las
If you have many
.las files you can do
for file in ????????.las; do python las2SM.py $file SM$file ; done
And port it to XYZ format
las2txt -i SMdata.las —parse xyz —delimiter “ “ -o SMdata.xyz
for file in SM????????.las; do las2txt -i $file --parse xyz --delimiter " " -o $file.xyz ; done
Build and feed
postgres -D /usr/local/var/postgres
Create a new database.
Enter to the new database
Activate the postgis extenstion (for spatial objects and functionality).
CREATE EXTENSION postgis;
CREATE extension hstore;
Create a table to hold our LIDar dataset. This creates the infrastructure to hold our x,y,z data. We are naming the table ‘elevation’ because that is primarily what we are interested with in regards to the lidar data.
CREATE TABLE elevation (x double precision, y double precision, z double precision);
Now we will add a geometry column. PostGIS documentation suggests to add a geometry column after table creation. If you create the geometry column in the initial table set up, you will have to manually register the geometry in the database. Lets avoid this and follow the documentation. ‘the_geom’ refers to the name of the column, and 32610 is the SRID. This specific SRID is UTM zone 10N. UTM zones are great for dealing with a small area and measurements, which we will be doing. Check out this esri blog for a short summary.
SELECT AddGeometryColumn ('elevation','the_geom',3857,'POINT',2);
Copy the data from our text file to our database. NOTE: if this does not work, re-type the single quotes inside of the terminal.
copy elevation(x,y,z) FROM '~/your/path/SMdata.xyz' DELIMITERS ' ';
Lets update the elevation table with the ‘geometry from text’ function. This will parse through the table we just copied and pull the coordinates from the x y values. This will locate the points on the planar dimension of our data. This is using the table we just copied in from the previous step. Again using the utm 10n SRID. This will take a couple of minutes.
UPDATE elevation SET the_geom = ST_GeomFromText('POINT(' || x || ' ' || y || ')',3857);
Now that we have our spatial data loaded, we definitely want to create a spatial index. Spatial indices create general envelopes around our geometries, and will allow for quicker query and function abilities later on. Read more about them here.
CREATE INDEX elev_gix ON elevation USING GIST (the_geom);
We will switch our focus over to our shapefile for a second. Let’s exit our database.
Finnally load the OSM data to the database:
osm2pgsql -s -E 3857 -d testDataBase --hstore -S osm2pgsql.style.txt data-file.osm.bz2