Postgis - How do i check the geometry type before i do an insert

i have a postgres database with millions of rows in it it has a column called geom which contains the boundary of a property.

using a python script i am extracting the information from this table and re-inserting it into a new table.

when i insert in the new table the script bugs out with the following:

Traceback (most recent call last):
  File "", line 258, in <module>
  File "", line 166, in main
    update_cursor.executemany("insert into parcels (par_id, street_add, title_no, proprietors, au_name, ua_name, geom) VALUES (%s, %s, %s, %s, %s, %s, %s)", inserts)
psycopg2.IntegrityError: new row for relation "parcels" violates check constraint "enforce_geotype_geom"

The new table has a check constraint enforce_geotype_geom = ((geometrytype(geom) = 'POLYGON'::text) OR (geom IS NULL)) whereas the old table does not, so im guessing theres dud data or non polygon (perhaps multipolygon data?) in the old table. i want to keep the new data as polygon so want to not insert anything else.

Initially i tried wrapping the query with standard python error handling with the hope that the dud geom rows would fail but the script would keep running , but the script has been written to commit at the end not each row so it doesnt work.

I think what i need to do is iterate through the old table geom rows and check what type of geometry they are so i can establish whether or not i want to keep it or throw it away before i insert into the new table

Whats the best way of going about this?

Asked by: Miranda481 | Posted: 06-12-2021

Answer 1

This astonishingly useful bit of PostGIS SQL should help you figure it out... there are many geometry type tests in here:

-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-- $Id: cleanGeometry.sql 2008-04-24 10:30Z Dr. Horst Duester $
-- cleanGeometry - remove self- and ring-selfintersections from 
--                 input Polygon geometries 
-- Copyright 2008 SO!GIS Koordination, Kanton Solothurn, Switzerland
-- Version 1.0
-- contact: horst dot duester at bd dot so dot ch
-- This is free software; you can redistribute and/or modify it under
-- the terms of the GNU General Public Licence. See the COPYING file.
-- This software is without any warrenty and you use it at your own risk
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

CREATE OR REPLACE FUNCTION cleanGeometry(geometry)
  RETURNS geometry AS
  inGeom ALIAS for $1;
  outGeom geometry;
  tmpLinestring geometry;


  outGeom := NULL;

-- Clean Process for Polygon 
  IF (GeometryType(inGeom) = 'POLYGON' OR GeometryType(inGeom) = 'MULTIPOLYGON') THEN

-- Only process if geometry is not valid, 
-- otherwise put out without change
    if not isValid(inGeom) THEN

-- create nodes at all self-intersecting lines by union the polygon boundaries
-- with the startingpoint of the boundary.  
      tmpLinestring := st_union(st_multi(st_boundary(inGeom)),st_pointn(boundary(inGeom),1));
      outGeom = buildarea(tmpLinestring);      
      IF (GeometryType(inGeom) = 'MULTIPOLYGON') THEN      
        RETURN st_multi(outGeom);
        RETURN outGeom;
      END IF;
      RETURN inGeom;
    END IF;

-- Clean Process for LINESTRINGS, self-intersecting parts of linestrings 
-- will be divided into multiparts of the mentioned linestring 
  ELSIF (GeometryType(inGeom) = 'LINESTRING') THEN

-- create nodes at all self-intersecting lines by union the linestrings
-- with the startingpoint of the linestring.  
    outGeom := st_union(st_multi(inGeom),st_pointn(inGeom,1));
    RETURN outGeom;
  ELSIF (GeometryType(inGeom) = 'MULTILINESTRING') THEN 
    outGeom := multi(st_union(st_multi(inGeom),st_pointn(inGeom,1)));
    RETURN outGeom;
  ELSIF (GeometryType(inGeom) = '<NULL>' OR GeometryType(inGeom) = 'GEOMETRYCOLLECTION') THEN 
    RAISE NOTICE 'The input type % is not supported %',GeometryType(inGeom),st_summary(inGeom);
    RETURN inGeom;
  END IF;     

Answered by: Jack413 | Posted: 07-01-2022

Answer 2

Option 1 is to create a savepoint before each insert and roll back to that safepoint if an INSERT fails.

