1

I'm trying to aggregate data from a tabular data source and then join it to a geometry file (sl_adm2) in PostGIS using a primary key. The geometry file was imported using shp2psql.

CREATE TABLE example_geom AS 
SELECT sl_adm2_sums_alias.* , adm2_code as district_pcode 
FROM gis_schema.sl_adm2
JOIN (
    SELECT st_collect(sl_adm2.geom) AS geometry, sum(sl_sos.HEALTH) AS health_sum,
        sum(sl_sos.PROTECTION) AS protection_sum, sum(sl_sos.WASH) AS wash_sum,
        sum(sl_sos.cccm) AS cccm_sum, sum(sl_sos.food) AS food_sum, 
        sum(sl_sos.livelihood) AS livelihood_sum, 
        sum(sl_sos.education) AS education_sum, 
        sum(sl_sos.agriculture) AS agriculture_sum, 
        sum(sl_sos.HUMANITARIAN_ACCESS) AS humaccess_sum, 
        sum(sl_sos.logistics) AS logistics_sum, 
        sum(sl_sos.SHELTER_AND_NFI) AS shelter_and_wifi_sum, sl_sos.district
    FROM sl_sos 
    JOIN gis_schema.sl_adm2 ON sl_sos.DISTRICT_PCODE = sl_adm2.adm2_code
    GROUP BY sl_sos.district)
sl_adm2_sums_alias ON
sl_adm2.adm2_name = sl_adm2_sums_alias.district;

and then

ALTER TABLE example_geom
ADD PRIMARY KEY (district_pcode);

The aggregation and join work fine but the geometry column in the resulting table does not appear to work. Further investigation reveals that the geometry column is registering as a GeometryCollection. I tried to change this with the following code.

ALTER TABLE example_geom ALTER COLUMN geometry
SET DATA TYPE geometry(MultiPolygon) USING ST_Multi(geometry);

This gave me the below error

ERROR: Geometry type (GeometryCollection) does not match column type (MultiPolygon)
SQL state: 22023

I've tried to fix this or get around it with no success. Could anyone help?

I'm using Postgres 9.4.1 with PostGIS version 2.1.5 on a Mac OS X 10.10.2.

  • Try wrapping your st_collect in a st_multi, eg ... st_multi(st_collect(sl_adm2.geom)) .... When you say it does not appear to work, do you mean you can't view it? Geometry collections are valid, but some software won't deal with them. – MickyT Mar 1 '15 at 18:03
  • Are you happy to get rid of everything that is not a multipolygon. There might be some points or lines also. – John Powell Mar 2 '15 at 14:35
  • @JohnBarça This is a purely polygonal dataset so I have no interest in points or lines. – rbanick Mar 2 '15 at 18:18
  • @MickyT Using st_multi(st_collect(sl_adm2.geom)), the original SQL processes fine but the final command returns the same error as above. I cannot view the geometry in QGIS or Mapbox Studio, no. It may be that there are other issues but I can't troubleshoot until I have the file in polygon form. – rbanick Mar 2 '15 at 18:21
2

I'm guessing that you may have a few invalid geometries in your base table. You will need to do a bit of investigation.

Another thing to do would be select the well known text (ST_AsText) from your created table then look for anomalies, eg POINTs and LINEs in the collections. You might want to add a few examples to your question.

Firstly I would make sure that your base geometries are valid

UPDATE gis_schema.sl_adm2
SET geom = ST_MakeValid(geom)
WHERE ST_IsValid(geom) IS FALSE

Then when you build your table add a check in to make sure only polygons are brought back. You may want to consider using ST_Union rather than ST_Collect, but the depends on the result you are trying to achieve. I recall that MULTIPOLYGONs can't having touching or overlapping parts, but I can't find a reference to it at the moment. Even though it is slower ST_UNION will likely give you better results. Be aware that it will dissolve boundaries between touching polygons.

CREATE TABLE example_geom AS 
SELECT sl_adm2_sums_alias.* , adm2_code as district_pcode 
FROM gis_schema.sl_adm2
JOIN (
    SELECT ST_Multi(st_collect(sl_adm2.geom)) AS geometry, sum(sl_sos.HEALTH) AS health_sum,
        sum(sl_sos.PROTECTION) AS protection_sum, sum(sl_sos.WASH) AS wash_sum,
        sum(sl_sos.cccm) AS cccm_sum, sum(sl_sos.food) AS food_sum, 
        sum(sl_sos.livelihood) AS livelihood_sum, 
        sum(sl_sos.education) AS education_sum, 
        sum(sl_sos.agriculture) AS agriculture_sum, 
        sum(sl_sos.HUMANITARIAN_ACCESS) AS humaccess_sum, 
        sum(sl_sos.logistics) AS logistics_sum, 
        sum(sl_sos.SHELTER_AND_NFI) AS shelter_and_wifi_sum, sl_sos.district
    FROM sl_sos 
    JOIN gis_schema.sl_adm2 ON sl_sos.DISTRICT_PCODE = sl_adm2.adm2_code
    WHERE ST_GeometryType(geom) in ('ST_Polygon','ST_MultiPolygon')
    GROUP BY sl_sos.district)
sl_adm2_sums_alias ON
sl_adm2.adm2_name = sl_adm2_sums_alias.district;
  • This works for me, but ONLY if I use ST_Union. Obviously I'm very happy to have this working but still curious what could cause that. Odder still the source polygon (sl_adm2) loads just fine and shows up as MultiPolygon in various GIS GUIs. So it's not a problem with the base geometries as far as I can tell. – rbanick Mar 3 '15 at 12:11
  • @rbanick I recall that multipolygons can't have adjoining or overlapping parts. I will check it out and change my answer when I get by a computer – MickyT Mar 3 '15 at 16:35
  • I'm using this polygon layer: ebolageonode.org/layers/geonode%3Asle_admn_ad2_py_ocha_district – rbanick Mar 4 '15 at 9:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.