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 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;
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';