The disk that held the personal geodatabases of our contour datasets died a while back, but not before I loaded the contours into PostGIS and started serving them up. Our new intern is working on putting together some shapefiles layer groups in ArcGIS for map production, and asked for a missing one… . I don’t have it, but in principle can always extract it from PostGIS and serve it back as a shapefile.
So I performed an intersection:
CREATE TABLE public.contour_2_cut ASSELECT contour.elevation, (ST_Intersection(contour.the_geom), cutline) AS the_geom, gid FROM public.contour_2 AS contour, (ST_GeomFromText('POLYGON((2000000 860000, 2500000 860000, 2500000 460000, 2000000 460000,2000000 860000))', 3734)) AS cutline;
After adding a primary key:
ALTER TABLE public.contour_2_cut ADD PRIMARY KEY (gid);
I couldn’t display it in uDig, which complained that it was a geometry collection. Ah, of course. In the ESRI world (whence I hark) an intersection takes the form of the least of the geometries involved, so the intersection of a polyline and a polygon is a polyline. Not so in PostGIS. PostGIS supports geometry collections, so naturally the intersection includes both types of geometry. (Side note– the viewing in uDig didn’t fail until it got to the polygon– it first loaded all the contour polylines, then discovering the polygon, threw an error).
What to do? Well, I hoped pgsql2shp would bail me out:
pgsql2shp -f contours_block_08_01 -h cmac-srv-gis -u postgres CM public.contour_2_cut
Initializing... type 'GEOMETRYCOLLECTION' is not Supported at this time.
The DBF file will be created but not the shx or shp files.
You've found a bug! (pgsql2shp.c:2864)
But it doesn’t support GEOMETRYCOLLECTIONs either.
Time to dig deeper. So how do we tell what a geometry is?
SELECT ST_GeometryType(st_geom) FROM contour_2_cut;
and I see ST_LineString is the descriptor for our lines. So, let’s create a copy of the table with just linestrings (we could remove records from the existing table as well):
CREATE TABLE contour_2_cut_line AS
SELECT * FROM contour_2_cut WHERE
ST_GeometryType(the_geom) ='ST_LineString';
Smathermather's Weblog 
Hallo
In 1.5 a new function was introduced to solve this type of problems:
ST_CollectionExtract
It will pull out the geometry type you are asking for from a geometry collection.
Regards
Nicklas
Thanks Nicklas,
Boy, I can’t upgrade to 1.5 fast enough. Your additions alone make the effort worth it. Do you have a sense if there is a performance difference between ST_CollectionExtract and testing the geometries as I did above?
Hallo again
Sorry, I didn’t read your post properly.
First, I guess the syntax in the first query has got some typo. The parenthesis seems strange placed for ST_Intersection.
Then I was thinking about why you get a geometry collection. It could be a point and a linestring, but never a polygon. My guess is that your geoemtry collection is avtuelly a GEOMETRYCOLLECTION EMPTY.
In PostGIS there is only one emty geometry and that is geometry collection empty. You get that when there is no intersection beteween the linestrings and the polygon. You can test for example:
SELECT ST_AsText(ST_Intersection(‘POLYGON((1 1, 1 10, 10 10, 10 1,1 1))’::geometry, ‘LINESTRING(15 15, 20 20)’::geometry));
What you should do is just to avoid sending the cases where you have no intersection to the ST_Intersection function.
So, you can use:
CREATE TABLE public.contour_2_cut AS
SELECT contour.elevation, ST_Intersection(contour.the_geom, cutline) AS the_geom, gid FROM
public.contour_2 AS contour, (ST_GeomFromText(‘POLYGON((2000000 860000, 2500000 860000, 2500000 460000, 2000000 460000,2000000 860000))’, 3734)) AS cutline
WHERE ST_Intersects(contour.the_geom, cutline);
But you can still get different types of output because, as mentioned above the linestring might stop at the boundary of the cutline polygon which will give you a point.
Then you can use ST_Overlaps instead of ST_Intersects to sort away those that stops on the boundary. But ST_Overlaps will not give you those linestring that is completely within your box, so you should combine with st_Within.
I think the working query will look like:
CREATE TABLE public.contour_2_cut AS
SELECT contour.elevation, ST_Intersection(contour.the_geom, cutline) AS the_geom, gid FROM
public.contour_2 AS contour, (ST_GeomFromText(‘POLYGON((2000000 860000, 2500000 860000, 2500000 460000, 2000000 460000,2000000 860000))’, 3734)) AS cutline
WHERE ST_Overlaps(contour.the_geom, cutline) OR ST_Within(contour.the_geom, cutline);
What you will miss here is if you have a linestring exactly following the boundary of your polygon, because that linestring will give a linestring from ST_Intersection, but not meet overlaps or within. But in your case I guess that is most unlikely to have a contour following your cutline boundary.
BTW, CollectionExtract is Pauls creation 🙂
Regards
Nicklas
I could easily have a typo in the query– I closed my session before realizing I wanted to post, so I was hoping I could retype without re-running… . Caught.
As to GEOMETRYCOLLECTION EMPTY, that’s likely, but I could have sworn I saw ST_LineString and ST_Polygon (or maybe ST_Multipolygon, I can’t remember). You’re right that logically it shouldn’t. I’ve already dropped the tables, but it will be interesting to rerun with your modification, so I’ll take another run at this, and retry my own query to verify the output from that. Besides, I can get GeoServer to display the different elements of it with an SLD, which will be really interesting to look at, and good practice for building GEOMETRYCOLLECTIONs as composite layers for GeoServer to render.
I wouldn’t have thought of the combo of ST_Overlaps and ST_Within as a replacement for ST_Intersects. It’s actually a likely case that there are linestrings starting at the edge of the polygon and emanating from there, and therefore would result in intersections resulting in a point. Thanks! New tool in my kit.