skeenp / QGIS3-getWKT

QGIS plugin that permits you to see the WKT of a geometry
GNU General Public License v3.0
6 stars 2 forks source link

WKT for 2.5D geometries #5

Closed heidivanparys closed 3 years ago

heidivanparys commented 4 years ago

When getting the WKT for polgyon with z values, the format does not following the standard (see section 7.2.3): there should be a space between "polygon" and "z".

E.g. when getting the WKT for the following GML layer:

<ogr:FeatureCollection
     gml:id="aFeatureCollection"
     xmlns:ogr="http://ogr.maptools.org/"
     xmlns:gml="http://www.opengis.net/gml/3.2">
    <ogr:featureMember>
        <ogr:export gml:id="id.1.geom">
            <gml:boundedBy>
                <gml:Envelope srsDimension="3" srsName="http://www.opengis.net/def/crs/EPSG/0/7416">
                    <gml:lowerCorner>722120.48 6178866.39 24.28</gml:lowerCorner>
                    <gml:upperCorner>722184.83 6178905.66 25.01</gml:upperCorner>
                </gml:Envelope>
            </gml:boundedBy>
            <ogr:geometryProperty>
                <gml:Polygon srsName="http://www.opengis.net/def/crs/EPSG/0/7416" gml:id="id.1.geom">
                    <gml:exterior>
                        <gml:LinearRing>
                            <gml:posList srsDimension="3">722178.52 6178866.39 24.61 722184.83 6178880.97 24.61 722137.52 6178901.05 24.7 722135.5 6178901.91 25.01 722126.66 6178905.66 24.61 722120.48 6178891.19 24.61 722126.3 6178888.7 24.28 722178.52 6178866.39 24.61</gml:posList>
                        </gml:LinearRing>
                    </gml:exterior>
                </gml:Polygon>
            </ogr:geometryProperty>
        </ogr:export>
    </ogr:featureMember>
</ogr:FeatureCollection>

The WKT retrieved with the plugin is PolygonZ ((722178.52000000001862645 6178866.38999999966472387 24.60999999999999943, 722184.82999999995809048 6178880.96999999973922968 24.60999999999999943, 722137.52000000001862645 6178901.04999999981373549 24.69999999999999929, 722135.5 6178901.91000000014901161 25.01000000000000156, 722126.66000000003259629 6178905.66000000014901161 24.60999999999999943, 722120.47999999998137355 6178891.19000000040978193 24.60999999999999943, 722126.30000000004656613 6178888.70000000018626451 24.28000000000000114, 722178.52000000001862645 6178866.38999999966472387 24.60999999999999943))

It should start with Polygon Z, not PolygonZ.

There is an additional issue: the plugin also returns too many decimals, these are not specified in the GML. But that may be a QGIS issue?

When using ogrinfo:

ogrinfo export.gml -wkt_format WKT2 -al -q

the result is

Layer name: export
OGRFeature(export):0
  POLYGON Z ((722178.52 6178866.39 24.61,722184.83 6178880.97 24.61,722137.52 6178901.05 24.7,722135.5 6178901.91 25.01,722126.66 6178905.66 24.61,722120.48 6178891.19 24.61,722126.3 6178888.7 24.28,722178.52 6178866.39 24.61))

And when using ogr2ogr:

ogr2ogr -f CSV export.csv export.gml -lco GEOMETRY=AS_WKT

the result in the CSV file is also

POLYGON Z ((722178.52 6178866.39 24.61,722184.83 6178880.97 24.61,722137.52 6178901.05 24.7,722135.5 6178901.91 25.01,722126.66 6178905.66 24.61,722120.48 6178891.19 24.61,722126.3 6178888.7 24.28,722178.52 6178866.39 24.61))

So ogr writes the WKT according to the standard, including a space between polygon and z.

IBM does it according to the standard too, see https://www.ibm.com/support/knowledgecenter/SSHRBY/com.ibm.swg.im.dashdb.spatial.doc/doc/rsbp4120.html?view=embed

skeenp commented 4 years ago

The plugin uses QGIS's inbuild asWkt() functionality to retrieve WKT representations of geometries.

there should be a space between "polygon" and "z".

Hmm, yes and no. This is where we could get into an argument of who's standard is better

