SQL and ogr to create 4D lines

A friend recently asked me to help him generate polyline shapefiles with Z and M values that he could deliver to a customer.  The problem was, the software he was using supported the import of Z and M values, but did not support the export of those files.  The other problem was, he has zillions of data tables that he needed to export!

Fortunately, PostGIS allows us to create Polyline ZM data rather easily.  The next part was to figure out how to get it exported to a shapefile.  So, here goes:

Assume you have a table with 4 columns: lineID, Z, M, and geometry.  The geometry field is a point feature that represent vertices on a line.  The lineID separates all those points by which line they are part of, and the Z and M values are, well, Z and M.

Make the Polyline ZM

The SQL command to create a 4D line is:

SELECT  lineid, 
ST_MakeLine(ST_MakePoint(ST_X(geometry), 
                         ST_Y(geometry), z, m)
            ) AS g
FROM zpts
GROUP BY lineid;

In this case, I am making a line (ST_MakeLine) as a series of X (ST_X), Y (ST_Y), Z, and M values.  Grouping the results of the query by the lineid allows me to create a line for each set of points.

Export the Shapefile

I tried using the latest OSGeo4W version of GDAL to export the shapefile but kept getting an error that ogr2ogr did not support 4D geometry.  If I created a Polyline Z, no problem.  But, adding the measure gave me problems.  So, I inserted the line into a PostGRES table, and issued:

pgsql2shp -f c:\temp\data\newfile.shp -h localhost -P postgres 
          -u postgres postgres public.zwkb

in the above, I created a table named zwkb.  But, I was still not satisfied.  After looking over lots of bug reports from the GDAL sites, I realized the the OsGeo4w version was rather old, but a newer version of GDAL fixed the ogr2ogr command.  So, after downloading the newer GDAL libraries, I was able to issue the following:

ogr2ogr -f "ESRI Shapefile" c:\temp\newfile.shp PG:"host=localhost 
        user=postgres  dbname=postgres password=postgres" "zwkb"

Now I was able to write out my shapefile with an XYZM.  But, I was still creating a table inside of PostGRES.  A better option is to issue the SQL from ogr2ogr and do everything in one shot:

ogr2ogr -f "ESRI Shapefile" c:\temp\newfile.shp PG:"host=localhost 
user=postgres dbname=postgres password=postgres" 
-sql "SELECT lineid, ST_MakeLine(ST_MakePoint(ST_X(geometry), ST_Y(geometry), z, m)) AS g FROM zpts GROUP BY lineid"

I highlighted the SQL in the ogr2ogr command so you can see how the query was formed.

This now created a really effective way to create and export lots of Polyline ZM files to a shapefile.

Want to learn more about spatial SQL?  Check out my online course Spatial SQL: A Language for Geographers.