PostGIS – Parsing Coordinates from st_isvalidreason() to a Geometry Column

postgispostgis-topology

With PostGIS, I'm using st_isvalidreason() function to check topologicals errors in my entities. When something is topologically wrong, it returns me something like:

Self-intersection[969880.222423769 6274973.28067924]

First I have the reason and after, surrounding by brackets I have coordinates of the error. You can notice that this brackets aren't a real ARRAY, because it's full text and there is no delimiter.

I wan't to use this coordinates to make a geometry column and display it in a GIS (like QGIS) like a point layer.

But how to parse this coordinates information to build a geometry column ?

Ideally I'd like it to look like st_point(969880.222423769,6274973.28067924).

Best Answer

The ST_IsValidDetail https://postgis.net/docs/ST_IsValidDetail.html function returns the location of the error as geometry. It can be converted into WKT with ST_AsText if needed.

select ST_IsValidReason(ST_GeomFromText('POLYGON (( 340 440, 460 440, 460 320, 400 320, 440 380, 360 380, 400 320, 340 320, 340 440 ))'))
as reason, 
ST_AsText(location(ST_IsValidDetail(
ST_GeomFromText('POLYGON (( 340 440, 460 440, 460 320, 400 320, 440 380, 360 380, 400 320, 340 320, 340 440 ))')
))) as place;

Result:
"Ring Self-intersection[400 320]"   "POINT(400 320)"

Related Solutions

[GIS] Adding Geometry column in PostGIS

First, as others have said, make sure you have PostGIS installed. Run this command to ensure that you have PostGIS installed,

SELECT postgis_version();

If you don't have postgis_version() then you need to install PostGIS on the database. If it's not installed, you can try this when connected to the specific database you want to install PostGIS on,

CREATE EXTENSION postgis;

If that fails than you need to install or build PostGIS which then becomes a question of what operating system you're using. If you've got the extension downloaded and installed, run this in psql

\dfS AddGeometryColumn

If that works, then you have the function installed too but you're providing it the wrong types. We need the actual function call. Provide us what you're calling AddGeometryColumn() with. Also, if your version of PostGIS is newer than 2.0, do not use AddGeometryColumn(). It is not needed. Instead use the ALTER TABLE ADD COLUMN with the subtype. Here is an example to add a geometry point in srs 4326.

ALTER TABLE some_table ADD COLUMN geom geometry(Point,4326);

PostGIS Geometry – Parsing Geometry WKB with OGR in PostGIS

Internally, PostGIS stores geometries in a binary specification, but it is queried and viewed outside as a hex-encoded string. There are two popular variations of well-known binary (WKB):

EWKB (via ST_AsEWKB) - an extended WKB specification designed by PostGIS.
OGC WKB (via ST_AsBinary) - specified by the OGC and ISO. For a while it was 2D-only, but later extended to support Z, M and ZM geometries.

The two specifications are the same for 2D geometries, but are different for higher-order geometries with Z, M and ZM coordinates.

Older versions of GDAL/OGR (1.x) only understand the EWKB for 3D geometries, so for these I recommend using ST_AsEWKB. (But if you only have 2D geometries, either format are fine). For example:

import psycopg2
from osgeo import ogr

ogr.UseExceptions()    
conn = psycopg2.connect('dbname=postgis user=postgres')
curs = conn.cursor()

curs.execute("select ST_AsEWKB('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex'))  # 0101000080000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToWkt())  # POINT (1 2 3)

curs.execute("select ST_AsBinary('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex'))  # 01e9030000000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
# RuntimeError: OGR Error: Unsupported geometry type

Also, note that older GDAL/OGR versions do not support M coordinates, and these will be parsed but ignored.

With GDAL 2.0 and more recent, ISO WKT/WKB is supported. This means that CreateGeometryFromWkb can read either WKB flavour (without specifying) and ExportToIsoWkt() shows output with a modern WKT syntax.

import psycopg2
from osgeo import ogr

ogr.UseExceptions()
conn = psycopg2.connect('dbname=postgis user=postgres')
curs = conn.cursor()

curs.execute("select ST_AsEWKB('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex'))  # 0101000080000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToIsoWkt())  # POINT Z (1 2 3)

curs.execute("select ST_AsBinary('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex'))  # 01e9030000000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToIsoWkt())  # POINT Z (1 2 3)

Additionally, GDAL 2.1 or later will create/export WKT/WKB with M or ZM coordinates as expected.

Best Answer

Related Solutions

[GIS] Adding Geometry column in PostGIS

PostGIS Geometry – Parsing Geometry WKB with OGR in PostGIS

Related Question