Thursday, March 03, 2011

Transformations in PostGIS

Lately I've been doing a lot of development centered around PostGIS (http://postgis.refractions.net/) and ran across a problem that has a relatively easy solution (for those well versed in GIS) but finding an answer proved to be somewhat difficult.

Here are the assumptions and the requirements before we look at the problem:

  • A point is defined by a latitude and longitude.
  • A line is a defined by a series of points.
  • A buffer is defined as a shape that envelopes a line. The width of the buffer is specified by the user in meters. The distance value the user supplies is a +/- distance from the supplied line.
  • We are assuming WGS 84 as the target map projection.

PostGIS provides a handy buffer function called ST_Buffer and in order to use this function the geometry supplied to the buffer function along with the desired buffer distance need to be in the same coordinate system. See the problem yet? My coordinates are in decimal degrees (i.e., latitude and longitude) and my buffer distance is in meters. So I need to transform the line into a different coordinate system that will allow me represent my geography based shape in metric coordinates. This transformation is non-trivial since converting meters to decimal degrees relies on the use of a map projection because the world isn't flat (sorry!). Before we go any further let's look at the query I've developed so far:

SELECT ST_Buffer(GeomFromText('LINESTRING(-76.543 42.567, -76.012 42.345, -75.890 42.445, -75.543 42.330)'), 300.0) AS buffer_shape

The query is valid but the buffer distance is thought to be in decimal degrees (i.e., WGS84) with respect to PostGIS and therefore produces a shape that pretty much covers the world. I tried guesstimating the conversion from decimal degrees to meters but quickly realized, at best, a guesstimate would be terribly wrong.

To solve this problem - I employed the Transform function from PostGIS. In order to get this function to work properly I had to determine what coordinate system to transform my geometry into...this is where my Google-fu fell relatively short. Luckily I went to http://gis.stackexchange.com and found a reference to the proj4js.org project. Using proj4js I was able to determine the SRID of the desired metric coordinate system (that SRID is 900913). This coordinate system bases coordinates off of meters instead of decimal degrees. So I changed my original query to be the following:

SELECT Transform(ST_Buffer(Transform(GeomFromText('LINESTRING(-76.543 42.567, -76.012 42.345, -75.890 42.445, -75.543 42.330)'), 900913), 300.0), 4326) AS buffer_shape

The above statement transforms the line from WGS 84 to a Mercator projection, performs the buffer operation, and then does another transformation to project the resulting buffer shape back into WGS 84. Voila! Again, nothing earth shattering here but if you don't speak GIS all day then performing the transformations may not be entirely obvious. In my opinion the PostGIS documentation falls short when mentioning the transformation details.

No comments: