Loading OSM and LIDar to PostGIS

postgisOSM-LAS.md

This tutorial is partialy based on Yuriy’s Czoli article ‘Processing LiDAR to extract building heights’.

Install PostGIS

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

1. Download the OSM Metro Extracts data in OSM XML format and the style file osm2pgsql.sty

2. Translate the LIDar data to epsg:3857 projection using las2SM.py

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

or

for file in SM????????.las; do las2txt -i $file --parse xyz --delimiter " " -o $file.xyz ; done

Build and feed

Start postgres.

postgres -D /usr/local/var/postgres  

Create a new database.

createdb testDataBase

Enter to the new database

psql testDataBase

Activate the postgis extenstion (for spatial objects and functionality).

CREATE EXTENSION postgis;

Enable hstore

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.

\q

Finnally load the OSM data to the database:

osm2pgsql -s -E 3857 -d testDataBase --hstore -S osm2pgsql.style.txt data-file.osm.bz2