As I mentioned in the previous post, there were a couple of handy PostGIS functions I found quite useful after importing prepped data via shp2pgsql tool. They involve working with Spatial Reference of geometries and working with Multi* geometries.
Even after figuring out the way to define Spatial Reference in QGIS – it’s still possible to miss providing that SRID information in the PostGIS Shapefile Import/Export Manager or shp2pgsql command line tool – I know one such person (hint, it’s me)… Oh well, let’s figure out how to do the same this time in PostGIS, so we don’t need to drop the table and re-import the shapefile. It indeed can be done quite easily via setting the SRID on the whole table’s geometry column:
ALTER TABLE mytable
ALTER COLUMN geom
TYPE Geometry(Polygon, 26916)
USING ST_SetSRID(geom, 26916);
Note, this will work if your geometry data is already properly projected and you just need to supply the projection name via EPSG code.
Also, note the difference with ST_Transform – in this case your existing spatial reference is already known and defined, but is not the one you need the data to be stored in; you need to actually change the way geometry is saved in the table. For example you may have multiple input shapefiles in different UTM zones that you’ve imported, data from which you want to store together in one PostGIS table in lat long.
ALTER TABLE mytable
ALTER COLUMN geom
TYPE Geometry(Polygon, 4326)
USING ST_Transform(geom, 4326);
Of course, these functions are very useful not just for modifying SRID of the whole table! They can be applied to separate features and record values as well. For example, in the above case, geometries from all different UTM zones could be transformed on the fly at the time of inserting records from separate tables to the combo table:
INSERT INTO combotable(id, geom)
SELECT gid, ST_Transform(geom, 4326)
Next, here is an example of using both functions when generating point geometries into ‘geom’ column from UTM zone coordinate values stored in the same table in columns ‘x’ and ‘y’. First a point is generated, then its spatial reference is defined as UTM zone 16, and finally it is reprojected to WGS84 and saved into the ‘geom’ column:
SET geom = ST_Transform(ST_SetSRID(ST_Point(x, y), 26916), 4326);
If you are not sure what SRID is defined on your table, if at all – another handy function ST_SRID will help you find out:
2. PostGIS: ST_Dump
Sometimes you may have situations where the data is stored via multi geometries, but you need to do operations on each single element of those objects without permanently splitting them. For example, I needed to know if simple polygon features of layer A covered at least some of the polygons making up multypolygon features of the layer B. ST_Dump turned out to be the perfect tool for that!
SELECT a.gid, b.gid, a.geom, b.geom
FROM simplepoly_table_a a, (select gid, (ST_Dump(geom)).geom from multipoly_table_b) b
WHERE st_within(b.geom, a.geom)
Note that ST_Dump returns a geometry_dump construct, containing geom and path to the feature as (geom, path). Thus, to access just the geometry, you need to call ST_Dump(geom) of the multi geometry, and then call the geometry of that as ST_Dump(geom)).geom. By using ST_Dump to build alias ‘b’ we are able to check each simple polygon geometry of multipoly_table_b for an overlay with each polygon of simplepoly_table_a.
This is it for this article, hope you found it useful. It is safe to say that PostGIS supplies a great wealth of diverse functionality for working with spatial data. Be sure to explore the reference if you want to learn more about performing similar tasks, or wonder if your unique need can be addressed!
Make sure to read Daily Discoveries of a FOSS4G User: QGIS, the Vector Doctor.