QGIS exports the WKT as polygonz in line with a number of standards. The standard in use probably dates back to the program's roots in PostGIS compatibility which uses one of the standards defining no space. ISO 13249-3, Information technology - Database languages - SQL Multimedia and Application Packages - Part 3: Spatial Simple feature access - Part 2: SQL option

But you are right, other standards exist that define the space and if its a problem its a problem.

There is an additional issue: the plugin also returns too many decimals, these are not specified in the GML. But that may be a QGIS issue?"""

You are correct, the extra decimals derive from how QGIS stores coordinates internally. Every coordinate system has different expectations for the number of decimal places in a coordinate (think Lat Long vs UTM). I believe QGIS uses Double type with a set precision regardless, as it probably doesn't want to take a guess and risk losing data. I don't know well enough to explain why the coordinates change very slightly, and its an absolute pain that I have banged my head against a wall about for years.

The Vertex Editor in QGIS seems to default to 4 decimal places for all coordinate systems. This is good for the UTMs, providing micrometre accuracy but terrible for LatLong up to 11 meters.

So what can I do?

I'd be happy to add in the option to configure the number of decimal places exported and specify the format for WKT type specifiers. It probably would not be the default, but could be put in a settings menu.

We could advance that and put some checks on the output coordinate system to determine what 'should' be the precision. Projected Coordinate Systems could be, for example, 4 for meters and 3 for feet (0.1-0.3 millimetre accuracy), Geographic Systems could be raw. The output decimal precision could be set in settings as auto, custom or raw.

I've stayed away from adding a settings menu to try to keep things simple, but if it's causing grief I'll be happy to get to work on it. Let me know what you think.

skeenp commented 4 years ago

Interestingly, PostGIS via OGR doesn't even specify its a 2.5 or 3D geometry in its type. Its got the data in the string, but doesn't explicitly specify z. How does you DB systems handle this?

ogr2ogr -f PGDump dump.sql export.gml -lco WRITE_EWKT_GEOM=ON

POLYGON ((722178.52 6178866.39 24.61,722184.83 6178880.97 24.61,722137.52 6178901.05 24.7,722135.5 6178901.91 25.01,722126.66 6178905.66 24.61,722120.48 6178891.19 24.61,722126.3 6178888.7 24.28,722178.52 6178866.39 24.61))

Could we just have an option to remove the z?

Interestingly in PostGIS all of the following return POINT Z (722178.52 6178866.39 24.61)

select st_astext(st_geomfromtext('POINTZ (722178.52 6178866.39 24.61)')); select st_astext(st_geomfromtext('POINT Z (722178.52 6178866.39 24.61)')); select st_astext(st_geomfromtext('POINT (722178.52 6178866.39 24.61)'));

heidivanparys commented 4 years ago

QGIS exports the WKT as polygonz in line with a number of standards. The standard in use probably dates back to the program's roots in PostGIS compatibility which uses one of the standards defining no space. ISO 13249-3, Information technology - Database languages - SQL Multimedia and Application Packages - Part 3: Spatial Simple feature access - Part 2: SQL option

But you are right, other standards exist that define the space and if its a problem its a problem.

The standard I referred to is the OpenGIS Implementation Standard for Geographic Information Simple feature access - Part 1: Common architecture , version 1.2.1, so we are both referring to the same standards (SQL/MM is based on the Simple feature access, and WKT is actually defined in Simple feature access part 1, not part 2).

From Simple feature access part 1:

image

From the SQL/MM document you referred to above:

image

image

I don't yet have access to the latest version of ISO 13249-3, but I would be surprised if the spaces have been removed in there.

Interestingly, PostGIS via OGR doesn't even specify its a 2.5 or 3D geometry in its type. Its got the data in the string, but doesn't explicitly specify z. How does you DB systems handle this?

ogr2ogr -f PGDump dump.sql export.gml -lco WRITE_EWKT_GEOM=ON

POLYGON ((722178.52 6178866.39 24.61,722184.83 6178880.97 24.61,722137.52 6178901.05 24.7,722135.5 6178901.91 25.01,722126.66 6178905.66 24.61,722120.48 6178891.19 24.61,722126.3 6178888.7 24.28,722178.52 6178866.39 24.61))

But that is EWKT (PostGIS specific format, see https://postgis.net/docs/ST_GeomFromEWKT.html, not WKT (OGC standard).

Could we just have an option to remove the z?

That would be against the standard, which seems to be supported, also for 3D, in quite some database systems). I guess the reason it must be there, is that else, it is not possible to differentiate between e.g. Polygon Z and Polygon M.

Interestingly in PostGIS all of the following return POINT Z (722178.52 6178866.39 24.61)

select st_astext(st_geomfromtext('POINTZ (722178.52 6178866.39 24.61)')); select st_astext(st_geomfromtext('POINT Z (722178.52 6178866.39 24.61)')); select st_astext(st_geomfromtext('POINT (722178.52 6178866.39 24.61)'));

I guess PostGIS applies the Robustness principle : "Be conservative in what you send, be liberal in what you accept".

I tried now with a Spatialite database in QGIS:

image

Some additional testing with Spatialite from within QGIS:

SELECT ST_AsText(ST_GeomFromText("POLYGON Z ((722178.52 6178866.39 24.61,722184.83 6178880.97 24.61,722137.52 6178901.05 24.7,722135.5 6178901.91 25.01,722126.66 6178905.66 24.61,722120.48 6178891.19 24.61,722126.3 6178888.7 24.28,722178.52 6178866.39 24.61))")) as wkt FROM "test"

and

SELECT ST_AsText(ST_GeomFromText("POLYGONZ ((722178.52 6178866.39 24.61,722184.83 6178880.97 24.61,722137.52 6178901.05 24.7,722135.5 6178901.91 25.01,722126.66 6178905.66 24.61,722120.48 6178891.19 24.61,722126.3 6178888.7 24.28,722178.52 6178866.39 24.61))")) as wkt FROM "test"

give

POLYGON Z((722178.52 6178866.39 24.61, 722184.83 6178880.97 24.61, 722137.52 6178901.05 24.7, 722135.5 6178901.91 25.01, 722126.66 6178905.66 24.61, 722120.48 6178891.19 24.61, 722126.3 6178888.7 24.28, 722178.52 6178866.39 24.61))

whereas

SELECT ST_AsText(ST_GeomFromText("POLYGON ((722178.52 6178866.39 24.61,722184.83 6178880.97 24.61,722137.52 6178901.05 24.7,722135.5 6178901.91 25.01,722126.66 6178905.66 24.61,722120.48 6178891.19 24.61,722126.3 6178888.7 24.28,722178.52 6178866.39 24.61))")) as wkt FROM "test"

gives

NULL

So the robustness principle seems to apply there as well, POLYGONZ and POLYGON Z are accepted as input, and the output is always POLYGON Z.

So what can I do?

I'd be happy to add in the option to configure the number of decimal places exported and specify the format for WKT type specifiers. It probably would not be the default, but could be put in a settings menu.

We could advance that and put some checks on the output coordinate system to determine what 'should' be the precision. Projected Coordinate Systems could be, for example, 4 for meters and 3 for feet (0.1-0.3 millimetre accuracy), Geographic Systems could be raw. The output decimal precision could be set in settings as auto, custom or raw.

I've stayed away from adding a settings menu to try to keep things simple, but if it's causing grief I'll be happy to get to work on it. Let me know what you think.

As for the format for the WKT, I would definitely suggest to stick to the standard (which, ideally, the QGIS source code should do).

As for the number of decimal places exported, the solution you propose sounds fine to me.

skeenp commented 4 years ago

Oops, there is the other problem with standards, me misunderstanding them 🙁.

I'll put in a modifier to change polygonz to polygon z, as well as the same for zm and m values. I agree a fix to qgis core would be better but I should be able to implement this alot quicker. Have you/did you want to raise this at qgis core as an issue? I'm happy to if you would like.

I will tackle the decimal places as well.

That's for your feedback and the time taken to explain the issues. Really appreciated.

heidivanparys commented 4 years ago

Have you/did you want to raise this at qgis core as an issue? I'm happy to if you would like.

I haven't raised this at QGIS core. If you could do that, that would be much appreciated, thanks!

skeenp commented 4 years ago

Added code changes to the repo, will do some more tests over the next couple of days. Will release over at qgis when tests are complete.

skeenp commented 3 years ago

Closed with submittion to repo on March 26, 2020