Option 2 is to attach the check constraint expression as a WHERE condition on the original query that produced the data to avoid selecting it at all.

The best answer depends on the size of the tables, the relative number of faulty rows, and how fast and often this is supposed to run.

Answered by: Chester713 | Posted: 07-01-2022

Answer 3

I think you can use ST_CollectionExtract — Given a (multi)geometry, returns a (multi)geometry consisting only of elements of the specified type.

I use it when inserting the results of an ST_Intersection, ST_Dump breaks any multi-polygon, collections into individual geometry. Then ST_CollectionExtract (theGeom, 3) discards anything but polygons:

ST_CollectionExtract((st_dump(ST_Intersection(data.polygon, grid.polygon))).geom, )::geometry(polygon, 4326)

The second parameter above 3 can be: 1 == POINT, 2 == LINESTRING, 3 == POLYGON

Answered by: Maria669 | Posted: 07-01-2022

Similar questions

python - How do I overlap widgets with the Tkinter pack geometry manager?

I want to put a Canvas with an image in my window, and then I want to pack widgets on top of it, so the Canvas acts as a background. Is it possible to have two states for the pack manager: one for one set of widgets and another for another set?

Good geometry library in python?

Closed. This question does not meet Stack Overflow guid...

python - How do I use Django to insert a Geometry Field into the database?

class LocationLog(models.Model): user = models.ForeignKey(User) utm = models.GeometryField(spatial_index=True) This is my database model. I would like to insert a row. I want to insert a circle at point -55, 333. With a radius of 10. How can I put this circle into the geometry field? Of course, then I would want to check which circles overlap a given circle. (my select stat...

Help with Windows Geometry in Python

Why are the commands to change the window position before and after sleep(3.00) being ignored? if self.selectedM.get() == 'Bump': W1 = GetSystemMetrics(1) + 200 print W1 w1.wm_geometry("+100+" + str(W1)) w2.wm_geometry("+100+" + str(W1)) w3.wm_geometry("+100+" + str(W1)) w4.wm_geometry("+100+" + str(W1)) self.rvar.set(0) self.rvar2.set(0) ...

python - Parsing MSDN Geometry Data Type

I have a database where one field gives spatial coordinates. I have learned the field is a serialised MSDN geometry Data Type ( I want to access this database from Python and was wandering if anyone knew the format of the Geometry Data Type, or any libraries capable of parsing ...

python - GUI layout using Tk Grid Geometry Manager

Building a small application for personal use with Python and thought I'd try my hand with a little GUI programming using Tkinter. This is the GUI I've created so far: Application doubts: How can I make sure that the three LableFrames - A, B and C in the screenshot - have the same width? (Or rather, have ...

geometry - Integral of a triangle function in python

I have just defined a function which is the integral of a given step function, and now I would like to integrate the triangle function obtained. jerk:1,-1,1-1 The step function correspond of a given list of successive jerk. For getting the triangle function, I have used the formula: a(t)=kt+ac with k the Jerk and ac a con...

python - Is there a GUI design app for the Tkinter / grid geometry?

Closed. This question does not meet Stack Overflow guid...

geometry - Python - Pygame Drawing Angle

I need to draw an angle(in Pygame), given it's 1. Measure of the angle (θ) 2. Endpoints of the base(A &amp; B) Here I know 1. Measure of θ (in radians and degrees) 2. (x,y) of A and B 3. Measure of BC My Question How do I calculate the position of the Co-ordinates(...

python - Convert Point Geometry to list

I have the following script that creates a point geometry. How can I convert this point geometry to a list only containing the coordiantes to look like [258432.79138201929, 1001957.4394514663]? &gt;&gt;&gt; import ogr &gt;&gt;&gt; driver = ogr.GetDriverByName('ESRI Shapefile') &gt;&gt;&gt; pointshp = driver.Open('U:/My Documents/Tool/shp/point.shp', 0) &gt;&gt;&gt; pointlyr = pointshp.GetLayer...

Still can't find your answer? Check out these communities...

PySlackers | Full Stack Python | NHS Python | Pythonist Cafe | Hacker Earth | Discord Python